Origami: Single-cell 3D shape dynamics oriented along the apico-basal axis of folding epithelia from fluorescence microscopy data

A common feature of morphogenesis is the formation of three-dimensional structures from the folding of two-dimensional epithelial sheets, aided by cell shape changes at the cellular-level. Changes in cell shape must be studied in the context of cell-polarised biomechanical processes within the epithelial sheet. In epithelia with highly curved surfaces, finding single-cell alignment along a biological axis can be difficult to automate in silico. We present ‘Origami’, a MATLAB-based image analysis pipeline to compute direction-variant cell shape features along the epithelial apico-basal axis. Our automated method accurately computed direction vectors denoting the apico-basal axis in regions with opposing curvature in synthetic epithelia and fluorescence images of zebrafish embryos. As proof of concept, we identified different cell shape signatures in the developing zebrafish inner ear, where the epithelium deforms in opposite orientations to form different structures. Origami is designed to be user-friendly and is generally applicable to fluorescence images of curved epithelia.


Introduction
Complex morphologies across taxa and tissue types are generated through the deformation of epithelial sheets [1][2][3]. In the embryo, many developing epithelia form highly curved surfaces. Epithelial folding processes are driven by polarised mechanical forces and involve threedimensional changes in shape at the cellular level [4,5]. Fluorescence imaging techniques have made it possible to follow such shape changes at cellular resolution, in vivo and in real-time [6][7][8]. These imaging advances have consequently driven the development of tools to quantify epithelial dynamics, especially cell shape changes.
Solving the orientation of individual cells relative to the known overall polarity of the epithelial sheet is critical, as cell-polarised biomechanical processes drive cell shape changes; constriction or expansion can occur along either the apical [24,25] or baso-lateral [26] cell surfaces and can be detected by any skew in mass distribution within a cell along an apico-basal axis of symmetry. Epithelial folding may be initiated or influenced by cell proliferation, cell death, cytoskeletal remodelling, or changes in cell surface properties [27,28]. These mechanisms can lead to changes in cell shape features, including cell height and width, volume, surface area and sphericity.
Cell orientation or polarity can be defined along the plane of the epithelium (planar cell polarity) or perpendicular to the epithelial plane, along the apico-basal axis of the cell. Existing automated methods for assigning polarity to segmented cells often rely on additional biochemical markers for polarity [29][30][31]. Including such additional markers in fluorescence imaging experiments increases the time taken to generate each image, potentially leading to phototoxicity, and the resulting larger volume of image data makes analysis computationally expensive. Moreover, producing the required animals carrying multiple transgenes for live imaging can be challenging and costly. Some image analysis methods compute direction vectors for individual cells by drawing normal vectors to polynomial functions, often ellipsoids, used to estimate the surface of the specimen-for example, entire embryos [15] or blastoderms [32] undergoing morphogenesis. These methods are specific to the geometry of the specimen and are unsuitable for analysing complex folded topologies at advanced morphologic developmental stages. A third method uses known features of cell shape to assign cell orientation, for example by applying principal component analysis (PCA) to compute the apico-basal axis in columnar cells in EDGE4D [33] and the anterior-posterior axis in zebrafish lateral line primordia using landmark-based geometric morphometrics [31], or orienting cells along their long axis in the zebrafish optic cup as in LongAxis [34]. These strategies will be applicable only if a dominant cell shape feature, for example cell height/width ratio, is known and remains unchanged over space and time.
We introduce a new automated and easy-to-use tool, Origami, for extracting direction-variant shape features along the apico-basal axis by reconstructing the epithelial surface using a triangular mesh (Fig 1). Origami applies to a wide range of geometries of specimens undergoing morphogenesis and automatically extracts direction vectors for individual cells aligned to the apico-basal axis of the epithelial sheet without requiring additional labels for polarity. Direction-variant shape features are calculated by computing the geometric moments for the volume enclosed by the polygon representation of each segmented cell [35]. We showcase the versatility of our method using data from an assortment of structures at a range of developmental stages within the otic vesicle (developing inner ear) of zebrafish embryos.

Design and implementation
The Origami pipeline is preceded by a membrane-based segmentation step. For this, we employed the open-source ACME segmentation software [14]. The segmented data are subjected to two main operations within Origami; epithelial polarity direction vector assignment ( Fig 1B) and extraction of shape features ( Fig 1C).

