MorphoGraphX 2.0: Providing context for biological image analysis with positional information

Positional information is a central concept in developmental biology. In developing organs, positional information can be idealized as a local coordinate system that arises from morphogen gradients controlled by organizers at key locations. This offers a plausible mechanism for the integration of the molecular networks operating in individual cells into the spatially-coordinated multicellular responses necessary for the organization of emergent forms. Understanding how positional cues guide morphogenesis requires the quantification of gene expression and growth dynamics in the context of their underlying coordinate systems. Here we present recent advances in the MorphoGraphX software (Barbier de Reuille et al. eLife 2015;4:e05864) that implement a generalized framework to annotate developing organs with local coordinate systems. These coordinate systems introduce an organ-centric spatial context to microscopy data, allowing gene expression and growth to be quantified and compared in the context of the positional information thought to control them.


Introduction
Many aspect of animal morphogenesis are thought to be controlled by positional information (Wolpert, 1969), where cells can sense their position in a developing organ and respond accordingly. This phenomenon may be even more pervasive in plants, as cells cannot relocate within organs, and must decide their fate based on their location. For example, root morphogenesis appears to be controlled by an organizing center at the root tip that provides founder cells and positional information to the growing structure (Scheres et al., 2002). Ablation of cortical cell initials in the root meristem causes the neighboring pericycle cells to divide and fill the available space, subsequently adopting the fate associated to their new location (van den Berg et al., 1995). A similar effect his has been demonstrated for a variety of cell types in the Arabidopsis root (Marhava et al., 2019). In leaves, development is thought to be coordinated by polarity fields oriented from leaf base to tip Kuchen et al., 2012). Over time organs can initiate new growth axes, such as when serrations or leaflets develop in more complex leaves (Barkoulas et al., 2008;Kierzkowski et al., 2019), or lateral roots emerge from the primary root (Scheres et al., 2002). In these cases information from several organizers must be integrated to direct cell response.
To understand how positional information controls morphogenesis, it is necessary to quantify cell shape, gene expression and morphogen concentration changes over time, preferably at the cellular level. This information then needs to be related to its position relative to the organizers controlling development within the organ. As computational power and imaging methods improve, new software packages for cell segmentation and lineage tracking are being developed (Sommer et al., 2011;Stegmaier et al., 2016), including many specialized for plants (Barbier de Reuille et al., 2015;Eschweiler et al., 2019;Fernandez et al., 2010;Schmidt et al., 2014;Wolny et al., 2020). This progress has enabled the segmentation of time-lapse data at increasingly higher resolution and throughput (Hervieux et al., 2016;Kierzkowski et al., 2019;Sapala et al., 2018;Willis et al., 2016)). Although this increase in data volume offers tremendous potential to understand how genes control form, the analysis of geometric data from thousands of cells is non trivial. Information about a cell's shape, gene expression and growth directions is of limited value when the cell's spatial context within the developing organ is unknown.
MorphoGraphX is a computer software platform that is specialized for image processing on surface layers of cells (Barbier de Reuille et al., 2015). It has proven especially useful for the analysis of confocal microscopy images from time-lapse data in order to quantify the cellular level dynamics of growth, cell division and gene expression (e.g. Bringmann & Bergmann, 2017;Feng et al., 2018;Hervieux et al., 2016;Hong et al., 2016;Kierzkowski et al., 2019;Louveaux et al., 2016;Sapala et al., 2018;Scheuring et al., 2016;Tsugawa et al., 2017;Vlad et al., 2014;Zhang et al., 2020;Zhu et al., 2020). Key to the approach taken in the software is the representation of cell layers as curved, triangulated surface meshes that capture the overall 3D shape of organs, which retains much of the simplicity of 2D segmentation and lineage tracking. These "2.5D" images contain the geometry of the sample at two scales. The global shape of the organ is captured by the mesh's geometry, while a cellular-scale representation is obtained from the confocal signal projected onto the mesh, which is segmented to extract the shape of individual cells on the surface ( Fig. 1A-C). When combined with time-lapse data acquisition and cell lineage tracking, MorphoGraphX allows cell growth and its relationship to gene expression to be quantified calculated (Fig. 1D,E; Kierzkowski et al., 2019;Sapala et al., 2018;Vlad et al., 2014). In addition to cell surface analysis, MorphoGraphX also supports the creation and analysis of full 3D meshes with volumetric cells (Fig. 1S1, Vijayan et al., 2021). Here we describe new methods we have developed in MorphoGraphX to understand these data by additionally annotating cells with positional information. Not unlike the annotation of sequence data, this allows cellular data to be given spatial context, and a frame of reference within the organ relative to its developmental axes and the organizers instructing morphogenesis.

Defining directions within an organ
The simplest method to provide positional information for the cells in a sample is by aligning the sample with a set of 3D coordinate axes (Fig 2A). For example, a developing root meristem can be aligned and positioned such that the organizing quiescent center is at the origin with the Y-axis increasing in the longitudinal direction of the root. Provided the sample is reasonably straight, this allows cellular measures to be compared with their distance from the quiescent center ( Fig  2B) (Schmidt et al., 2014).
However, for highly curved organs significant errors will occur, especially in more distal regions, further from the origin. It is possible to overcome this problem by placing a curve along the central axis (Fig 2C, Montenegro-Johnson et al., 2015). For this curve, MorphoGraphX uses Bezier splines defined by control points. These points can be positioned to create an axis that conforms to the curvature of the organ, using either interactive manipulation of the control points, or automatically from a selected file of cells. Distance can then be calculated along the line, and transferred to cells in the cross section perpendicular to the line. MorphoGraphX also allows a 2D Bezier surface to be positioned next to or within a sample, enabling two directions to be aligned with the natural curvature of the sample.
An alternative method applicable to curved organs with more complex shape, is to select one or more cells at a reference position, and calculate the distance relative to the selection ( Fig 2D, Movie 1). This offers an easy method to create a distance field, and greatly increases the variety of organs that can be accommodated. The distance is determined by computing the shortest path along cells through the tissue, causing it to naturally follow the curvature of the organ. As many organs have a layered cellular organization, a distance field normal to the sample's surface can be calculated. The surface of the sample is extracted and is used to annotate a full 3D segmentation of the interior cells. Fig. 3A shows a shoot meristem with the cells colored by distance to the surface, which can be used to classify cells ( Fig  3C). A manually defined Bezier curve (white) allows the assignment of accurate cell coordinates along a curved organ axis. (D) Side and top view of an A. thaliana sepal with a proximal-distal (PD) axis heat coloring. The cell coordinates were assigned by computing the distance to manually selected cells (outlined in red) at the organ base. This method allows organ coordinates to be assigned in highly curved tissues. (E) Side and top view of (D) with a heat map coloring based on cellular growth to the next time point. (F) Plot summarizing the growth data of (E) using the PD-axis coordinates from (D). See Figure 2S1 for the analysis of the complete time-lapse series. Scale bars: (A) 20 μm; (C) 100 μm; (D, E) 50 μm.  Fig 2D). The distance of the maximum of the growth zone from the base or the organ is relatively constant. However, organ length is increases about 10x between the first and last time points, making a comparison of the different curves difficult. (D-F) When plotting the same data with normalized cell distance values averaged using 20 bins along the PD-axis it becomes more apparent that the growth zone moves from the proximal to the distal regions over the course of development (D). The trend of lower and more proximal maxima (highlighted with arrows) is even clearer when proliferation is plotted in the same way (E). (F) Cell area data plotted as in (D) and (E). Average cell areas increase mainly at the distal end during later time points. Scale bars: (A, B) 100 μm.

