Quantification of Nematic Cell Polarity in Three-dimensional Tissues

TL;DR

This paper introduces a new framework for quantifying complex nematic cell polarity patterns in 3D tissues, revealing local coordination and alignment with tissue structures, and models liver tissue as a biaxial liquid crystal.

Contribution

It develops a novel algorithmic approach using multipole expansion and co-orientational order parameters to analyze nematic cell polarity in 3D tissues.

Findings

Nematic cell polarity axes show local coordination in liver tissue.

Polarity axes align with the sinusoidal blood transport network.

Liver tissue can be modeled as a biaxial liquid crystal.

Abstract

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-Navarette et al., eLife 8:e44860, 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 or in-vitro organoids.

Tables

Table 1. Table S1: List of symbols used in the main text.

Symbol	Description
$\rho \left(𝐱\right)$	scalar area density, e.g., density of polarity marker on cell surface
$f\left(𝐱\right)$	scalar area density on unit sphere ${𝒮}^2$; projection of $\rho \left(𝐱\right)$
$F_l\left(𝐱\right)$	$l$-th mode of spherical Fourier transform of $f\left(𝐱\right)$, see Eq. (1)
$f_l^m$	$m$-th expansion coefficient of $F_l\left(𝐱\right)$, for expansion into spherical harmonics
${𝒮}^2$	unit sphere
${\int }_{{𝒮}^2}{d}^2𝐱$	integral over unit sphere, using standard Euclidean measure
$S_{ff}\left(l\right)$	spherical power: $L^2$-norm of $l$-th mode $F_l\left(𝐱\right)$ of spherical Fourier transform of $f\left(𝐱\right)$
$𝑨$	nematic tensor associated to surface density $f\left(𝐱\right)$, see Eq. (3)
${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$	eigenvalues of $𝑨$
${𝐚}_1$, ${𝐚}_2$, ${𝐚}_3$	eigenvectors of $𝑨$ corresponding to ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ with ${\alpha }_1\ge {\alpha }_3\ge {\alpha }_2$; we refer to ${𝐚}_1$ as bipolar axis and ${𝐚}_2$ as ring axis
$𝐧$, $𝐦$, $𝐥$	principal axes; the tripod of ortho-normal vectors $𝐧$, $𝐦$, $𝐥$ represents a permutation of ${𝐚}_1$, ${𝐚}_2$, ${𝐚}_3$: $𝐧={𝐚}_{\pi \left(2\right)}$, $𝐦={𝐚}_{\pi \left(3\right)}$, $𝐧={𝐚}_{\pi \left(1\right)}$ with $𝐧$ first principal axis, $𝐦$ second principal axis, $𝐥$ third principal axis
$\pi ,\rho \in S_3$	permutations of the indices $(1,2,3)$
$𝑸$, $𝑩$	traceless tensors, which characterize the tripod $𝐧$, $𝐦$, $𝐥$, see Eq. (4)
${𝑹}_Q$, ${𝑹}_B$	rotation matrices that diagonalize $𝑸$ and $𝑩$, respectively, see Eqs. (8,12)
$𝐰$, $𝐯$, $𝐮$	reference axes; derived from either a common eigenframe ${𝐞}_1$, ${𝐞}_2$, ${𝐞}_3$ of the tensors $𝑸$ and $𝑩$, $𝐰={𝐞}_{\rho \left(2\right)}$, $𝐯={𝐞}_{\rho \left(3\right)}$, $𝐮={𝐞}_{\rho \left(1\right)}$ (OOP), or from a second set of nematic tensors $𝑬$ with eigenvalues ${\epsilon }_1$, ${\epsilon }_2$, ${\epsilon }_3$ that are ordered such that ${\epsilon }_1\ge {\epsilon }_3\ge {\epsilon }_2$, and corresponding eigenvectors ${𝐞}_1$, ${𝐞}_2$, ${𝐞}_3$: $𝐰={𝐞}_2$, $𝐯={𝐞}_3$, $𝐮={𝐞}_1$ (COOP); in both cases, we refer to the orthonormal vectors $𝐰$, $𝐯$, $𝐮$ as $𝐰$ first reference axis, $𝐯$ second reference axis, $𝐮$ third reference axis
$S$, $P$, $D$, $C$	orientational order parameters (OOP) that characterize biaxial order in an ensemble of tripods ${𝐚}_1^{\left(i\right)}$, ${𝐚}_2^{\left(i\right)}$, ${𝐚}_3^{\left(i\right)}$ (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. (13)
$\mathrm{co}-S$, $\mathrm{co}-P$, $\mathrm{co}-D$, $\mathrm{co}-C$	co-orientational order parameters (COOP), introduced here, corresponding to a fixed ordering of principal axes, $𝐧={𝐚}_2$, $𝐦={𝐚}_3$, $𝐥={𝐚}_1$, derived from the eigenvectors ${𝐚}_1$, ${𝐚}_2$, ${𝐚}_3$ of a nematic tensor $𝑨$ with corresponding eigenvalues ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ that are ordered such that ${\alpha }_1\ge {\alpha }_3\ge {\alpha }_2$, and a fixed ordering of reference axes $𝐰={𝐞}_2$, $𝐯={𝐞}_3$, $𝐮={𝐞}_1$ derived from the eigenvectors ${𝐞}_1$, ${𝐞}_3$, ${𝐞}_2$ of a second nematic tensor $𝑬$ with corresponding eigenvalues ${\epsilon }_1$, ${\epsilon }_2$, ${\epsilon }_3$ that are ordered such that ${\epsilon }_1\ge {\epsilon }_3\ge {\epsilon }_2$, see Eq. (16)
$𝑺$	nematic tensor of local anisotropy of sinusoidal network, see Eq. (17)
${𝐬}_1$, ${𝐬}_2$, ${𝐬}_3$	eigenvectors of $𝑺$, corresponding to eigenvalues ${\epsilon }_1$, ${\epsilon }_2$, ${\epsilon }_3$ with ${\epsilon }_1\ge {\epsilon }_3\ge {\epsilon }_2$
$\mathrm{𝟙}$	identity tensor with components ${\mathrm{𝟙}}_{\alpha \beta }={\delta }_{\alpha \beta }$
$H$	dimensionless Hamiltonian of minimal interaction model, see Eq. (18)
$\lambda$	effective interaction parameter in $H$
$𝑰$	moments-of-inertia tensor, see Eq. (S3)
$a$, $b$, $c$	side-lengths of equivalent cuboid for the visualization of nematic tensors $𝑨$, see SI text S4; the convention ${\alpha }_1\ge {\alpha }_3\ge {\alpha }_2$ for the eigenvalues ${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$ of $𝑨$ implies $a\ge b\ge 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 2. Table S2: Transformation of the orientational order parameters S 𝑆 S , P 𝑃 P , D 𝐷 D , C 𝐶 C , under a permutation π ∈ S 3 𝜋 subscript 𝑆 3 \pi\in S_{3} of the principal axes, or permutation ρ ∈ S 3 𝜌 subscript 𝑆 3 \rho\in S_{3} of the reference axes. Permutations are shown in one-line notation (i.e., second row of Cauchy’s two-line notation).

Equations

$$f(\mathbf{x})=\sum _{l=0}^{\infty }{F}_l(\mathbf{x}).$$

$$S_{ff}(l)=\frac{1}{4\pi }{\int }_{{\mathcal{S}}^2}{\mathrm{d}}^2\mathbf{x}|F_l{|}^2.$$

$$\mathbf{A}=\frac{1}{2}{\int }_{{\mathcal{S}}^2}{\mathrm{d}}^2\mathbf{x}f(\mathbf{x})\left(3\mathbf{x}\otimes \mathbf{x}-\text{mathbbm}1\right),$$

$$\mathbf{Q}=\frac{1}{2}\left(3\mathbf{n}\otimes \mathbf{n}-\text{mathbbm}1\right),\mathbf{B}=\frac{3}{2}\left(\mathbf{l}\otimes \mathbf{l}-\mathbf{m}\otimes \mathbf{m}\right).$$

$${\mathbf{R}}_Q^T⟨\mathbf{Q}⟩{\mathbf{R}}_Q$$

$${\mathbf{R}}_B^T⟨\mathbf{B}⟩{\mathbf{R}}_B$$

$$S$$

$$P$$

$$D$$

$$C$$

$$\mathbf{n}={\mathbf{a}}_{\pi (2)},\mathbf{m}={\mathbf{a}}_{\pi (1)},\mathbf{l}={\mathbf{a}}_{\pi (3)},\pi \in S_3,$$

$$|S|\text{ is maximal, }P\ge 0,\text{ and }C\ge 0.$$

$$\mathrm{co}\text{mathchar}45S$$

$$\mathrm{co}\text{mathchar}45P$$

$$\mathrm{co}\text{mathchar}45D$$

$$\mathrm{co}\text{mathchar}45C$$

$$\mathbf{S}({\mathbf{x}}^{(i)})=\sum _{k}w({\mathbf{x}}_k-{\mathbf{x}}^{(i)})l_k\left({\mathbf{d}}_k\otimes {\mathbf{d}}_k-\frac{1}{3}\text{mathbbm}1\right).$$

$$H$$

$$S=\sqrt{\frac{4\pi }{5}}{f}_2^{0}P=\sqrt{6\frac{4\pi }{5}}{f}_2^2$$

$$\displaystyle S_{f\!f}(2)=\frac{1}{4\pi}\sum_{m=-2}^{2}|f_{2}^{m}|^{2}={\color[rgb]{0,0,0}\frac{5}{(4\pi)^{2}}}\left(S^{2}+\frac{1}{3}P^{2}\right)\quad.$$

$$\mathbf{I}={\int }_{{\mathcal{S}}^2}{\mathrm{d}}^2\mathbf{x}\left(\text{mathbbm}1-\mathbf{x}\otimes \mathbf{x}\right)\rho (\mathbf{x}),$$

$$\mathbf{I}=\frac{2}{3}\left(F_0\text{mathbbm}1-\mathbf{A}\right),$$

$$a^2$$

$$b^2$$

$$c^2$$

$$⟨\mathbf{A}{⟩}_{\mathrm{loc}}(\mathbf{x})=\sum _{i\ne j}\frac{1}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{|{\mathbf{x}}^{(i)}-\mathbf{x}{|}^2}{2{\sigma }^2}\right){\mathbf{A}}^{(i)}.$$

$$\text{model variant (a): }H$$

$$\text{model variant (b): }H$$

$$\mathbf{T}$$

$$\mathbf{V}$$

$$I_1$$

$$I_4$$

$$I_1$$

$$I_4$$

$$\mathbf{A}={\xi }_0\mathbf{Q}+{\xi }_1\mathbf{B},$$

$$\mathbf{V}={\zeta }_0⟨\mathbf{Q}⟩+{\zeta }_1⟨\mathbf{B}⟩+{\zeta }_c\text{mathbbm}1,$$

$$\begin{array}[]{ccccccc}2\,\mu_{1}&=&&-\xi_{0}S&+\xi_{0}P&\phantom{2}-\xi_{1}D&+3\,\xi_{1}C\ ,\\ 2\,\mu_{2}&=&&-\xi_{0}S&-\xi_{0}P&\phantom{2}-\xi_{1}D&-3\,\xi_{1}C\ ,\\ 2\,\mu_{3}&=&&2\xi_{0}S&&+2\xi_{1}D&\ ,\\[2.84526pt] 2\,\nu_{1}&=&{\color[rgb]{0,0,0}2}\zeta_{c}&\phantom{2}-\zeta_{0}S&+\zeta_{0}P&\phantom{2}-\zeta_{1}D&+3\,\zeta_{1}C\ ,\\ 2\,\nu_{2}&=&{\color[rgb]{0,0,0}2}\zeta_{c}&\phantom{2}-\zeta_{0}S&-\zeta_{0}P&\phantom{2}-\zeta_{1}D&-3\,\zeta_{1}C\ ,\\ 2\,\nu_{3}&=&{\color[rgb]{0,0,0}2}\zeta_{c}&+2\zeta_{0}S&&+2\zeta_{1}D&\quad.\end{array}$$

Full text

**Quantification of Nematic Cell Polarity in Three-dimensional Tissues

** André Scholich1, Simon Syga1,2, Hernán Morales-Navarrete3, Fabián Segovia-Miranda3, Hidenori Nonaka3, Kirstin Meyer3, Walter de Back2,4, Lutz Brusch2, Yannis Kalaidzidis3, Marino Zerial3,5,6, Frank Jülicher1,5,6, Benjamin M. Friedrich5,6,7*

1 Max Planck Institute for the Physics of Complex Systems, Dresden, Germany

2 Centre for Information Services and High Performance Computing, TU Dresden, Dresden, Germany

3 Max Planck Institute of Molecular Cell Biology and Genetics, Dresden, Germany

4 Institute for Medical Informatics and Biometry, Faculty of Medicine Carl Gustav Carus, TU Dresden, Dresden, Germany

5 Center for Advancing Electronics Dresden, TU Dresden, Germany

6 Cluster of Excellence Physics of Life, TU Dresden, Germany

7 Institute for Theoretical Physics, TU Dresden, Germany

• benjamin.m.friedrich@tu-dresden.de

Abstract

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-Navarette et al., eLife 8:e44860, 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 or in-vitro organoids.

Author Summary

Cell polarity enables cells to carry out specific functions. Cell polarity is characterized by the formation of different plasma membrane domains, each with specific composition of proteins, phospholipids and cytoskeletal components. In simple epithelial sheets, or tube-like tissues such as kidney, epithelial cells are known to display a single apical domain, facing a lumenal cavity, and a single basal domain on the opposite side of the cell, facing a basal layer of extracellular matrix. This apico-basal polarity defines a vector of cell polarity, which provides a direction of fluid transport, e.g., from the basal side of the sheet to the lumen-facing side. In more complex, three-dimensional epithelial tissues, such as liver tissue with its complex network of blood-transporting sinusoids, the membrane domains of hepatocyte cells display more intricate patterns, including rings and antipodal pairs of apical membrane. Here, we develop a mathematical framework to precisely characterize and quantify complex polarity patterns. Thereby, we reveal ordered patterns of cell polarity that span across a liver lobule. Our new method builds on physical concepts originally developed for ordered phases of liquid crystals. It provides a versatile tool to characterize the spatial organization of a complex three-dimensional tissue.

Introduction

In multi-cellular organisms, almost all tissue cells are spatially asymmetric to serve their function inside their host tissue [1]. This cell polarity can be realized by different kinds of physical anisotropies, including cell shape, the structural polarity of their cytoskeleton [2], or the protein and lipid composition within the cell membrane [3, 4].

Here, we focus on the anisotropic distribution of functional membrane domains on the surface of cells, and use the term cell polarity specifically for this important case. A prototypical example is the distribution of polarity-specific apical and basal membrane proteins on the surface of epithelial cells [1].

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. Apical domains are separated by tight junctions from other membrane domains, such as lateral and basal domains. Lateral domains provide cell-cell adhesion, while basal domains form the interface with the basement membrane and extracellular matrix [3, 4]. This structural asymmetry of apical and basal domains in simple epithelia defines a vectorial cell polarity (sometimes referred to as columnar polarity [4]), see Fig. 1A. This vectorial cell polarity sets a direction for the directed transport of macromolecules.

However, there are also epithelial tissues with a more complex, three-dimensional architecture, such as liver tissue [3, 4, 7, 5, 6]. The functional unit of the liver, the liver lobule, is organized around a central and a portal vein, which are connected by a dense, 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 with a volume fraction of approximately 80% [7]. Each hepatocyte is in contact with the sinusoidal network at multiple basal membrane domains, which facilitate the exchange of metabolites with the blood stream [4]. The sinusoidal network was proposed to provide orientational cues to hepatocytes [8, 9]. In addition to the basal contacts, each hepatocyte possesses multiple apical membrane domains that form narrow lumina with adjacent cells, into which bile is excreted [4, 10, 11]. 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 hepatocytes [12, 13]. This model has recently been revisited using high-resolution imaging data of mouse liver tissue [14]. This study showed that the structure of liver tissue is intermediate between an amorphous structure and a perfect crystal, best described as a liquid crystal with orientational but not positional order. The quantification of cell polarity in three-dimensional reconstructions of such high-resolution data prompts new analysis methods to infer the coordination of cell polarity at the tissue level.

We will characterize the distribution of apical membrane domains on the surface of hepatocytes by a tripod of nematic axes. Intuitively, a nematic axis can be thought of as a double-headed arrow that specifies an axis, but does not single out any of the two directions parallel to that axis, see Fig. 1B. 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 [15]. Here, we extend the analysis in [15] to the biaxial case of a full tripod of 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 non-polarized mesenchymal 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 found 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, in line with previous observations [15]. We demonstrate that the coordination of cell polarity is biaxial, i.e., its description requires two local reference axes. Additionally, we find that nematic cell polarity of hepatocytes is correlated with the local biaxial anisotropy of the sinusoidal network. A minimal interaction model conceptualizes the co-alignment of hepatocyte cell polarity and the local anisotropy of the sinusoidal network.

The co-orientational order parameters (COOP) introduced here to characterize the structure of three-dimensional liver tissue naturally generalize previous work on order in effectively two-dimensional cells and tissues. Drew et al. introduced COOPs for two-dimensional systems, and applied this analytical metric to quantify the co-alignment of cytoskeletal structures in muscle cells [2]. Other authors addressed planar cell polarity [16, 17], or nematic alignment of cell shape elongation [18]. The COOP introduced here provide a unified framework to characterize such cellular anisotropies also in three space dimensions.

Results

Nematic cell polarity

We present a 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 can be found in SI text S1.

Let $f(\mathbf{x})$ with $\mathbf{x}\in\mathcal{S}^{2}$ represent an area density on the surface of the unit sphere $\mathcal{S}^{2}$. Similar to the two-dimensional Fourier transform for functions defined on a plane, we decompose the density $f(\mathbf{x})$ into orthogonal modes

[TABLE]

Here, the mode $F_{l}(\mathbf{x})$ of degree $l$ is given by $F_{l}(\mathbf{x})=\sum_{m=-l}^{l}f_{l}^{m}Y_{l}^{m}(\mathbf{x})$, where $Y_{l}^{m}(\mathbf{x})$ 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_{l}^{m}$ are given by $f_{l}^{m}=\int_{\mathcal{S}^{2}}\mathrm{d}^{2}\mathbf{x}\,f(\mathbf{x})\,Y_{l}^{m*}(\mathbf{x})$. Here, integration is performed over the unit sphere $\mathcal{S}^{2}$ (with respect to the standard Euclidean measure); 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(\mathbf{x})$. The first mode $F_{1}(\mathbf{x})$ can be represented by a vector that points to the spherical average of the surface distribution[19]. The second mode $F_{2}(\mathbf{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_{f\!f}(l)$ of each spherical mode $F_{l}(\mathbf{x})$ as its $L^{2}$-norm (normalized by the surface area of the unit sphere)

[TABLE]

This defines the spherical power spectrum, for which a generalized Parseval’s theorem holds, $\sum_{l}S_{f\!f}(l)=\int\mathrm{d}^{2}\mathbf{x}\,|f(\mathbf{x})|^{2}$.

Fig. 1D and Fig. 1E show prototypical vectorial and nematic distributions and their respective spherical power spectra. We visualize surface patterns also as Mollweide projections, an equal-area, pseudocylindrical geographic projection [20]. 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 ring-like pattern as shown in Fig. 1E (and, analogously, for a bipolar pattern with two antipodal caps, see 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 these distributions 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 (i.e., we require that shape can be represented as distance from a common center, which is taken to be the origin, as a single-valued function of the solid angle). This method allows to analyze the anisotropy of surface patterns independent of any anisotropy of cell shape. For hepatocytes, the correlation between cell shape and apico-basal cell polarity is weak, for a discussion see [15].

Specifically, we radially project from the star-convex shape to a unit sphere concentric with the shape, retaining the nominal value of the original distribution, see SI text S2 for additional details. To each cell with surface distribution $\rho(\mathbf{x})$ of (apical) polarity proteins, we associate the projection $f(\mathbf{x})$ of this distribution on the unit sphere $\mathcal{S}^{2}$. 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. Apical plasma membrane domains of these cells from kidney and liver tissue, respectively, were identified by staining fixed cells with anti-CD13 (Novus, cat NB100-64843, rat, 1/500) as reported previously [14, 15]. 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 [4]. This observation from a typical kidney cell is confirmed for an ensemble of cells ($n=286$); note that the relative magnitudes of spherical power modes exhibit a consistent pattern, whereas their absolute magnitude scales with the square of the projected area fraction of apical membrane.

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. Note that the large ensemble of cells analyzed here may contain segmentation errors; thus the ensemble-averaged spherical power spectra presented in Fig. 1FG represent a lower bound.

This analysis of spherical power spectra 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 ${\boldsymbol{A}}$ of the spherical distribution $f(\mathbf{x})$,

[TABLE]

where $\mathbbm{1}$ denotes the identity tensor with components $\mathbbm{1}_{\alpha\beta}=\delta_{\alpha\beta}$. The nematic tensor ${\boldsymbol{A}}$ encodes the same information as the second multipole $F_{2}(\mathbf{x})$. More generally, there is a formal link between the spherical modes of order $l$ and the reduced Cartesian multipole moments [21], see also SI text S3. The nematic tensor ${\boldsymbol{A}}$ is closely related to a moments-of-inertia tensor, see SI test S4. We order the eigenvalues $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ of ${\boldsymbol{A}}$ such that $\alpha_{1}\geq\alpha_{3}\geq\alpha_{2}$ holds and denote the eigenvectors corresponding to $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ by $\mathbf{a}_{1}$, $\mathbf{a}_{2}$, $\mathbf{a}_{3}$ Motivated by Fig. 2AB, we will refer to $\mathbf{a}_{2}$ as the ring axis and $\mathbf{a}_{1}$ as the bipolar axis. Below, the axis $\mathbf{a}_{2}$ will represent an example of a first principal axis used to define co-orientational order parameters, while the axis $\mathbf{a}_{1}$ will represent the second principal axis. The numbering of axes $\mathbf{a}_{2}$ and $\mathbf{a}_{1}$ was chosen to be consistent with [15] (there $\alpha_{i}=\sigma_{i}$, $i=1,2,3$).

Cuboid representation 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 traceless moments-of-inertia tensor of the spherical distribution $f(\mathbf{x})$, see SI text S4 for details. Briefly, the edges of the cuboid are parallel to the eigenvectors of the nematic tensor ${\boldsymbol{A}}$ associated to $f(\mathbf{x})$, while the side-lengths of the cuboid depend on the eigenvalues of ${\boldsymbol{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 $\mathbf{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 $\mathbf{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 $\mathbf{a}_{1}$, and blue to the ring axis $\mathbf{a}_{2}$. Fig. 2C shows the cuboid representation of a typical hepatocyte (with apical membrane distribution shown in green).

Using this cuboidal representation, we can visualize biaxial nematic cell polarity with respect to apical membrane distribution 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 (orange) and the central vein (cyan). In top view, most of the polarity cuboids are faced with their blue side up (indicating an approximately parallel alignment of the ring axis $\mathbf{a}_{2}$ of hepatocyte polarity with the large veins; for the chosen tissue sample these are approximately parallel to the $z$ axis). 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.

Order parameters of nematic cell polarity

We can quantify orientational order of nematic cell polarity within a tissue in terms of orientational order parameters (OOP), $S$, $P$, $D$, $C$. These order parameters were originally developed for the study of biaxial order in liquid crystals [22, 23]. We briefly review their definition and provide illustrative examples to convey their geometric meaning.

Before we present the formal definition of $S$, $P$, $D$, $C$, we want to motivate the different roles played by the four OOPs. For an ensemble of uniaxial objects, each characterized by a single principal axis, $S$ quantifies how well the ensemble of these principal axes are aligned to a common mean direction. The mean direction of this nematic alignment defines one of several reference axes introduced below. In the theory of simple liquid crystals, $S$ is the most widely used OOP (and probably the most important one). Like $S$, the OOP $P$ is also defined already for an ensemble of uniaxial objects (with a single principal axis each). While $S$ characterizes the strength of nematic alignment of the single principal axis in the ensemble, $P$ characterizes anisotropic fluctuations of this principal axis, i.e., $P$ becomes non-zero if the deviations of the principal axes from the mean direction are skewed in a particular direction. The direction of anisotropic fluctuations defines a second reference axis. The two other OOPs, $D$ and $C$, are only defined for biaxial objects that are characterized by two principal axes (a third principal axis can be deduced from the two). The OOP $C$ characterizes nematic alignment of the second principal axis towards a second reference axis. Finally, $D$ quantifies an ‘unexpected relation’ between the fluctuations of the first and second principal axis, respectively, with respect to their alignment to a first reference axis of mean alignment. For a maximum-entropy distribution of orientations, constrained to specific numerical values of $S$, $P$, $C$, the OOP $D$ would be zero.

In liquid crystals, built up by an ensemble of anisotropic molecules, a non-zero value of $S$ can arise from the interactions of uniaxial molecules (e.g., molecules approximated as a rod) with a uniaxial external field. Here, the axis of the external field sets the direction of mean alignment, i.e., the first reference axis. A non-zero value of $P$, however, additionally requires either a second external field, orthogonal to the first one, or boundary conditions that break rotational symmetry for rotations around the first reference axis. Non-zero values of $D$ and $C$ obviously require objects that are intrinsically biaxial, i.e., that do not possess rotational symmetry around their first principal axis. While a non-zero value of $D$ can result already from interactions of biaxial objects with an external uniaxial field, a non-zero value of $C$ requires either an external field that is biaxial, boundary conditions that break uniaxial symmetry, or, possibly, biaxial inter-molecular interactions.

We now proceed to the formal definition of OOPs, which will later be generalized to the case of co-orientational order parameters and applied to quantitatively characterize design principles of liver tissue. We consider an ensemble of cuboids, each characterized by a tripod of principal axes, which we characterize by orthonormal vectors $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$. These unit vectors are only defined up to sign, thus any meaningful physical quantity should be invariant under sign flips $\mathbf{n}\rightarrow-\mathbf{n}$, $\mathbf{m}\rightarrow-\mathbf{m}$, $\mathbf{l}\rightarrow-\mathbf{l}$. Formally, each cuboid is said to possess so-called $D_{2h}$-symmetry, i.e., it is invariant under line reflections at its principal axes. We distinguish a first principal axis $\mathbf{n}$, as well as a second principal axis $\mathbf{m}$, and a third principal axis $\mathbf{l}$. It is convenient to introduce, for each cuboid, two traceless tensors ${\boldsymbol{Q}}$ and ${\boldsymbol{B}}$ that characterize its tripod of axes [24, 23]

[TABLE]

Here, $\mathbbm{1}$ is the identity tensor and $\otimes$ denotes the outer product.

If the principal axes $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$ of each tripod are given by the normalized eigenvectors of a nematic tensor ${\boldsymbol{A}}$ (e.g., the nematic tensor of a projected surface distribution $f(\mathbf{x})$ of membrane proteins), we can recover ${\boldsymbol{A}}$ as linear superposition of the two tensors ${\boldsymbol{Q}}$ and ${\boldsymbol{B}}$, see SI text S8. The mathematical advantage of the traceless tensors ${\boldsymbol{Q}}$ and ${\boldsymbol{B}}$ is that they conveniently allow to compute ensemble averages, $\langle{\boldsymbol{Q}}\rangle$ and $\langle{\boldsymbol{B}}\rangle$. The eigenvalues of these averaged tensors provide important invariants of orientational order

[TABLE]

Here, ${\boldsymbol{R}}_{Q}$ and ${\boldsymbol{R}}_{B}$ are rotation matrices that diagonalize $\left\langle{\boldsymbol{Q}}\right\rangle$ and $\left\langle{\boldsymbol{B}}\right\rangle$, respectively.

In principle, the rotation matrices ${\boldsymbol{R}}_{Q}$ and ${\boldsymbol{R}}_{B}$ might be different. However, for important special cases, e.g., ensembles of biaxial molecules interacting with simple external fields, both rotation matrices can be chosen equal, ${\boldsymbol{R}}={\boldsymbol{R}}_{Q}={\boldsymbol{R}}_{B}$. In this case, the ensemble-averaged tensors $\langle{\boldsymbol{Q}}\rangle$ and $\langle{\boldsymbol{B}}\rangle$ possess the same eigenvector basis, given by three mutually orthogonal unit vectors $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$. The physical meaning of $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$ is that these vectors define mutually orthogonal symmetry axes, such that the statistics of the ensemble of cuboids is invariant under line reflections at these axes. The ensemble of cuboids is said to possess $D_{2h}$-symmetry in this case. In the case of $D_{2h}$-symmetry, the common eigenvectors $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$ of $\left\langle{\boldsymbol{Q}}\right\rangle$ and $\left\langle{\boldsymbol{B}}\right\rangle$ define a director frame of reference axes of the ensemble. Note ${\boldsymbol{R}}_{Q}={\boldsymbol{R}}_{B}=[\mathbf{u},\mathbf{v},\mathbf{w}]^{T}$ (with row-vectors $\mathbf{u}^{T}$, $\mathbf{v}^{T}$, $\mathbf{w}^{T}$).

Using Eq. (8) and Eq. (12), we can rewrite $S$, $P$, $D$, $C$ as averaged direction cosines [23]

[TABLE]

where $\langle\cdot\rangle$ averages over an ensemble of nematic axes $\mathbf{n}^{(i)}$, $\mathbf{m}^{(i)}$, $\mathbf{l}^{(i)}$ indexed by $i$, using a fixed tripod of reference axes $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$.

In Eq. (13), there is ambiguity regarding the ordering of the various axes, $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$, and $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$. For our specific case, the mapping from the nematic cell polarity axes $\mathbf{a}_{1}$, $\mathbf{a}_{2}$, $\mathbf{a}_{3}$ to the principal axes $\mathbf{n},\mathbf{m},\mathbf{l}$ is only defined up to a constant permutation

[TABLE]

where $S_{3}$ denotes the group of all permutations of the indices $(1,2,3)$. Note that the same permutation $\pi$ must be used for all tripods of nematic axes $\mathbf{a}_{1}^{(i)}$, $\mathbf{a}_{2}^{(i)}$, $\mathbf{a}_{3}^{(i)}$ of the ensemble. Similarly, if $\mathbf{e}_{1}$, $\mathbf{e}_{2}$, $\mathbf{e}_{3}$ denote the common eigenvectors of $\langle{\boldsymbol{Q}}\rangle$ and $\langle{\boldsymbol{B}}\rangle$ (with some fixed ordering), we distinguish a first reference axis $\mathbf{w}=\mathbf{e}_{\rho(2)}$ from a second reference axis $\mathbf{v}=\mathbf{e}_{\rho(1)}$, and a third reference axis $\mathbf{u}=\mathbf{e}_{\rho(3)}$, where $\rho\in S_{3}$ denotes a permutation of reference axes. (The axes $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$ are also called director axes [23].)

We distinguish two different choices for $\pi$ and $\rho$, which give rise to orientational order parameters (OOP) as commonly used in the theory of liquid crystals [23], and co-orientational order parameters (COOP) introduced here.

A common choice, put forward, e.g., by Zannoni et al. [23] in the field of biaxial nematics, is to chose the permutations $\pi$ and $\rho$ of principal and reference axes such that

[TABLE]

This condition specifies an ordering of both principal and reference axes. The tensor $\left\langle{\boldsymbol{Q}}\right\rangle$ and the scalar order parameters $S$ and $P$ quantify alignment of the first principal axis $\mathbf{n}$, whereas the tensor $\left\langle{\boldsymbol{B}}\right\rangle$ and the scalar order parameters $D$ and $C$ quantify the alignment of the second and third principal axes, $\mathbf{m}$ and $\mathbf{l}$. We will refer to the values of $S$, $P$, $D$, $C$ corresponding to the ordering of nematic axes specified by Eq. (15) as order parameters (OOP) without further specification. Note that different normalization conventions for OOP are in use, an overview can be found in [21]. While OOPs are always well-defined, they have a crucial disadvantage: OOPs may change discontinuously if system parameters are smoothly varied due to abrupt changes of either $\pi$ or $\rho$, see SI text S8.

We propose an alternative choice, where the ordering of axes is directly determined by the properties of a nematic tensor ${\boldsymbol{A}}$. In the case of surface distributions on a sphere considered here, we take $\mathbf{n}$ to point in the direction of the ring axis $\mathbf{a}_{2}$, and $\mathbf{m}$ to point in the direction of the bipolar axis $\mathbf{a}_{1}$.

We consider the general case, where for each nematic tensor ${\boldsymbol{A}}^{(i)}$ from an ensemble of tensors indexed by $i$, we additionally have a second nematic tensor ${\boldsymbol{E}}^{(i)}$ for each $i$. Below, we discuss two natural cases of such reference tensors. Let $\varepsilon_{1}$, $\varepsilon_{2}$, $\varepsilon_{3}$ be the eigenvalues of one of the ${\boldsymbol{E}}^{(i)}$, ordered such that $\varepsilon_{1}\geq\varepsilon_{3}\geq\varepsilon_{2}$, and $\mathbf{e}_{1}$, $\mathbf{e}_{2}$, $\mathbf{e}_{3}$ be the corresponding (normalized) eigenvectors. We introduce a tripod of reference axes for each index $i$ as $\mathbf{w}^{(i)}=\mathbf{e}_{2}$, $\mathbf{v}^{(i)}=\mathbf{e}_{1}$, $\mathbf{u}^{(i)}=\mathbf{e}_{3}$, i.e., we chose $\rho$ as the identify permutation, $\rho=\mathrm{id}$. We define co-orientational order parameters (COOP) by generalizing Eq. (13) to this general case, where the reference axes $\mathbf{u}^{(i)}$, $\mathbf{v}^{(i)}$, $\mathbf{w}^{(i)}$ are derived from a set of reference tensors ${\boldsymbol{E}}^{(i)}$

[TABLE]

where $\langle\cdot\rangle$ averages over the ensemble of pairs of tripods indexed by $i$. We propose a scheme to compute reference tensors ${\boldsymbol{E}}^{(i)}$ for the important case, where the tensors ${\boldsymbol{A}}^{(i)}={\boldsymbol{A}}(\mathbf{x}^{(i)})$ depend on spatial position $\mathbf{x}^{(i)}$. For each position $\mathbf{x}^{(i)}$, we define ${\boldsymbol{E}}^{(i)}=\langle{\boldsymbol{A}}(\mathbf{x}^{(i)})\rangle_{\mathrm{loc}}$ using a local average with a “punctured” three-dimensional Gaussian kernel centered at $\mathbf{x}^{(i)}$ (excluding the tensor ${\boldsymbol{A}}^{(i)}$ at the central position $\mathbf{x}^{(i)}$), see SI text S5 for details. 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. 2DE.

Below, we additionally consider a variation of this theme, where the tripod of reference axes $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$ 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, $\mathrm{co}\mathchar 45\relax S$, $\mathrm{co}\mathchar 45\relax P$, $\mathrm{co}\mathchar 45\relax D$, $\mathrm{co}\mathchar 45\relax C$, is that the permutation $\pi$ of principal axes, and $\rho$ of reference axes is determined by the orientational order of the ensemble itself for OOP, but prescribed by the eigenvalues of a nematic tensor for COOP. This apparently small change in the mathematical definition renders co-orientational order parameters (COOP) a robust metric that is applicable also in the case of curved director fields, where classical order parameters (OOP) may change discontinuously.

Below, we will apply these biaxial co-orientational order parameters to quantify the alignment of nematic cell polarity of hepatocytes in liver tissue.

Geometric meaning of order parameters.

We illustrate the geometric meaning of the orientational order parameters introduced in Eq. (13), see Fig. 3. The case of co-orientational order parameters defined in Eq. (16) is analogous if principal axes $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$ are plotted relative to the reference axes $\mathbf{w}$, $\mathbf{v}$, $\mathbf{u}$.

When $S>0$ and all other order parameters vanish, as in panel 3A, the ensemble is said to possess uniaxial prolate order (also called cluster-type order [25]). Such uniaxial orderings are axially symmetric around their first reference axis $\mathbf{w}$. If fluctuations of the first principal axis $\mathbf{n}$ are anisotropic, as shown in panel 3B, the ensemble is said to possess phase-biaxial order. This is quantified by the magnitude of the order parameter $P$. In panel 3C, an axially-symmetric distribution with $S<0$ is shown, termed uniaxial oblate order, where the first principal axis $\mathbf{n}$ scatters close to the ‘equator’ (with north-pole south-pole axis set by $\mathbf{w}$). This type of uniaxial order is occasionally called girdle order [25].

So far, we only examined the distribution of the first principal axis $\mathbf{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, $\mathbf{n}$, $\mathbf{m}$, $\mathbf{l}$. In panels 3D, E, and F, we show examples of an additional ordering of a second principal axis $\mathbf{m}$, which are quantified by the other two order parameters $D$ and $C$. We illustrate distributions of the second principal axis $\mathbf{m}$ by antipodal pairs of red points on the sphere. Panel D shows the reference case of an uniaxial prolate distribution as in panel 3A. In absence of any additional ordering, the axis $\mathbf{m}$ displays uniaxial oblate order, as it is must be perpendicular to the first principal axis $\mathbf{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 $\mathbf{m}$. In panel 3E, $\mathbf{m}$ aligns towards the second reference axis $\mathbf{v}$ (green). This breaks axial symmetry around $\mathbf{w}$ for the second principal axis $\mathbf{m}$ (red), but not for the first principal axis $\mathbf{n}$ (blue). Correspondingly, the order parameter $P$ describing the phase biaxiality of the first principal axis $\mathbf{n}$ remains zero, but the molecular biaxiality parameter $C$ becomes nonzero. This parameter thus describes the deviation from axial symmetry with respect to the first reference axis $\mathbf{w}$ of the ensemble of second principal axes $\mathbf{m}$. In contrast, both the first and second principal axis, $\mathbf{n}$ and $\mathbf{m}$, compete for the same reference axis $\mathbf{w}$ in panel F. Correspondingly, their respective distributions remain axially symmetric around $\mathbf{w}$. In this case, both $P$ and $C$ are zero, yet the molecular ordering parameter $D$ is non-zero.

Application to liver tissue

We now apply the framework of biaxial order parameters to quantify lobule-level 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 nematic polarity axes of each individual hepatocyte (introduced in Fig. 2D) with a local reference frame, given by a local average of the tensors ${\boldsymbol{A}}$ with a punctured Gaussian kernel (illustrated in Fig. 2F). This provides reference axes $\mathbf{w}=\langle\mathbf{a}_{2}\rangle_{\mathrm{loc}}$, $\mathbf{v}=\langle\mathbf{a}_{1}\rangle_{\mathrm{loc}}$, and $\mathbf{u}=\mathbf{v}\times\mathbf{w}$ at each hepatocyte position. (Mathematically, $\langle\cdot\rangle_{\mathrm{loc}}$ is defined using nematic tensors, see SI text S5 for details.) We choose the first principal axis $\mathbf{n}$ to point in the direction of the ring axis $\mathbf{a}_{2}$, and the second principal axis $\mathbf{m}$ to point in the direction of the bipolar axis $\mathbf{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 $\mathbf{a}_{2}$ (blue dots) is clustered around the first reference axis $\mathbf{w}=\langle\mathbf{a}_{2}\rangle_{\mathrm{loc}}$. Correspondingly, the scalar order parameter $\mathrm{co}\mathchar 45\relax S$ of uniaxial nematic order is larger than zero. Additionally, we find a statistically significant phase biaxiality with $\mathrm{co}\mathchar 45\relax P>0$, revealing that fluctuations of the ring axis $\mathbf{n}=\mathbf{a}_{2}$ are biased away from the average bipolar axis $\mathbf{v}=\langle\mathbf{a}_{1}\rangle_{\mathrm{loc}}$. This phase biaxiality is also visible in the distribution plot on the sphere in Fig. 3H. The second principal axis $\mathbf{m}$ (bipolar axis $\mathbf{a}_{1}$, red dots) also exhibits a weak ordering with a bias towards $\mathbf{v}$ (averaged bipolar axis $\langle\mathbf{a}_{1}\rangle_{\mathrm{loc}}$, green) and away from $\mathbf{w}$, reflected by positive values of $\mathrm{co}\mathchar 45\relax C$ and $\mathrm{co}\mathchar 45\relax D$, respectively. 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 [15, 26]. Sinusoids are specialized blood vessels forming a network within the liver lobule [10]. 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 $\mathbf{d}_{k}$ are unit vectors parallel to straight network segments, $\mathbf{x}_{k}$ their midpoint positions and $l_{k}$ their respective lengths, we define nematic tensors ${\boldsymbol{S}}$ at each hepatocyte position $\mathbf{x}^{(i)}$

[TABLE]

Here, $w(\mathbf{x})$ is a weighting function normalized as $\sum_{k}w(\mathbf{x}_{k}-\mathbf{x}^{(i)})l_{k}=1$. We choose $w(\mathbf{x})$ as a binary cutoff with fixed radius of 20 $\mu$m around the center of each hepatocyte. The geometric meaning of ${\boldsymbol{S}}$ can be understood as follows: The eigenvector $\mathbf{s}_{1}$, corresponding to the largest eigenvalue, characterizes the direction of preferred sinusoid orientation and will be referred to as preferred axis. The eigenvector $\mathbf{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. The biaxial anisotropy of the sinusoidal network with a distinguished plane axis is indicative of a local layered order, where $\mathbf{s}_{2}$ represents the normal vector of a stack of parallel layers, which characterizes this layered order. Approximately, $\mathbf{s}_{2}$ is parallel to both the centerline of the portal vein, and the centerline of the central vein. 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 $\mathbf{a}_{2}$ of apical cell polarity is well-aligned with the plane axis $\mathbf{s}_{2}$ of the local sinusoid anisotropy. For our choice of axes, this is quantified by the order parameter $\mathrm{co}\mathchar 45\relax S$. We also find phase-biaxiality of this axis, reflected by a non-zero value of $\mathrm{co}\mathchar 45\relax P$. The other co-orientational order parameters $\mathrm{co}\mathchar 45\relax D$ and $\mathrm{co}\mathchar 45\relax C$ are close to zero, i.e., we do not find a particular ordering of the bipolar cell polarity axis $\mathbf{a}_{1}$ relative to ${\boldsymbol{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 $\mathbf{a}_{2}$ of hepatocytes, the leading order term of an effective interaction energy is dictated by symmetry and reads

[TABLE]

Here, ${\boldsymbol{A}}:{\boldsymbol{E}}$ denotes the contraction of two tensors ${\boldsymbol{A}}$ and ${\boldsymbol{E}}$. Thus, we treat the hepatocytes as uniaxial objects, whereas we retain the biaxial anisotropy of the sinusoid network. The choice of Eq. (18) is motivated by a general Landau theory of liquid crystals, see [24, 23]. We calculate the order parameters of an ensemble of axes according to the Boltzmann distribution following this Hamiltonian, using inverse sampling. The control parameter $\lambda$ is measured in units of an effective temperature that mimics dynamic processes that reduce spatial order [27]. We emphasize that we do not consider liver tissue to represent a thermodynamic equilibrium, despite our use of Eq. (18). Instead, Eq. (18) represents a phenomenological model that addresses a competition between dynamic processes that either generate or reduce spatial order, respectively.

Figure 5A displays computed COOPs, 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 $\lambda$ (shaded gray in Fig. 5A), where the minimal model adequately accounts for the experimental observed values of the co-orientational order parameters. Thus, the interaction of the ring axis $\mathbf{a}_{2}$ of hepatocytes with the local anisotropy of the sinusoidal network is sufficient to account for the observed biaxial co-orientation. Intriguingly, alternative models assuming either an interaction between ${\boldsymbol{S}}$ and the bipolar axis $\mathbf{a}_{1}$, or the full tensor ${\boldsymbol{A}}$, did not reproduce the observed co-orientational order, see SI text S6.

This finding suggests the cartoon picture of sinusoid-hepatocyte co-alignment in liver tissue shown in Fig. 5B. We propose that the ring axis $\mathbf{a}_{2}$ of hepatocytes preferentially aligns parallel to the plane axis $\mathbf{s}_{2}$ of the local sinusoidal network. Fluctuations of $\mathbf{a}_{2}$ break axial symmetry and are biased away from the preferred axis $\mathbf{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 [14]. 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 [15]. 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 [15].

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 [28, 29, 23], and generalized these as co-orientational order parameters, which we then apply to quantify structural order in living matter. Traditionally, OOP are used to quantify the partial alignment of anisotropic molecules, where each molecules characterized by a tripod of nematic axes [22, 31, 30, 32]; here, we use COOP to quantify the partial alignment of nematic cell polarity axes in a tissue.

The co-orientational order parameters (COOP) introduced here have 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 liquid-crystal order of nematic cell polarity of hepatocytes is biaxial, which is a rare finding even in inanimate matter [33]. Furthermore, we found co-alignment between nematic cell polarity of hepatocytes and the local anisotropy of the sinusoidal network. This mutual alignment is not uniaxial, but of phase-biaxial type, i.e., its description requires two reference axes, a preferred axis of the sinusoidal network (approximately parallel to the direction of blood flow [34, 26]), and a plane axis, which characterizes local layered order of the sinusoidal network. We conceptualized this biaxial order using a minimal interaction model, which quantitatively reproduces the COOP observed in the experimental data. The results from our minimal model provide insight into the previous observation of lobule-level spatial order of the bipolar axis ($\mathbf{a}_{1}$) of apical cell polarity of hepatocytes [15]. In fact, our minimal model suggests that the observed order of the bipolar axes can be considered a consequence of the phase-biaxial co-alignment of the second cell polarity axis, the ring axis ($\mathbf{a}_{2}$), relative to the local anisotropy of the sinusoidal network.

Our findings hint at a close interplay between hepatocyte polarity and the local anisotropy of the sinusoidal network. A recent study involving some of the authors showed that interference with communication from sinusoids to hepatocytes disrupts the liquid-crystal order of both hepatocyte cell polarity and the sinusoidal network [15]. 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.

Table S1 contains a list of symbols used in the main text.

S2

Spherical projection of membrane protein density.

We discuss two possible methods to project surface distributions $\rho(\mathbf{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 $\rho(\mathbf{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 $\rho(\mathbf{x})$ on the cell surface yields a projected distribution $f(\mathbf{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 changing the projection method almost did not change computed cell polarity axes for most hepatocytes in liver tissue.

S3

Relation between second mode spherical spectral power and order parameters.

We can consider a distribution $p(\mathbf{n})$ of nematic axes $\mathbf{n}$, represented by antipodal pairs of points, also as a surface distribution on the unit sphere, which is symmetric with respect to a point reflection at the center. Under this correspondence, the uniaxial order parameters $S$ and $P$ for $\mathbf{n}$ as first principal axis are intimately linked to the expansion of $p(\mathbf{n})$ into spherical harmonics, Eq. (1). We assume that $p(\mathbf{n})$ possesses $D_{2h}$-symmetry. Without loss of generality, the reference axes $\mathbf{u}$, $\mathbf{v}$, $\mathbf{w}$ shall be aligned with the $x$, $y$, $z$ axis, respectively. Then, all coefficients $f_{2}^{m}$ of the second spherical mode $F_{2}(\mathbf{x})=\sum_{m=-2}^{2}\,f_{2}^{m}Y_{2}^{m}(\mathbf{x})$ that correspond to functions that are odd relative to a line reflection at the $z$-axis vanish, i.e., $f_{2}^{-1}=f_{2}^{1}=0$ and $\mathrm{Im}f_{2}^{-2}=\mathrm{Im}f_{2}^{2}=0$. Moreover, $f_{2}^{2}=f_{2}^{-2}$, since $p(\mathbf{n})$ is real.

We thus find [21]

[TABLE]

Conversely, the spherical power in the second mode can be expressed in terms of $S$ and $P$

[TABLE]

S4

Cuboid visualization of nematic cell polarity.

In the main text, we present a method to visualize nematic tensors ${\boldsymbol{A}}$ by colored cuboids as shown in Fig. 2. We provide additional details on this method. For a surface density $f(\mathbf{x})$ on the unit sphere $\mathcal{S}^{2}$, the moments-of-inertia tensor ${\boldsymbol{I}}$ reads

[TABLE]

i.e.,

[TABLE]

where ${\boldsymbol{A}}$ is the nematic tensor associated to $f(\mathbf{x})$, see Eq. (3), and $F_{0}=\int_{\mathcal{S}^{2}}\mathrm{d}^{2}\mathbf{x}\,\rho(\mathbf{x})$. Both tensors diagonalize in the same eigenframe. The eigenvalues $\iota_{1}$, $\iota_{2}$, $\iota_{3}$ of ${\boldsymbol{I}}$ (called principal moments of inertia), and the eigenvalues $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$ of ${\boldsymbol{A}}$ are related by $\iota_{i}=(2/3)F_{0}-\alpha_{i}$, $i=1,2,3$ (for a suitable ordering of $\iota_{i}$).

In turn, the principal moments of inertia $\iota_{i}$ of a solid cuboid with side lengths $a$, $b$, $c$ are given by $\iota_{1}=(b^{2}+c^{2})/12$, $\iota_{2}=(a^{2}+b^{2})/12$, $\iota_{3}=(a^{2}+c^{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 ${\boldsymbol{I}}$

[TABLE]

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 nematic tensors ${\boldsymbol{A}}$ of individual hepatocytes by averaging with a Gaussian kernel. Specifically, given nematic tensors ${\boldsymbol{A}}^{(i)}$ at cell center locations $\mathbf{x}^{(i)}$, the coarse-grained tensor $\langle{\boldsymbol{A}}\rangle_{\mathrm{loc}}(\mathbf{x})$ at location $\mathbf{x}$ is calculated by

[TABLE]

Here, $\sigma$ 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 ${\boldsymbol{A}}$, each derived from an individual surface distribution $f^{(i)}(\mathbf{x})$, we could have equivalently averaged the surface distributions first, and then computed $\langle{\boldsymbol{A}}\rangle_{\mathrm{loc}}$ as the nematic tensor of an averaged surface distribution $\langle f(\mathbf{x})\rangle_{\mathrm{loc}}$.

For the principal axes of the averaged tensor $\langle{\boldsymbol{A}}\rangle_{\mathrm{loc}}$, we write $\langle\mathbf{a}_{1}\rangle_{\mathrm{loc}}$, $\langle\mathbf{a}_{2}\rangle_{\mathrm{loc}}$, $\langle\mathbf{a}_{3}\rangle_{\mathrm{loc}}$, for short.

S6

Nematic interaction models.

In addition to the interaction proposed in Eq. (18), two model variants are conceivable: (a) the bipolar axis $\mathbf{a}_{1}$ of hepatocytes could be coupled to the local anisotropy tensor ${\boldsymbol{S}}$ of the sinusoidal network, or (b) the full nematic tensor ${\boldsymbol{A}}$ of hepatocyte polarity, which comprises both the ring and the bipolar axes, could couple to ${\boldsymbol{S}}$, which correspond to respective effective interaction energies

[TABLE]

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. (13) change under a permutation $\pi\in S_{3}$ of the principal axes, $\mathbf{n}=\mathbf{a}_{\pi(2)}$, $\mathbf{m}=\mathbf{a}_{\pi(1)}$, $\mathbf{l}=a_{\pi(3)}$, as well as under a permutation $\rho\in S_{3}$ of the reference axes, $\mathbf{w}=\mathbf{e}_{\rho(2)}$, $\mathbf{v}=\mathbf{e}_{\rho(1)}$, $\mathbf{u}=\mathbf{e}_{\rho(3)}$. The action of the direct product of both permutation groups, $G=S_{3}\times 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 $G$ 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$.

S8

Relation between biaxial order parameters and invariants of moment tensors.

We quantified orientational order of nematic tensors by the four classical order parameters $S$, $P$, $D$, $C$, as well as by co-orientational order parameters $\mathrm{co}\mathchar 45\relax S$, $\mathrm{co}\mathchar 45\relax P$, $\mathrm{co}\mathchar 45\relax D$, $\mathrm{co}\mathchar 45\relax C$. We present yet a third variant to quantify biaxial order: invariants of moment tensors [35].

We consider the first two moments, ${\boldsymbol{T}}$ and ${\boldsymbol{V}}$, of a distribution of nematic tensors ${\boldsymbol{A}}$

[TABLE]

From these averaged tensors, we obtain scalar invariants by tensor contraction

[TABLE]

Note that since ${\boldsymbol{A}}$ is traceless, all non-zero contractions of the rank-4 super-tensor $\langle A_{\alpha\gamma}A_{\delta\beta}\rangle$ can already be derived from the rank-2 tensor ${\boldsymbol{V}}$.

If the ensemble of tensors ${\boldsymbol{A}}^{(i)}$ exhibits $D_{2h}$-symmetry, the moment tensors ${\boldsymbol{T}}$ and ${\boldsymbol{V}}$ diagonalize in a common eigenframe [36]. The invariants $I_{1},\ldots,I_{6}$ can then be expressed in terms of symmetric polynomials in the eigenvalues of these tensors. Specifically, we denote the eigenvalues of ${\boldsymbol{T}}$ by $\mu_{1}$, $\mu_{2}$ and $\mu_{3}$, and the eigenvalues of ${\boldsymbol{V}}$ by $\nu_{1}$, $\nu_{2}$, $\nu_{3}$. Then the invariants of tensor moments are given as

[TABLE]

Conversely, given the invariants $I_{1}$, $I_{2}$, $\dots$, $I_{6}$, we can compute the eigenvalues $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, and $\nu_{1}$, $\nu_{2}$, $\nu_{3}$, yet only up to a permutation, by solving the polynomial system of equations, Eq. (S12).

We now show how the tensor invariants are related to the classical order parameters $S$, $P$, $D$, $C$. We make the simplifying assumption that all nematic tensors ${\boldsymbol{A}}$ of the ensemble have identical eigenvalues $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$. As usual, we assume $\alpha_{1}\geq\alpha_{3}\geq\alpha_{2}$, and denote the corresponding eigenvectors as $\mathbf{n}=\mathbf{a}_{2}$, $\mathbf{m}=\mathbf{a}_{1}$, $\mathbf{l}=\mathbf{a}_{3}$. By Eq. (4), each tensor ${\boldsymbol{A}}$ is now associated with tensors ${\boldsymbol{Q}}$ and ${\boldsymbol{B}}$. We can write

[TABLE]

with weights $\xi_{0}$, $\xi_{1}$ that satisfy $\alpha_{2}=\xi_{0}$, $\alpha_{1}=-(\xi_{0}+3\xi_{1})/2$, $\alpha_{3}=-(\xi_{0}-3\xi_{1})/2$. Note that a permutation of eigenvalues also changes the weights $\xi_{0}$ and $\xi_{1}$. Our usual ordering of eigenvalues with $\alpha_{1}\geq\alpha_{3}\geq\alpha_{2}$ corresponds to $0\geq\xi_{0}/3\geq\xi_{1}$. In the general case, where the eigenvalues of ${\boldsymbol{A}}$ vary within the ensemble, the invariants $I_{j}$ will depend on both the orientational order of the principal axes of ${\boldsymbol{A}}$, as well as on the distribution of weights.

In the case of constant weights $\xi_{0}$, $\xi_{1}$, it follows ${\boldsymbol{T}}=\langle{\boldsymbol{A}}\rangle=\xi_{0}\,\langle{\boldsymbol{Q}}\rangle+\xi_{1}\,\langle{\boldsymbol{B}}\rangle$. Likewise, the second moment ${\boldsymbol{V}}$ can be expressed as a linear superposition of $\langle{\boldsymbol{Q}}\rangle$, $\langle{\boldsymbol{B}}\rangle$, and $\mathbbm{1}$ as

[TABLE]

where $\zeta_{0}=(\xi_{0}^{2}-3\xi_{1}^{2})/2$, $\zeta_{1}=-\xi_{0}\xi_{1}$, $\zeta_{c}=(\xi_{0}^{2}+3\xi_{1}^{2})/2$.

We thus have a direct correspondence between the eigenvalues of the moment tensors and the orientational order parameters $S$, $P$, $D$, $C$

[TABLE]

Note $\mu_{1}+\mu_{2}+\mu_{3}=0$, while $\nu_{1}+\nu_{2}+\nu_{3}=3\zeta_{c}$.

Together, Eq. (S12) and Eq. (S15) allow to compute the orientational order parameters $S$, $P$, $D$, $C$ from the invariants of tensor moments. Note that the first tensor moment is not sufficient to determine the OOP, but that at least the second moment is needed. (In the non-generic case

$\xi_{0}=-1$ and $\xi_{1}=\pm 1$ for which $\zeta_{0}=\xi_{0}$ and $\zeta_{1}=\xi_{1}$,

also a third tensor moment needs to be taken into account.)

We emphasize that the values of the invariants $I_{1}$, $I_{2}$, $\dots$, $I_{6}$ are independent of any ordering of axes, whereas the values of the orientational order parameters $S,P,D,C$ depend on the ordering of both the principal and the reference axes. This is reflected in Eq. (S12) by the fact that $I_{1},I_{2},\dots,I_{6}$ do not change under neither a permutation of the eigenvalues $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, nor a permutation of the eigenvalues $\nu_{1}$, $\nu_{2}$, $\nu_{3}$. In contrast, Eq. (S15) shows that $S$, $P$, $D$, $C$ depend on the ordering of eigenvalues. As a consequence, the orientational order parameters $S$, $P$, $D$, $C$ can change discontinuously if system parameters are smoothly varied, while the invariants $I_{1}$, $I_{2}$, $\dots$, $I_{6}$ do not, see Fig. S4. Despite this desirable property of the invariants $I_{1}$, $I_{2}$, $\dots$, $I_{6}$, the invariants lack the intuitive geometric interpretation of the orientational order parameters $S$, $P$, $D$, $C$. The co-orientational order parameters $\mathrm{co}\mathchar 45\relax S$, $\mathrm{co}\mathchar 45\relax P$, $\mathrm{co}\mathchar 45\relax D$, $\mathrm{co}\mathchar 45\relax C$ introduced in the main text combine the advantageous property of a smooth dependence on system parameters with intuitive geometric interpretation.

Acknowledgments

We thank Samuel Safran for stimulating discussions. AS acknowledges support from the Max-Planck society. Parts of this work were supported by the German Federal Ministry of Education and Research (BMBF) within the research network Systems Medicine of the Liver (LiSyM) [grant 031L0033 to LB, grant 031L0038 to MZ]. Parts of this work were supported by the European Research Council (ERC) (grant number 695646 to MZ). BMF, FJ and MZ acknowledge support by the DFG through the Excellence Initiative by the German Federal and State Governments (Clusters of Excellence cfaed EXC 1056 and PoL EXC 2068).