Assigning polarity direction to individual cells
To compute directionally variant cell shape features, such as skewness (asymmetry in cell mass), and longitudinal and transversal spread, the positioning of segmented cells must be found in 3D space along a biologically relevant axis-we chose the known apico-basal axis of the cell. The folding epithelium was reconstructed in silico as a thin 'crust'-an open surface mesh that triangulates the centroids of the segmented cells in 3D space using the Crust algorithm [36,37] (Fig 1B). The Crust method computes a surface mesh from unorganised points -cell centroids in our case, using the Voronoi diagram of the cell centroids.
Following this, our automated method corrects imperfections in the estimated surface mesh that can cause errors in the resulting direction vectors. The mesh is refined by removing duplications (in vertices or triangular faces computed) and any self-intersecting triangular faces. Non-manifold edges, that is, those edges shared by more than two triangular faces, are re-meshed as a manifold mesh using the ball-pivoting algorithm [38,39].
The triangular faces of the refined mesh are ordered, and so by applying the right-hand rule when generating normal vectors to the surface mesh, these vectors all point to the same side of the mesh representation of the epithelial surface (Fig 1B). At this point, there are still two possible opposing orientations for each computed vector-facing the apical or the basal face of the epithelium, with a difference only in sign. In the developing zebrafish otic vesicle, the apical surface of the epithelium faces the fluid-filled lumen of the vesicle [2,8,40]. We used this prior knowledge to inform the orientation of the vectors by setting the value of a binary orientationdetermining parameter to 'in' so that they point to a convergent point which falls on the side Blue arrows mark the direction of apicobasal polarity (pointing towards the apical side). b. Polarity assignment on segmented data; ROI surrounding the anterior projection was segmented (here overlaid on the MIP) using ACME, centroids were generated for each segmented cell and a triangular surface mesh was produced from these centroids. Normal vectors (blue arrows) to this surface mesh represent the apico-basal axis. c. Cell shape features were computed concerning the assigned apico-basal axis; here, three example cells are highlighted, alongside a 3D rendering showing their position in the anterior projection and the corresponding shape metrics in a table.
of the curved surface mesh that corresponds to the apical (lumenal) side of the epithelium at each cell. When a structure folds multiple times in opposing orientations, such as in the synthetic data generated for this study, 'in' sets the polarity direction vectors to point towards a global convergent point, in our case determined by the curvature of the whole synthetic epithelium.
Under-segmentation can cause missing regions or unwanted holes in the triangular mesh, introducing errors when ordering the triangular faces. Our pipeline attempts to repair these holes by detecting and then remeshing them where possible. Holes, when detected, are flagged as a warning to users about potential errors in the output. Normal vectors to the reconstructed surface represent the apico-basal axis of the epithelium and are generated for each segmented cell at their centroid position (Fig 1B and 1C).

Computing shape features using 3D geometric moments
The shape of an object can be characterised using central geometric moments [41]. Geometric moments are widely used in object recognition and classification problems [42,43] since they (i) are simple to compute, (ii) organise features in orders of increasing detail, and (iii) can be extended to n dimensions. Each moment, G ðVÞ ijk , is defined by the integral over the object (in our case, each segmented cell), of the Cartesian coordinates monomial x i y j z k , where i, j, k � 0, with the origin of coordinates at the centroid.
In our analysis pipeline, 3D geometric moments were computed using the algorithm introduced in [35]. The defining continuous integrals are exactly computed within the triangular surface mesh generated for each individual segmented cell, split into a sum: where each tetrahedron T c is defined by a triangle in the surface mesh and the origin (cell centroid). The determinant gives the oriented volume of this tetrahedron, Considering its sign, the determinant allows the algorithm to be applied to shapes of any complexity and topology. The integral in each T c is given by a closed formula involving only the Cartesian coordinates of the triangular vertices.
The geometric moments of low orders have simple, intuitive interpretations. The zero th order moment G ðVÞ 000 provides the volume of the object, here an individual cell. For central moments, the first order moments are trivially null: The secondorder moments correspond to the spread (covariance tensor) of the distribution. So, the projection of the mass of each cell along the corresponding polarity vector represents the 'spread' as variance in mass 'longitudinally' (along the apico-basal axis) and 'transversally' (along the epithelial plane). This allowed us to identify if cells were more or less columnar (tall cells) or squamous (flat cells) in shape. The third-order moments represent 'skewness', which is the deviation from symmetry. In our pipeline, skewness was measured along the polarity direction vector in the apico-basal direction, with positive skewness values indicating apical cell constriction and/or basal relaxation and negative values indicating basal cell constriction and/or apical expansion. A value of zero indicated no skew. Additionally, the sphericity of each cell was computed as the ratio of the cell surface area to the surface area of a sphere with the same volume as the cell [44], from 0 for a highly irregularly-shaped cell to 1 for a perfect sphere.