Combining directions
In 2D or on 2.5D surfaces, local directions can be fully defined by a single distance measure, by taking one direction aligned with the gradient of the distance field or a Bezier curve, and the other perpendicular to the first. This is similar to methods used to specify directions in developmental modeling in plants (Green et al., 2010;Kennaway & Coen, 2019;Kierzkowski et al., 2019;Kuchen et al., 2012;Whitewoods et al., 2020), and thus facilitates direct comparison between models and experimentally observed patterns of growth and gene expression.
In 3D, a third direction must be defined (Kennaway & Coen, 2019;Whitewoods et al., 2020). In MorphoGraphX this can be done by combining the directions defined by different distance measures. The 3D Cell Atlas add-on for MorphoGraphX (Montenegro-Johnson et al., 2015) combines several distance measures for radially symmetric structures such as root and hypocotyls. A Bezier curve is placed along the center in the longitudinal direction and combined with a surface mesh to obtain radial directions ( Figure 3B). This also puts bounds on the radial direction (and also implicitly on the circumferential direction) which allows relative coordinates to be assigned to cells in addition to absolute values. For less regular organs, like the ovule, coordinates derived from a surface mesh (Fig 4A) can be combined with a Bezier curve along the surface (Fig 4C) to define a coordinate system for the region of interest in full 3D. In both of these cases, relative coordinates facilitate the classification of cells into layers. When annotating a root using the 3D Cell Atlas coordinate system, the relative radial coordinate will follow the layer as the organ narrows towards the tip. The use of relative coordinates also makes it possible to pool data from multiple samples (Vijayan et al., 2021;Zhang et al., 2020) and to compare data from different genotypes Montenegro-Johnson et al., 2019). For volumetric tissues often a single direction is not enough to capture the geometry of the organ. Different methods can be combined such as a Bezier curve (white dashed line) with a surface mesh (grey) to create a heat map of the relative radial distance of cells in the A. thaliana root. (C-D) Organ coordinates can be used to assign cell type labels as demonstrated in the 3D Cell Atlas plugin for meristem and root. See also Figure 3S1. (E-H) Different methods to create cell type labellings. (E) A. thaliana gynoecium (fruit epidermis) surface segmentation with a heat map of the length of the minor cell axis as obtained from a PCA on the cells' triangles. The heat values can be thresholded to assign two cell types. (F) The same principle can be used on organ coordinates which results in a clean separation of replum (green) and valve tissue (blue). (G) We generalized the 2D clustering approach of 3D Cell Atlas (see Figure 3S1) so that it can be used for any measure pair and on subset selections of cells. Shown is a 2D plot of the minor axis length (xcoord) and cell signal intensity (y-coord) on the valve tissue in (F). Manually assigning clusters can separate the stomata, which are typically smaller with higher signal values (yellow) and the remaining valve cells (blue) efficiently. See also Figure 3S1