Evaluation of computed cell polarity direction vector
To evaluate the computed direction vectors denoting cell polarity, we generated 3D synthetic data representing curved, folding epithelia with varying degrees of curvature and height of folded peak in two opposing orientations (S1 Text and Fig 2A). To reflect real-world in vivo fluorescence imaging conditions, these synthetic data were corrupted with three incremental levels of Gaussian and Poisson noise (S1 Text and Fig 2A). Using the synthetic data, two types of error in computed polarity direction vectors were assessed: (1) an orientation flipping error, measured as the percentage of polarity vectors with an opposing orientation (opposite sign) to the polarity ground truth (S1 Text), and (2) direction accuracy, measured as the mean deviation angle between the polarity vectors produced by Origami, correctly oriented, and the polarity ground truth.
Of the two aspects of surface geometry analysed, height of folded peak (in two opposing directions) did not contribute significantly to orientation flipping errors (Linear Regression; p = 0.86, R 2 = -0.04). However, a larger radius of curvature of epithelium (a flatter epithelial sheet) did correlate with orientation flipping errors-albeit with a small effect of 0.08% increase for every 1 μm (5 pixels) increase in radius of curvature (Linear Regression; p = 0.042, R 2 = 0.12, effect), and a lower quality of segmentation output from ACME (Linear Regression p < 0.001, R 2 = 0.46; Fig 2B) computed as a Dice score. This meant a 0.2% reduction in Dice score for every 1 μm (5 pixels) increase in the radius of curvature. This correlation may be attributed to the reduced ability of ACME to segment flat, squamous cells in an epithelium oriented mostly along the lateral (xy) plane in data with anisotropic voxel resolution (here modelled using an anisotropic point spread function (PSF)). We found a correlation between noise applied to the synthetic images and errors in both polarity orientation flipping (ANOVA: p � 0.001; Tukey's contrasts showed 11.3% increase in errors at highest noise level compared with the lowest noise level applied: p = 0.0039) and segmentation output (ANOVA: p < 0.01; Tukey's contrasts showed 16.3% reduction in Dice score at highest noise level from the lowest noise level applied: p = 0.0045). Segmentation quality, in turn, influenced polarity orientation flipping, with errors below 1.5% at Dice scores above 0.8, but increasing with further decrease in Dice scores (Polynomial Regression; first-order: p < 0.001, Effect size = -28.78; secondorder: p < 0.01, Effect size = 16.26; Fig 2C). Comparisons of many available segmentation algorithms when validating with fluorescent images from non-folded structures such as earlystage nematode embryos [16] or plant roots [18] have been shown to give Dice scores above 80%, suggesting a good performance under real experimental conditions. Quantitative direction accuracy was evaluated in the synthetic data, for which, in contrast to data from real fluorescence images, a reliable ground truth could be generated from the known underlying surface functions. Compared to the polarity ground truth data, an overall offset of 10.6˚± 15.5˚(mean ± std) was measured from our entire synthetic dataset. Just as for the polarity orientation flipping error, height of folded peak did not influence polarity direction accuracy (Linear Regression; p = 0.39, R 2 = -0.01), but there was a small effect of curvature of the epithelium with an additional 0.06˚offset for every 1 μm (5 pixels) increase in radius of curvature of the epithelium (Linear Regression; p = 0.005, R 2 = 0.24). At the highest level of noise applied, errors in polarity orientation had a 6.6˚greater offset than at the lowest noise level applied (Tukey's contrasts; p = 0.003). There was also a negative linear effect of segmentation quality with a 2.9˚offset predicted for every 10% reduction in Dice score (Linear Regression; p < 0.0001, R 2 = 0.53; Fig 2C).
We further tested the effect of such errors in direction accuracy on the direction-variant shape metrics computed by applying directional noise-with a mean equal to the measured mean error above-to polarity vectors of three example cells showing extreme shape features from the synthetic dataset and computed direction-variant shape metrics for each new displaced polarity vector (n = 50; Fig 2E). The resulting computed shape metrics could still successfully differentiate between the three cells, showing that direction accuracy errors (excluding orientation flipping errors) should not adversely affect the shape metrics computed. On the other hand, orientation flipping errors will affect shape metrics, but as shown above, these errors are predicted to be small for a well-segmented image volume and can be easily identified by visual inspection and corrected if needed using the Origami pipeline.
Additionally, orientation flipping errors were quantified from real light-sheet fluorescence microscopy data from structures in the developing zebrafish otic vesicle (Figs 1 and 3). For this, cells assigned the wrong orientation along the apico-basal axis-that is, facing the basal surface instead of the apical surface, were identified by visual assessment in the Origami pipeline, showing errors in 3.65% of n = 949 cells analysed (Fig 2D).

Proof of principle: Insights into cell shape dynamics during epithelial morphogenesis within the zebrafish inner ear
To further validate our method, we used Origami to characterise cell shape dynamics involved in the formation of different structures in the otic vesicle of the zebrafish embryo (Figs 1 and  3). We analysed light-sheet fluorescence image data from the anterior epithelial projection (AP) for the developing semicircular canal system, together with the endolymphatic sac (ES), at three developmental time points: 42.5 hours post fertilisation (hpf) (time point 1), 44.5 hpf (time point 2) and 50.5 hpf (time point 3), using three different fish for each time point. We also analysed the posterior epithelial projection (PP), a similar structure to the AP, but which develops later [40], at developmentally equivalent time points to that of the AP (46.5 hpf, 50.5 hpf and 60.5 hpf). The AP and PP are finger-like projections of the epithelium that move into the lumen of the vesicle, with the apical side of the cell on the outside of the curved projection surface [40]. By contrast, the ES forms as an invagination from dorsal otic epithelium, with the constricted apical surfaces of the cells lining the narrow lumen of the resultant short duct [8,45,46]. As the ES is formed through deformation of the epithelial sheet with opposite polarity to that of the epithelial projections, we expect cells in the ES and the projections to show significant differences in cell shape. Conversely, we do not expect significant differences in cell shape between the AP and PP cells, which form equivalent structures in the developing ear.
For each structure, the following shape attributes were computed at the single-cell level: surface area, sphericity, longitudinal spread, transversal spread and skewness. Since volume and surface area show high collinearity within our data (Pearson correlation coefficient = 0.98, 95% confidence intervals = [0.977, 0.984]), cell volume was excluded from further analysis. Although images included cells in the non-folding epithelium around the developing structures of interest, only cells from the folding epithelium were analysed. A multivariate analysis ({MANOVA.RM} package [47] (Table 1 and Fig 3).
The differences in surface area are likely to be attributed to differences in sphericity between the cells in the three structures, but not in dimensions, as the transversal and longitudinal spread showed no significant differences.

Availability and future directions
Origami is free to download from: https://github.com/cistib/origami. It is implemented within MATLAB (compatibility with version 2018b onwards) and includes additional tools for visualising cell shape metrics from complex folding epithelia at the single-cell level. Instructions for installation and use are included with the software.
Our software can accept pre-segmented data, making it compatible with segmentation algorithms of the user's choice, potentially allowing for data acquired using other 3D imaging techniques such as tomography to be analysed. Segmented data must represent cell shape accurately, and so the choice of imaging technique that can faithfully detect 3D cell shape alongside membrane or cytoplasm-based segmentation is critical.
We used a priori knowledge of the otic epithelium organisation to inform the orientation of the apico-basal axis of the epithelial sheet to face the lumen of the otic vesicle [2,8,40]. It is essential to know the organisation of the apico-basal axis of cells within any new structure studied to apply Origami-wherein, the orientation-determining parameter can then be set to always be 'in' or 'out' depending on if the polarity direction vector is required to point towards the inside or outside face of a curved structure. We also assumed that individual cells do not violate this organisation, as this cannot be detected without additional polarity-specific labels. In such a case, polarity vectors from our analysis can be complemented with information from polarity-specific labelling to track such behaviour. Moreover, to compute shape features along an alternative axis of polarity, the pipeline can accept pre-assigned polarity as a cell-specific vector-list.
We expect Origami to be applied to studying a wide range of morphogenetic processes and to contribute to our understanding of the biomechanical processes underpinning them. Supporting information S1 Software Code. Zip file containing MATLAB scripts and instructions for installing and running Origami software. Requires MATLAB (v 2018b onwards). (7Z) S1 Text. Text file detailing methodology used for collection of fluorescence microscopy data, generation of synthetic membranes and parameters used for membrane segmentation. (DOCX)