Using positional information for data analysis
Once cells have been annotated with positional information, it can be used to analyze cell-level data, such as growth, cell proliferation, cell shape and gene expression. Using the distance measure to define the proximal-distal axis (Fig 2D), geometric measures can be plotted against the local coordinate system. In Fig 2F cell area extension was plotted against distance from the base of the sepal. On the full 7 day sepal time-lapse shown in Fig 2S1 ( Hervieux et al., 2016), it can be seen that initially cell growth is more distal, with a band of high growth progressing towards the base of the sepal. By time point 6, the growth has slowed and become more uniform as the organ differentiates. Proliferation is initially more uniform, but otherwise follows a pattern similar to growth, progressing basally as the organ matures. The data can be indexed by position and visualized in graphs, showing how growth and proliferation vary along a developmental axis ( Fig. 2S1C-F). This also makes it possible to group data collected from multiple samples and to compare different genotypes Zhang et al., 2020).
In addition to scalar information such as areal growth rate or cell volume, MorphoGraphX can also quantify directional information, such as the Principal Directions of Growth (PDGs) that represent the maximal and minimal directions of growth for each cell ( Figure 5A-C). A common problem with the interpretation of such directional information is the tendency for directions to be locally heterogeneous when growth is nearly isotropic. This happens because the maximal and minimal growth amounts are almost the same, and the displayed directions become arbitrary, and heavily influenced by noise. This can make the comparison of growth directions between neighboring cells difficult. A more informative approach is to look at growth with respect to the directions of the developmental axes. This can be done by first setting up an axis defining the positional information for the leaf, for example by using a distance field ( Fig 5B). The growth directions are then projected onto this developmental axis, and separated into components that are parallel (Fig. 5D) and perpendicular ( Fig 5E) to the axis. In Fig 5D a gradient of growth along the proximal-distal axis can be seen, with an increase in lateral growth around the forming serration ( Fig 5E). This is not immediately apparent in the original PDG visualization ( Figure 5C). Another benefit of looking at PDGs in the context of a local coordinate system is that it can provide a more direct comparison to the outputs of computational simulations. Developmental models of emergent organ shape often use morphogens that are thought to specifically control growth in the different directions in relation to a developmental axis Kuchen et al., 2012;Whitewoods et  al., 2020). By projecting the PDGs onto this axis, it is possibly to directly compare model growth rates in the different directions to experiments. Since MorphoGraphX can load a wide variety of mesh formats, this allows the direct comparison of similar quantifications made on templates extracted from model simulations from various sources. Figure 5F-K shows a similar growth quantification for the tomato meristem, where suitable local organ coordinates were created using cell distance measures around each emerging leaf primordium with directions pointing towards (radial) and around (circumferential) their respective center. Additionally, the signal intensity of the auxin reporter DR5 was quantified in the same sample. For both primordia we found radial growth to have a high negative correlation with DR5 signal intensity, whereas circumferential growth was more or less constant. The DR5 maximum tends to be on the adaxial side of the initiating leaf, whereas growth is much higher on the opposing abaxial side. This suggests that auxin acts as a trigger for primordium initiation, rather than via controlling growth rates directly.
The mature ovule in Arabidopsis shows a more complicated structure than a root or sepal, with four layers of integument cells encapsulating the nucellus that contains the embryo sac (Schneitz et al., 1995;Vijayan et al., 2021, Fig1S1, Fig4). After segmentation and 3D mesh extraction in MorphoGraphX, directions normal to the surface ( Fig. 4A) were combined with a Bezier curve computed from a userselected cell file, and used to construct a coordinate system for the outermost cell layer (Fig. 4B,C). The Bezier curve defined the longitudinal direction whereas the surface directions obtained from the organ surface mesh were used to compute perpendicular width and depth axes and distances. Cell volume and geometry acquired from 3D segmentation and mesh extraction ( Fig 1S1) were calculated along the various directions of the organ axes and analyzed ( Fig. 4D-G).
Measuring the length, width and depth of cells along the cell layer and surface directions revealed the underlying cause for differences in cell volume between different proximal-distal regions of the outermost integuement layer. Moving along the proximal-distal axis, we found variations in cell volume with a clear minimum at around 100 μm and a steady increase towards the distal end (Fig. 4G). Cell length at the proximal and distal ends showed substantial differences, not seen in cell width and cell depth. The increase in cell length at the distal end was accompanied by enhanced cell anisotropy. These findings suggest that the overall steady rise in cell volume at the distal end is mainly due to differential cell length (Fig. 4G).  (G). (H) To analyze the data we defined a circumferential coordinate system with its axes directions (white lines) around the primordium and initiation center (not shown), and aligned them towards the meristem center. (I) Heat maps of cell distance, area extension, radial and circumferential growth and normalized DR5 signal intensity of the aligned primordium and initiation site. (J) Plotting the data of (I) reveals a negative correlation of the DR5 signal intensity and radial growth around the developing primordium. (K) Detailed plots of radial (red) and circumferential growth (orange) as well as the normalized DR5 signal intensity of the primordium and initiation site. Scale bars: (A) 50 μm; (B -H) 20 μm; (I) 10 μm.

Using positional information for automatic cell type classification
Plant organs typically emerge as primordia consisting of undifferentiated tissue. Cells subsequently differentiate, acquiring a unique identity that depends on their location within the organ, via genetic processes that integrate spatial and environmental cues. Although cell differentiation is ultimately controlled by differential gene expression, it is often the case that cell fate can be predicted by geometry, even at very early stages (Yoshida et al., 2014). It is rare that cells with different cell type have identical morphological features. Although most workflows in MorphoGraphX begin by segmenting 3D image stacks into cells, this is primarily an initial step to enable further analysis. As such, recent method development for MorphoGraphX has largely focused on the downstream processing of segmented data. The software now supports a large variety of measures to quantify different features of cell morphology, including: simple geometric measures (area, volume, perimeter, surface area, min and max axis), shape quantifiers (convexity, circularity, lobeyness, largest empty circle, aspect ratio), neighborhood measures (number of neighbors, variability), gene expression (average, total, near boundary), and cell network measures (betweenness centrality, betweenness current flow). Most measures can be used on time-lapse data to quantify changes over time (growth rates, gene expression changes, cell proliferation). For a complete list of the measures implemented in MorphoGraphX, see Supplemental tables 1, 2. The modular architecture of MorphoGraphX also allows custom measures tailored to specific problems to be easily added through its plug-in interface. More sophisticated calculations, for example the averaging of data over multiple samples, can be calculated externally in packages such as R and imported back into MorphoGraphX for visualization on segmented meshes. Here, the development of more complex data-flows is enabled by the use of a standardized attribute system to store and visualize cell data, for both scalar values and tensor (directional) information.
During the segmentation process MorphoGraphX assigns a unique label to each cell. A secondary cell label is provided (parent label), which is often used for lineage tracking. Other secondary labelings are also possible, for example cell type, cell layer, or zones within an organ. MorphoGraphX has several methods to assign these labels to classify cell types and layers. These labels can be assigned manually by interactively selecting cells, or by employing a number of processes that use heat map measure data to assign secondary labels ( Since cells of a common biological cell type have similarities in one or more geometrical, positional or gene expression attributes, the values of these attributes will often form a cluster, facilitating their automatic classification. An example can be seen in the 3D Cell Atlas addon for MorphoGraphX (Fig. 3S1 A-C; T. D. Montenegro-Johnson et al., 2015) that clusters cells by relative radial distance and cell size to classify the cell layers of the root, hypocotyl, mature embryo, or other radially symmetric plant organs. To aid in optimizing cell clustering, MorphoGraphX offers a two dimensional interactive heat map, where information from two independent measures can be visualized, and clusters selected ( Figure  3S1 B, 3G). These methods can be used repeatedly on sub-sets of cells to enable a classification of cells that differ across more than 2 features.
Multi-feature classifications tasks can be solved automatically by machine learning approaches (Cortes & Vapnik, 1995) when provided with sufficient training data. Of particular relevance are Support Vector Machines (SVMs) which have been used to classify cell types based on geometrical features of plant cells (de Reuille & Ragni, 2017;Sankar et al., 2014) . MorphoGraphX provides a simple interface to the libSVM support vector machines library (Chang & Lin, 2011). Cells can be selected and classified into different cell types for use as training data (Fig. 3H, Movie 2). Any cell attribute or measure that can be quantified in MorphoGraphX can be used by the classifier. These include all the morphological and gene expression measures, time-lapse measures, the positional information created via distance maps or other coordinate systems, and even custom measures created via plugins or calculated externally with R or MATLAB. Once trained on a small group of cells with the desired measures, the classifier can be used to classify all the cells in a sample (Fig. 3H). After manual curation, the classification can then be used as additional training data, improving the power of the classifier. Figure  3G,H and 3S1 F-G shows the cell types of the Arabidopsis gynoecium, with stomata homogeneously distributed within the valve, consistent with the uniform growth and differentiation of this tissue (Eldridge et al., 2016;Ripoll et al., 2019).

Mapping positional information through time
In the analysis of morphogenesis, many key quantifications such as growth depend on the ability to track samples through time. In MorphoGraphX this can be done following cell segmentation by manually assigning parent labels to the second time point, a process that has been highly streamlined in the user interface for 2.5D surfaces. However for full 3D samples, or large 2.5D samples with multiple time points, this method can be cumbersome. One method to address this problem is to find a non-linear coordinate transformation or deformation that maps all the points from one time point onto the next. Parent labeling can then be determined by mapping cell centers of the later time point to an earlier one, allowing the cell they came from to be identified. This can be used to directly assign lineage, or to seed algorithms that use more involved methods, such as the minimization of the total distances between the mapped cells (Fernandez et al., 2010). In MorphoGraphX, to define a mapping between the meshes of two successive time points (Fig 6A) we have implemented a 3D deformation function based on scattered data point interpolation (using cubic radial-basis functions; Duchon, 1977;Turk & O'Brien, 1999). An initial transformation is computed based on a few preassigned landmarks, by matching several cells with their parents in the previous time point (Fig 6C). A deformation field is then calculated which provides a mapping for all points in 3D. This is then used to assign parent labels in the second time point based on their closest match in the previous time point. Close to the landmark points, this mapping will be very accurate, with accuracy decreasing with distance from the landmarks. The decrease in accuracy away from landmarks is larger if the deformation between the time points is highly non-uniform. After assigning all the cells their closest parent, the mapping is then verified by checking the correspondence of neighborhoods between each cell and its parent. The labeling for cells which do not match is cleared, and the process is repeated. This causes correctly labeled regions to "grow" out from the initially placed landmarks, until the entire sample is correctly labeled (Fig 6C, Movie 3). At each step, only correctly labeled cells remain. Sometimes the iterative cell-labeling process can get stuck in highly proliferative areas where cells have divided repeatedly between time points. In this case a few additional landmarks can be manually added at trouble spots. One significant advantage of the method is that incorrect cells remain unlabeled, making manual curation straightforward. Once all of the parents are assigned and have passed the neighborhood correspondence check, one can be assured that both the lineage and the underlying segmentations are correct.
Deformation functions can also be used to create animations of organ development from 2.5D or 3D time-lapse data. This requires two or more time points of a segmented mesh with corresponding cell lineages. The cell centers and/or junctions are used as the landmarks defining the deformation function that maps one mesh onto the other. Interpolating mesh vertices between stages creates a smooth animation with as many intermediary steps as desired (Fig. 6S1, Movie 4). MorphoGraphX has a user-friendly pipeline to record animations directly from the GUI with options to adjust the camera angle and to visualize cell lineages, heat maps and cell outlines during the animation. Temporal smoothing of morphing animations created from more than two time points is achieved using Catmull-Rom splines to interpolate the position of mesh vertices over time (Catmull & Rom, 1974). Heat and signal values in the mesh, such as cell area, growth rates or gene expression can also be interpolated along with vertex positions.
In large cells, growth can vary significantly within the same cell (Armour et al., 2015;Elsner et al., 2018). As the deformation function provides a smooth mapping between time-points, its gradient can be used to create a continuous growth map at any point on a mesh. This enables the approximation of areal expansion and PDGs at a sub-cellular level, where the quality of the approximation is limited by the number and placement of landmarks (junctions). It is also possible to apply the process to subcellular landmarks, such as those obtained by tracking microbeads, as done previously for 2D images (Armour et al., 2015;Elsner et al., 2018). Our 3D implementation of this method has been used to compute growth directions on curved surface meshes ( Fig 6E) and volumetric meshes (Fig 7I).
A comparison of the areal growth and PDGs calculated with deformation functions vs the cell-based method is shown in Figure 6D, E. The deformation function captures differences in growth within single cells, as is often apparent in larger cells that straddle areas of fast and slow growth. Figure 7A-I shows a 3D timelapse of the Arabidopsis root where the deformation functions have been used to perform lineage tracking in 3D. It can be seen in Fig 7D that the 4 tissue types, epidermis, cortex, endodermis and pericyle all show a similar growth pattern, with slow growth in the meristem and faster growth near the transition zone. When the growth is displayed as a function of the distance from the root tip (Fig 7 E-H), again it can be seen that the layers are almost the same, reflecting the almost completely anisotropic growth of this system. . P cells showed a larger volume (E), which was caused by a greater cell length (F), an observation which has been made before by (Andersen et al., 2018). In contrast, X and E cells were smaller in volume due to different reasons: While X cells were the shortest (E), E cells showed a lower cell width with increasing distance from the QC (G). The time-lapse analysis confirmed above observations: While volume change was similar across the cell type (I), P cells showed a lower proliferation rate (H), whereas E cells showed the smallest extension of cell width (J).Scale bars: (A, B) 20 μm; (C, D) 50 μm.

Advanced geometric analysis
While MorphoGraphX was created to work with 2.5D surface projections, it now supports a complete set of tools for full 3D image processing, and in many samples advantages can be gained from combining both techniques. Additional tools, such as the automated cell lineage tracking on surfaces, have also been extended to 3D to facilitate growth analysis in full 3D. This is a much harder task than the analysis of surface images, as 3D cellular meshes lack the relatively easy to identify junctions that serve as material points for surfaces. However in many cases entire organs are well defined by their surfaces meshes, allowing landmarks on the surface to be used to construct a 3D deformation function to aid lineage tracking in full 3D. Surface landmarks can also be combined with 3D cell centers and/or the face centers as material points to improve the internal resolution of the deformation functions for full 3D samples. These techniques allow the methods used to calculate growth rates and PDGs in 2.5D to be extended to full 3D (Fig. 7I). Cell proliferation and most of the other measures can also be quantified from 3D timelapse data (Fig 7S1). In addition to the automated tools, improved manual 3D parent labeling and the ability to relabel cells so that adjacent cells are always a different color, aid in the manual curation of 3D lineage maps.
One of the more advanced quantifications from time-lapse data is the analysis of cell division. As plant cells cannot move, cell division and growth are the main determinants of morphogenesis. MorphoGraphX has processes to identify dividing cells from time-lapse data, and quantify the orientation of the division wall in both 2.5D and 3D (Fig 8A-I, Fig. 8S1). In 2.5D the best fit line to the division wall is calculated (Fig 8A,D), whereas in 3D the best fit plane is used (Fig 8B,C,G). There are also measures to determine the asymmetry of the daughter cells. The use of positional information to give organ context is even more important for directional information. For instance, quantifying the orientation of the division plane is of little use without knowing how it relates to the developmental axes. If the organ shape is simple and aligned with one of the axes, then cell division orientations can be used directly (Fig 8S1A). When more complex shapes are involved, orientations can be computed with respect to the axes of a local coordinate system defined for the organ, along with its associated positional information (Fig. 8E,I). It is also possible to quantify how close cell divisions are to common division rules proposed in the literature, such as the shortest wall through the center of the cell including local minima (Fig 8H; Besson & Dumais, 2011), along the principal directions of growth (Hejnowicz, 1984), or rules based on network measures . The organization of cells in organs may be analyzed through the extraction of cell connectivity networks from 2.5D or 3D segmented data. The physical associations between cells (cell-cell wall areas) can be extracted and converted into networks where they are analyzed using network measures (Fig. 8J-K). Local measures such as the number of immediate neighbors (degree) can be calculated, along with more global measures, such as betweenness centrality based on the number of shortest paths cells lie upon, or random walk centrality (Fig. 8K). These global measures are central to understanding how information flows within tissues (Jackson et al., 2017. The use of these measures uncovered the presence of a global property in cellular organization within the Arabidopsis SAM . Namely, the length of paths between cells is maximized, whereby cells which lie upon more shortest paths have a great propensity to divide, and the orientation of this division tends to leave the two daughter cells on the fewest number of shortest paths. Using this approach, the local geometric properties of cells can be related to the emergent global organization of cellular arrangements. Perturbation of cell shape in the katanin1 mutant led to alterations in path length in the SAM, which correlated with defects in phyllotaxis . Reporter signals, such as proteins tagged with fluorescent proteins, can also be quantified in MorphoGraphX. After segmenting a surface into cells by using a cell wall stain or marker line, a signal collected in a second channel can be projected onto the surface mesh, and the abundance, orientation and polarity of signals can be computed. Examples are the PIN-FORMED (PIN) auxin transporter report line (Benková et al., 2003) or the GFP:MBD (Van Bruaene et al., 2004) line that tags microtubules (MTs), ( Figure 8L-N, 8S2). For the quantification of cell polarity the projected signal along the cell border is binned based on its position in relation to the the cell center to obtain its predominant direction and its intensity. Figure 8L shows an example of the PIN1 polarity quantification at the cell wall of surface segmented cells in the SAM. A similar quantification can be performed for 3D meshes as shown in Figure 8N and Figure 8S2A,B where we computed the PIN2 polarity in epidermis and cortex cells of an Arabidopsis root. Again, this directional information can be combined with the organ coordinates to compute the angle between cell polarity and the organ axis ( Figure 8S2C). Another example is the quantification of MT alignment using an implementation of Fibril tool (Boudaoud et al., 2014) that has been adapted for processing on surfaces. After projecting the MT signal onto the surface, the alignment direction and strength of the signal can be quantified at the sub-cellular level ( Figure 8M) or for entire cells ( Figure 8S2E). Again, this information can be interpreted using organ coordinates, as we demonstrate on cells of a SAM which tend to have their MTs aligned circumferentially from the meristem center ( Figure 8S2E).   Figure 7C) that was used for the division plane analysis in Figure 8E. Cells are shown semi-transparently (grey) with their longitudinal organ axis (cyan) obtained from the analysis using 3D Cell Atlas in Figure 4B. Planar approximations of the division planes between cells that divided between the 2 time points are shown as red circles. Consistent with the quantitative analysis in Figure 8E

3D visualization and interactive tools
MorphoGraphX has a flexible rendering engine that can handle meshes containing millions of triangles. It supports the independent rotation and translation of different stacks and meshes in the same world space and the ability to render both voxel and geometric data together with blending and transparency. It has adjustable clipping plane pairs and a bendable Bezier cutting surface that can be used to look inside 3D samples, and an interface to support the creation of animations (Movie 5). However, visualizing and interacting with 3D data on a 2D computer screen still remains a challenge. A particular problem is the validation and correction of 3D segmentations of organs, as internal cells are obscured by outer layers. Exploded views are a commonly used method to visualize the internal structure of multi-component 3D objects (Fig 8S1 C;Li et al., 2008), which can be used on mesh data in MorphoGraphX. To add biological meaning to the exploded visualization, cells can be bundled by their parent or cell type label to visualize key aspects of biological data sets such as cell divisions or to separate organs into cell layers or by cell type (Fig 7D, 8S1C). These processes can also aid user interactions, making cell selection and annotation more straightforward when users are curating 3D cellular segmentations and lineage maps.
In addition to mesh editing tools, MorphoGraphX has several tools to edit voxel data. The simplest is an eraser tool which can be used to remove portions of the stack that would otherwise interfere with processing. An example is the digital deletion of the peripodial membrane overlying the Drosophila wing disc, which needs to be removed to allow for the extraction of the organ's surface (Aegerter-Wilmsen et al., 2012). A typical workflow for 3D segmentation starts with a 3D image of a cell boundary marker. This is then pre-processed with operations such as blurring to reduce noise or background removal filters, before segmentation with algorithms such as the ITK (www.itk.org) morphological watershed. More recently, deep learning methods with convolutional neural networks (CNN) have been developed to predict cell boundaries, such as the 3D U-Net model (Çiçek et al., 2016), that can improve the stacks for downstream segmentation. The modular structure of MorphoGraphX has allowed us to interface an implementation of the 3D U-Net model developed by Eschweiler et al. (2019). This enables the interactive use of the CNN boundary prediction tool from within MorphoGraphX (Movie 6), and simplifies experimentation with different networks or downstream segmentation strategies. It also avoids the requirement to set up a full Python environment, and the associated installation issues. Although some tools have packaged 3D U-Net prediction with a choice of several segmentation algorithms  (Wolny et al., 2020), the tools are typically written as pipelines without intermediate visualization of the data, making it cumbersome to experiment with the different methods and parameters of the individual data processing steps in the pipeline.
Once the data is segmented, it often requires some manual correction before it is ready for final analysis in a chosen coordinate system. MorphoGraphX has interactive tools that operate on voxel data both to combine and split labels (cells), although typically it is easier to over-segment and combine, rather than undersegment and split (Movie 6). This can be used to correct segmentations, which can then be used to help train deep learning networks to further improve automatic segmentation (Vijayan et al., 2021). In this context MorphoGraphX has been used to segment and curate Arabidopsis ovule data to create ground truth for confocal prediction networks (Wolny et al., 2020).

Conclusion
Similar to sequencing data, geometric data on the shape and sizes of hundreds or thousands of cells is of limited value without annotation. For many developmental questions, the spatial context for information on cell shape, division and gene expression is paramount. However it is not enough to know the 3D position of cells, but rather their position in a coordinate system relative to the developing tissue or organism. These coordinates would typically reflect the developmental axes of the organism or tissue, allowing the direct comparison of cell and organ shape changes with the gene expression controlling their morphology. In addition to putting data in a mechanistic context, organ-centric relative coordinates can be used to compare samples with different morphologies (Thompson, 1942), such as different mutants, or even in different species . This also applies to changes in morphology over time, where organ coordinate systems can be used to determine the correspondence of cells at different stages of development. Several tools have been used to successfully harness organ-centric coordinates for specific problems (Montenegro-Johnson et al., 2015Schmidt et al., 2014), and MorphoGraphX provides a generalized framework to make coordinate systems customized to the particular organ or organism of interest.
In addition to annotation with organ-centric positional information, MorphoGraphX has a comprehensive toolkit for cell shape analysis, growth tracking, cell division, and the quantification of polarity markers, both on 2.5D image meshes, and for full 3D. All of these measures can be calculated and stored within the mesh, or exported to files for further processing with other software. Custom measures can also be calculated and imported for visualization within MorphoGraphX. Cell shape measures in combination with positional information provide a powerful framework for cell type classification, both with machine learning methods (Fig 3G-H)  A key strength of MorphoGraphX is that it offers both manual and automatic tools for segmentation, lineage tracking and data analysis. Although fully automated methods are improving, streamlined methods for manual and semi-automatic segmentation and analysis provide a path to completion for many samples where the automatic methods are "almost" good enough. For example, the automatic lineage tracking now available benefits from the streamlined tools we developed previously to do the process manually, as these are now used to correct and fill in missing portions when the automatic segmentation is incomplete. This reflects the interactive nature of MorphoGraphX and its focus on low-throughput but highquality data sets. Nevertheless, portions of a workflow can be fully automated by using python scripts to process many files. Furthermore, operations performed in MorphoGraphX are written to a python log file, allowing easy cut-paste script creation.
The internal architecture of MorphoGraphX has been designed to make it easily extendable, while retaining the speed of the fully compiled, statically typed language C++. The relatively small visualization and data management core is augmented with processes that are loaded dynamically at startup and provide almost all of the software's functionality. MorphoGraphX has grown to provide a wealth of custom processes for 2.5D image processing and coordinate system creation. Additionally, it has become a platform to integrate published tools and methods that have no visual interface to the data of their own, which increases their accessibility and ease of exploration for biologists. Examples include the previously mentioned CNN process, and the interface to the XPIWIT software (Bartschat et al., 2016) which provides a tool to develop ITK image processing pipelines interactively without any C++ programming required. These pipelines can be packaged into plugins and called directly from MorphoGraphX, allowing the exploration of different ITK filters. One XPIWIT pipeline we have prepackaged in a process is for the TWANG method (Stegmaier et al., 2014), a fast parallel algorithm for nuclear segmentation. Another example of software integration is the processes we have developed to interface with R (Movie 7) that provides plots for basic statistical analysis on attributes created in MorphoGraphX, including positional information provided by organ-centric coordinate systems. This simplifies the creation of the most commonly used plots without the need for export files.
As more and more imaging datasets are becoming available for community use, their annotation with positional information and gene expression data will be critical to help understand how the cell-level action of different genes and genetic networks is translated into the 3D forms of tissues and organs of different species. In this context MorphoGraphX provides a tool set to help maximize the attainable information from these datasets, in an accessible platform tailored to the experimental biologist.

Software Availability
MorphoGraphX is available for Linux and Windows, although we recommend Linux as some add-ons are only available for Linux, and some, such at the CNN add-on require Cuda. For Linux we provide a Cuda version for machines with a compatible nVidia graphics card, and a non-Cuda version for those without. Currently there is only a non-Cuda version for Windows. Although there is no Mac version, some have had success running it in a virtual machine. The software and documentation are available at www.MorphoGraphX.org.

Arabidopsis Mature Embryo (Fig2C)
Arabidopsis thaliana Col-0 seeds were sterilized in 70% ethanol with Tween20 for two minutes, replaced with 95% ethanol for 1 minute and left until dry. Seeds were placed on the Petri plates containing 1/2 MS medium including vitamins (at pH 5.6) with 1.5% agar and stratified at 4°C for 3 days in darkness. Next, seeds were imbibed for 3 hours, and the mature seed embryo was isolated from the seed coat. For live imaging, the embryos were stained with propidium iodide 0.1% (Sigma-Aldrich) for 3 minutes and imaged with Leica SP8 laser scanning confocal microscope with a water immersion objective (x20). Excitation wavelengths and emission windows were 535 nm and 617 nm. Confocal stacks were acquired at 1024x1024 resolution, with 0.5-μm distance in Z-dimension. Images were acquired at 48 hours intervals and samples were kept in a growth chamber under long-day condition (22°C, with 16 h of light per day) between imaging. From more than 10 replicates a sample with curved overall shape was selected for the demonstration of the organ coordinates using a Bezier curve. To quantify the cell area, change and anisotropy, the fluorescence signal was segmented, and cells were parent labeled manually between two successive time points. Heat-maps are displayed on the later time-point (after 48 hours of growth). Scale bars are displayed on the image.

Marchantia time-lapse (Movie 3)
Marchantia polymorpha gemmaling Cam-1 PM::GFP reporter line (Shani et al., 2010) were transferred on a Petri plate containing 1/2 Gamborg's B5 medium including vitamins (pH 5.5) with 1.2% agar and grown for 24 hours. For live imaging, the gemmaling were imaged with Leica SP5 laser scanning confocal microscope with a water immersion objective (x25/0.95). Excitation wavelengths and emission windows were 488 nm and 510 nm. Confocal stacks were acquired at 1024x1024 resolution, with 0.5-μm distance in Z-dimension. Images were acquired at 24 hours intervals and samples were kept in a growth chamber under constant light between imaging. For the move we selected a representative sample from 6 total replicates. To quantify the cell area, change and anisotropy, the fluorescence signal was segmented and semi-automated parent labelling as performed to couple  (Fig2D,E, Fig2S1A,B, Fig6D,E, Fig6S1B, Fig8A,D) Data previously published in (Hervieux et al., 2016). (Fig5A-E, Fig6S1A, Fig8J,K) Data previously published in .

Arabidopsis SAM (Fig3A,C, Fig3S1D)
Data previously published in (Montenegro-Johnson et al., 2019). (Fig 6B,C) pUBQ10::acyl:YFP has been described previously (Willis et al., 2016). Plants were cultivated on soil under the long-day condition (16 h light/ 8h dark), and 20°C ± 2°C. Flowers at post-anthesis stage from 5 weeks-old plants were dissected with fine tweezers to remove sepals and stamens to expose gynoecium and mounted on the 60 mm plastic dish filled with 1.5% agar. Confocal imaging was performed with a Zeiss LSM800 upright confocal microscope, equipped with a long workingdistance water immersion objective 40X (1 NA, Apochromat). Excitation was performed using a diode laser with 488 nm for YFP and the signal was collected between 500 and 600 nm. For both organs images of 3 replicates each were obtained. (Fig8B,C,F Data previously published in (Yoshida et al., 2014). (Fig8L,M, Fig8S2A-C) pUBQ10::acyl:TDT (Segonzac et al., 2012) and GFP:MBD (Van Bruaene et al., 2004) were crossed. F3 double homozygote line was used for imaging.

Arabidopsis Meristems for PIN1 and MT
Floral organs were removed with fine tweezers about 21 days after germination to expose inflorescence meristem. Meristems were mounted on the 60 mm plastic dish filled with 1.5% agar and imaged with a Zeiss LSM800 upright confocal microscope, equipped with a long working-distance water immersion objective 60X (1 NA, Apochromat). Excitation was performed using a diode laser with 488 nm for GFP and 561 nm for TDT. Signal was collected at 500-550 and 600-660 nm, respectively. Images of 3 replicates were obtained. (Fig8N-O, Fig8S2E) pPIN2::PIN2:GFP was previously described (Xu & Scheres, 2005). The seeds were stratified for 2 days at 4°C, grown on 1/2 Murashige and Skoog medium with 1% sucrose under the long-day condition (16 h light/ 8h dark) at 20C° ± 2°C. The roots were stained by 10mM propidium iodide (Sigma-Aldrich), and observed by Zeiss LSM780 with two-photon laser (excitation 990nm) with a band pass filter 500/550nm for GFP and 575-600nm for PI. Images of 3 replicates were obtained.

Data Analysis
For the data analysis examples in the paper we computed all necessary cellular data within MorphoGraphX and exported them as csv files. Those files were imported to RStudio for further processing or directly plotted using ggplot2 (R Core RStudio Team, 2020;Wickham, 2016).
Arabidopsis Flower meristem (Fig. 1) We selected one sample for segmentation and further analysis. The segmentation, cell lineages and heat maps were generated following the standard workflow as described in  and in the MGX user guide. Arabidopsis Ovule (Figures 1S1, Fig. 4) Segmentation was obtained by back blending the raw images to CNN boundary predictions (Wolny et al., 2020) as described in (Vijayan et al., 2021).
For the analysis in Fig. 4 we selected one sample of the published data.The trimmed surface mesh was used to extract the outermost layer using the method described in (Montenegro-Johnson et al., 2019). The organ axis was defined by a Bezier curve obtained from a manually selected central cell file using the processes "Misc/Bezier/Bezier From Cell File" and "Mesh/Cell Axis 3D/Custom/Create Bezier Grid Directions". Moreover, for each cell the direction towards the surface and the orthogonal direction of the Bezier and surface direction was computed ("Create Surface Direction"; "Create Orthogonal Direction"), resulting in 3 orthogonal organ axes. The cell sizes were quantified by first doing a PCA on the voxels of cells on the segmented stack ("Mesh/Cell Axis 3D/Shape Analysis /Compute Shape Analysis 3D") and finally computing the component of the PCA's tensor aligned with the axes of interest ("Mesh/Cell Axis 3D/Shape Analysis/Display Shape Axis 3D" with the appropriate "Custom" heat option), with the Bezier direction corresponding to cell length, the surface direction to depth and the orthogonal direction to width.
Shape Anisotropy was defined using the equation: (max -0.5*min -0.5*min)/(max+mid+min), with max, mid and min defined by the length of the PCA axes.
To create the plots, cells with a distance >40μm from the central cell file were ignored. The data of the remaining 213 cells was plotted. (Figures 2A, 4B,D From the 3 replicates we selected the sample with the best segmentation quality for further analysis. The two time points of the analyzed root data sample were segmented using ITK segmentation processes in MGX (see also  and the MGX user guide). From the segmented stack a surface mesh and volumetric cell mesh were obtained. For the axis alignment analysis (Fig 2A), the organ was manually aligned with the y-axis and the coordinates of the cell centroids were computed. In Fig 2B the  For the time-lapse analysis (Fig 7A-I, Fig 7S1) the cell lineages were determined semi-automatically using the pipeline introduced in this paper (Fig 6) followed by a manual error correction. PDGs in 3D were derived from the deformation function mapping the first onto the second time point using parent labelled cell centroids and cell wall centers.

Arabidopsis Root
For the analysis of the cell types in the endodermis (Fig 7S1), xylem cells in the stele and their neighboring pericycle cells were automatically identified by their circumferential coordinate. Endodermis cells touching two xylem-associated pericycle cells were determined as xylem pole. The two phloem poles in the endodermis were shifted by 90 degrees (or 2 cells) from the xylem pole. In total, we used 26 xylem pole, 21 phloem pole and 43 rest endodermis cells for the analysis.
For the analysis in Fig 8E and Fig 8S1A only the second time point was used. We computed the proliferation to the previous time point, extracted the vertices on each division plane between cells that have divided exactly once (proliferation = 2, n=249 mother cells that divided) and computed a PCA on each set of division plane vertices to extract the normals of the planes. Then we computed the angle between the longitudinal axis of the organ as extracted by 3D Cell Atlas and the division planes and exported to data. (Figures 2D,E, 2S1A,B, 6D,E, 6S1B, 8A,D)

Arabidopsis Sepal
For the sepal analysis one replicate of the data from (Hervieux et al., 2016) consisting of 7 time points was used (see Fig 2S1A,B). Fig 2D-F and Fig2S1 for each time point we manually determined the organ base based on the cell lineages from the first time point. Cells at the organ base were selected and used to compute the Euclidean cell distance measure. Finally, cell distances, growth, proliferation and cell sizes were exported.

For the analysis in
For the cell division analysis in Fig 8A,B we analyzed the divisions that occurred between the time point T4 and T5 (n=84). We computed the proliferation between these time points, extracted the vertices on each division plane between cells that have divided exactly once (proliferation = 2) and computed a PCA on each set of division plane vertices to extract the normals of the planes. Next, we computed the PD-axis direction of the organ using the Euclidean cell distance from the base. Finally we computed the angle between PD-axis and the division planes and exported to data.
For the growth analysis in Fig 6D we computed the PDGs from time point T4 to time point T5, visualized on the earlier time point. Fig 6E shows the same time point, but here growth was computed using the gradient of the deformation function obtained from the cells' junctions.

Arabidopsis Shoot Apical Meristem (Figures 3A,C, 3S1D)
We selected one sample of the published data for analysis. Cell type labels were determined using the methods described in 3D Cell Atlas Meristem (Montenegro-Johnson et al., 2019).

Arabidopsis Leaf (Figures 5A-E, 6S1A, 8J-K)
The Arabidopsis leaf data was previously published in . One replicate of a time-lapse series consisting of 7 time points was selected for analysis, but only time points T2 and T5 were used. The cell distance was computed similar to the Sepal example ( Fig 2D) as distance from the organ base. Additionally we computed the heat map gradient of the cell distance heat map ("Mesh/Cell Axis/Custom/Create Heatmap Directions") to obtain custom direction along the proximal distal (PD) axis and orthogonal to that the medial-lateral (ML) axis of the organ for each cell. PDGs were computed and used to determine the amount of growth along the previously computed PD and ML-axis.
For the cell network analysis in Fig 8J-K we computed the cell connectivity network of all cells in T5 weighted by the inverse of the length of the cell walls to determine the betweenness centrality ("Mesh/Heat Map/Measures/Network/Betweenness Centrality"), .

Tomato Shoot Apical Meristem (Figure 5F-I)
For the growth and DR5 signal analysis on the shoot apical meristem we used one replicate of the previously published data of (Kierzkowski et al., 2012). To objectively find the center of the meristem, the primoridium and the initiation site the curvature of the cells was computed. The resulting heat map was smoothed across neighboring cells for two rounds and resulting local maxima were identified as centers. Meshes were manually aligned along the x-axis with respect to the meristem center to compute circumferential coordinates ("Mesh/Heat Map/Measures/Location/Polar Coord") around primordium and initiation center. For the analysis only cells in the vicinity of the primordium and initiation centers were considered. These cells were determined using a threshold in the cell distance measure from those centers. Furthermore, the gradients of the Euclidean cell distance heat maps from both centers were used to compute custom directions along the heat (=radial) and orthogonal to the heat (=circumferential). Now the growth analysis was done similar to the leaf with computing PDGs and growth along the custom axis.

Arabidopsis Shoot Apical Meristems (Figures 8L,M, 8S2E,F)
For the MT analysis we selected one sample for segmentation and analysis. We determined the center of organ based on a smoothed curvature heat map. The center cell was selected and the Euclidean cell distance to the remaining cells was computed. The circumferential direction around the cell center was obtained from the orthogonal direction of the heat map directions. Cells were then binned by their Euclidean distance to the center.

Arabidopsis Embryo (Figures 8F,G, 8S1B-E)
The data for the 3D division analysis in A. thaliana embryos was previously published in (Yoshida et al., 2014). From this data set we chose one wild-type and one inducible bdl (pRPS5a>>bdl) sample at the 16 cell stage.
A surface mesh was generated from the cells in the embryo and the cells were parent labelled according to their predicted mother cell. Then the process "Mesh/Division Analysis/Analysis 3D/Division Analysis Multi" performed the following steps on all of the parent labelled cells (n=16 cells or 8 divisions in each genotype): First a planar approximation of the actual division plane was computed by performing a PCA on the vertex positions of the shared wall between the two daughter cells. Then 1000 equally distributed division planes were simulated on the combined mother cell and different measures were quantified. The actual and the best planes were visualized using "Mesh/Division Analysis/Display and Filter Planes". Figures 8N, 8S2A,B) For the analysis of the PIN directions in the A. thaliana root we selected one sample for segmentation. Next, we defined the main organ axis using a Bezier curve through the center of the organ. Then we computed the PIN2 polarity direction and determined the angle between the polarity direction and the Bezier line.