Abstract

It is a dynamic process that occurs at multiple spatial scales from macromolecular dynamics, to cell deformations, mitosis and apoptosis, to coordinated cell rearrangements that lead to global changes of tissue shape. Using time lapse imaging, it is possible to observe these events at a system level. However, to investigate morphogenetic events it is necessary to develop computational tools to extract quantitative information from the time lapse data. Toward this goal, we developed an image-based computational pipeline to preprocess, segment and track epithelial cells in 4D confocal microscopy data. The computational pipeline we developed, for the first time, detects the adherens junctions of epithelial cells in 3D, without the need to first detect cell nuclei. We accentuate and detect cell outlines in a series of steps, symbolically describe the cells and their connectivity, and employ this information to track the cells. We validated the performance of the pipeline for its ability to detect vertices and cell-cell contacts, track cells, and identify mitosis and apoptosis in surface epithelia of Drosophila imaginal discs. We demonstrate the utility of the pipeline to extract key quantitative features of cell behavior with which to elucidate the dynamics and biomechanical control of epithelial tissue morphogenesis. We have made our methods and data available as an open-source multiplatform software tool called 'I'I'T (http://github.com/morganrcu/TTT)

Author Summary

Understanding processes that make them acquire their final shape has implications to pathologies such as cancer progression and birth defects such as spina bifida. During development, epithelial tissues are remodeled by mechanical forces applied at the Adherens Junctions (AIS). The AIS form a belt-like structure below the apical surface that functions to both mechanically link epithelial cells and enable cells to remodel their shape and contacts with their neighbors. In order to study epithelial morphogenesis in a quantitative and systematic way, it is necessary to measure the changes in the shape of the AIs over time. To this end we have built a complete computational pipeline to process image volumes generated by laser scanning confocal microscopy of epithelial tissues where the AIs have been marked with A] proteins tagged with GFP. The system transforms input voxel intensity values into a symbolic description of the cells in the tissue, their connectivity and their temporal evolution, including the discovery of mitosis and apoptosis. As a proof of concept, we employed the data generated by our system to study aspects of morphogenesis of the Drosophila notum.

Introduction

The formation of structurally and functionally distinct embryonic structures requires that epithelial tissues change shape during development in a process called epithelial morphogenesis. A range of cellular behaviors drives these epithelial tissue shape changes, including cell shape change, rearrangements of cell-cell contacts, migration, proliferation, and programmed cell death. These behaviors, in turn, depend on intracellular molecular dynamics that allow cells to generate and transmit mechanical forces to one another, while maintaining epithelial cohesion [1, 2]. This dual requirement is fulfilled by the adherens junction (AI), a specialized protein complex that links epithelial cells together. The AIs form a planar belt-like structure below the apical surface of the epithelium composed primarily of the single pass adhesion protein E-cadherin (E-cad) and associated proteins. The extracellular domain of E-cad forms trans homo-dimers to promote cell adhesion. The intracellular domain of E-cad associates with the force-generating actomyosin cytoskeleton and functions as a site for the transmission of mechanical forces that can remodel cell-cell contacts and cell shape by influencing the dynamics of the AIs themselves [3].

Specifically, cell shape may be controlled by either autonomous or non-autonomous behaviors or forces, while the interaction between such local dynamics can lead to emergent effects on cell or tissue morphology. Live imaging of cell and molecular dynamics using fluorescently-tagged proteins is a key method to investigate these processes [4]. However, to fully leverage these experimental methods, quantitative approaches to automatically identify, track and interrelate molecular, cell and tissue level dynamics are required. The analysis of this quantitative information could then suggest molecular, cellular, and tissue level mechanisms that drive morphogenetic processes, and guide experimental approaches to test these possible mechanisms [5—9].

These methods are based on the detection of the AIs in projections of 3D information into 2D planes. These methods provide approximations of epithelial shape but often lead to inaccurate representations of cell shape, especially in curved regions of epithelial sheets. The projection of an image volume into a 2D plane also increases the image noise, which may interfere with image preprocessing and quality of segmentation. Therefore, to circumvent these problems and provide a more accurate representation of epithelial cells, methods to segment and track the AIs in 3D need to be developed.

These methods consider objects as particles moving independently without restrictions imposed by cell packing in an epithelial sheet, and thereby do not model cell shape and cell-cell contact dynamics. Only a few methods have been developed for the segmentation and tracking of cells in tissues. A common approach is to perform a dual analysis of multichannel time lapse images labeled for both nuclei and membrane. Locations of cell nuclei are first identified and later employed to find the location of the cell membrane. Using this approach Bourgine et al. [15, 16] have built a Partial Differential Equation framework to filter image noise, segment cell nuclei, locate the cell membrane and model the evolution of cell shapes to track the cells. A similar scheme has been developed by Luengo-Oroz et al. who employed a 4D morphological structuring element framework to denoise image volumes and locate cell nuclei, and a viscous watershed algorithm to infer cell nuclei locations as seeds for detecting cell extent properties [17]. This process is performed in 4D, implicitly tracking the cells. The same group has developed an alternative 3D based scheme to first identify nuclei, employ the viscous watershed method to determine cell extents and then track the cells and identify mitotic events [18].

Mosaliganti et al. have built the ACME system to detect the planar structure of the cell membrane, and employed it to segment the spatial extent of cells in zebrafish embryos [19]. The EDGE system has been developed to segment and track epithelial cells in Drosophila embryos using a spatial alignment of 2D slices of a membrane marker to reconstruct the structure of the plasma membrane. The 2D slices are first denoised and thresholded, and then polygonal approximations are fitted to the segmented regions and stacked in 3D. Finally, the cells are tracked based on a custom polygon matching method [8]. A new version of the software called EDGE4D was recently published [20], employing fluores-cently labeled nuclei in the cell segmentation.

However, our final goal is to study the role of molecular forces remodeling AIs so we sought to develop methods to segment and track cells in epithelial tissues using only A] markers. Although the AIs arrange in a belt-like structure, methods to detect this structure in 3D have not been reported. To the best of our knowledge, segmentation and tracking of AIs have only been done using 2D max-intensity projections of confocal time lapse data. Watershed algorithms are employed by Packing Analyzer [6] and by SeedWater Segmenter [21] to segment the AIs in 2D. A 2D simplification of the 3D Partial Differential Equation framework introduced in [15] has been reported in [22]. AIs assume the role of membranes when computing cell extent properties from cell nuclei. All the above methods are not easily generalized to 3D because the AIs do not surround the cells in 3D as they do in 2D.

In morphogenetic epithelia, E-cad molecules are constitutively endocytosed and recycled back to the cell surface providing epithelial cells with the plasticity to dynamically rearrange cell shape and cell-cell contacts [23—25]. As a result, E-cad is highly enriched in endocytic vesicles that can produce false positive membrane detections. At the cell surface, the distribution of E-cad is not always uniform increasing the complexity of detecting the AIs as a continuous structure. Certain cell types such as sensory cells accumulate high levels of E-cad thus increasing the complexity of segmenting these cells. Membrane detection methods have to distinguish between membrane and cytoplasmic signal, bridge discontinuities in signal distribution, and tolerate high signal intensity in certain locations.

Cell tracking methods in epithelial tissues need to identify cells as they change shape and contacts with their neighbors. To ensure cell track consistency it is also important to detect dividing cells and properly label the two daughter cells of mitotic events instead of creating a new track for one of them. Similarly, when a cell dies, the event should be properly identified and distinguished from a cell that leaves the scene. Tissue drift causes new cells to enter the scene and others to leave it. The cells that enter and leave the scene should be distinguished from cells that are generated by mitosis and lost by apoptosis, respectively. Tissue drift also introduces a global motion that should be compensated for.

We chose these tissues for their unique morphological and developmental characteristics. The leg is a tube-like epithelium, which narrows and so elongates dramatically during development, while the epithelium of presumptive joints invaginates by apical constriction [26]. By contrast, the notum is a more planar epithelium that undergoes more subtle morphological changes [27]. The proposed approach to detect the AIs is independent of a prior detection of nuclei. Instead, cells are detected and tracked based on segmentation of the AIs and the determination of the connectivity between cells in the tissue. The methods we developed: 1) directly detect the vertices of the AIs where three or more cells meet in image volumes; 2) compute a planar graph approximation to the AIs network structure; 3) compute cell locations in the tissue and spatial associations between cells; 4) track cells and discover mitotic and apoptotic events. Additionally, we show how the computed motion of cell centroids can be used to describe the global dynamic behavior of a tissue. Finally, we provide a framework to assess the performance of the methods over a range of parameter values.

Results

Segmentation and Tracking of Adherens Junctions in 3D

We have developed a system to segment the AIs of epithelial cells highlighted by A] markers such as E-cadzzGFP and find the correspondence among them in adjacent temporal frames. The system is able to segment and track from just a feW to hundreds of cells in 3D under different imaging conditions, including planar epithelial tissues such the Drosophila notum or tube-like epithelial tissues such the Drosophila leg.

First, confocal time lapses of epithelial tissue morphogenesis (Fig 1B) are acquired. To enhance the representation of the Als, timelapses are deconvoluted to eliminate both out-of-focus signal, and the Poisson noise introduced by the photodetector during imaging (Fig 1C). If required, the images are locally equalized to obtain an even intensity of the E-cad::GFP signal employing Contrast Limited Adaptive Histogram Equalization (CLAHE) [28].

A pair of filters based on the second order spatial derivatives of the image intensity are employed to measure the likelihood of each voxel being part either of a verteX or an edge [29] (Fig 1D). The vertices of the graph are initialized With the local maxima of the output of the verteX detection filter. A level-set based region growing algorithm is employed to establish the connectivity among vertices in the input volume, adding an edge to the A] graph for each pair of adjacent vertices (Fig 1E). See materials and methods for details.

The A] graph is planar and has a face for each cell. The Cell graph (Fig 1F) is obtained as the dual of the A] graph, With a verteX for each cell and an edge for each pair of adjacent cells. Each cell is associated With a list of vertices in the A] graph defining its spatial extent. The spatial moments of this polygonal representation provide a set of spatial features for each cell such as perimeter, area, centroid, size or rotation.

To this end we employ a min-cost maX-flow cell tracking framework capable of detecting cellular mitosis and apoptosis events [10] (Fig 2). We consider that the detection of these events is as important as the correct association of cells among frames, although their relatively low frequency makes their detection even harder. However, understanding where and when cells proliferate and die is as important as studying how cells change their shape. Our system is able to successfully detect cellular apoptosis as the area of the dying cells tends to vanish and to detect cellular mitosis as the spatial moments of the parent cell are similar to the moment of the union of the daughter cells.

We developed tools to fix errors in the output of the verteX location, edge segmentation and cell tracking algorithms (Screenshot shown in 811 Fig). Error correction is done after each step as errors at one step are highly magnified in the next. Error correction ensures that the data to be analyzed accurately reflects the observed data.

Case Study: Studying Morphogenesis of the Drosophila Notum

The timelapse captured the dynamics of cells in an area around the midline of the mid-scutum. We have employed the system to recover A] graphs and cell graphs (Fig 3B and 3C), identifying the cells in the tissue and establishing the temporal correspondence among them (Fig 3D). The output of the system included some error that we corrected manually with validation tools that we developed to obtain accurate data. A visual inspection of the recovered cell trajectories shown in Fig 3D reveals a velocity gradient, increasing from posterior (left) to anterior (right). To understand the reason for this, we have built a model to quantify the process that contributes to the global deformation of the tissue at this developmental stage.

Because the epithelium of the notum is relatively planar, we dropped the third dimension of the cell trajectories and we built a strain rate tensor for the tissue, taking the velocity vector field of the cell centroids. The strain rate tensor compo-3F. Their values confirm that the velocity vectors depend on the position of each cell in the tissue and their signs show they are higher the farther they are from the posterior scutellum. We decomposed the strain rate tensor at each frame into their symmetric and antisymmetric components to respectively compute the expansion coefficient 8 (Fig 3G) and the coefficient of rotation of the tissue 9 (Fig 3H). The expansion coefficient is always positive as the surface area of the tissue increases during the process. The coefficient of rotation remains almost constant, reporting a global rotation of the tissue of about one degree between frames,very likely resulting from a drift of the tissue in the medium. Prior to computing the strain rate tensors we have preprocessed the cell centroid trajectories to remove the noise in them employing a Kalman filter [32] (see Materials and Methods). The Kalman filter assumes that cells obey a linear motion model among frames, providing a smoothed estimation of the position and velocity of each cell at each instant.

Performance Evaluation

Performance metrics were computed by comparing the outputs of the system to a ground truth, which was established through manual correction of the system output. For the notum dataset we corrected (added, moved or deleted) about 100 vertices for a total of 824 vertices and 60 edges (added or deleted) for a total 1207 it contains. For the leg dataset we corrected 120 vertices for a total of 403 and 93 edges for a total of 586. The aim of this performance evaluation framework is not only to assess the quality of the proposed algorithms but also to provide a comparison framework to assess the quality of future alternatives.

We limited the assessment to the first frame of the notum and leg time lapses. A verteX in the ground truth is considered as detected if the algorithm marks a vertex location closer than x/i voxels from the manually annotated position, i.e., a vertex is detected if there is a detection in any of the voxels around it. The verteX location method depends on the interval of spatial scales employed to control the size of the A] vertices detected and a detection threshold controlling the strength of the detected vertices. The spatial scale interval is easy to set up visually to an acceptable value, so we only assess the performance regarding variations in the verteX detection threshold. See supplementary S4 Fig for guidelines on how to select the proper spatial scales.

Briefly, Precision measures the number of truly detected vertices found With a given threshold level, While Recall measures the number of vertices undetected for that threshold level. The curve for the Notum dataset has an anomalous behavior for low recall values that arises from the very high values of the Vertexness function at sensory bristle cell locations that are mistaken for A] locations. S6A Fig presents an example of this effect. Another common verteX location error that is produced at AIs are indentations between adjacent vertices as shown in 86B Fig. These are likely to arise by displacement of cell contacts by contractile forces generated by the actomyosin cytoskeleton. Vertex locations errors are more common in areas where cells are more parallel to the Z axis, where they are more difficult to locate due to voxel anisotropy and in areas with a low signal to noise ratio (Fig 4C and 4D).

The A] graph was initialized with the ground truth vertex locations curated for each frame in the time lapses. Fig 4E presents the Precision-Recall curves generated for the different data sets for different propagation thresholds T5 of the algorithm employed to detect edges. Fig 4F shows the variations of the performance according to the value of the threshold. The method achieves high detection scores detecting vertices in the different datasets as soon as the threshold is given an appropriate value.When the threshold is too low many cells are not be detected, while when the threshold is too high extra cells that don’t exist are detected. The optimal detection results are shown in Fig 4G and 4H. Note that areas surrounding bristles and areas with low signal-to-noise ratio produce more detection errors (Fig 4G and 4H). In addition, as the orientation of edges become more parallel to the Z-axis they become more difficult to detect.

To this end we employed all the timelapses. The method depends on 11 parameters to compute the association costs (see materials and methods for detail). As it is not practical to adjust all values for each new dataset, we searched for a combination of parameter values that works wells for the different tissues. We have exhaustively explored the parameter space and have found a combination that provides accurate tracking results across the different scenarios, with a mean AF 1 = 0.93. This represents a very high value as the AH metric highly penalizes failures in rare events such as mitosis and apoptosis. Fig 41, 4], 4K and 4L present cell trajectories recovered from the different timelapses. The position over time of each cell is projected in 2D. The variations of the performance according to the perturbation of the different parameters are provided as supplementary material (S7 and 88 Figs).

The performance drops to a mean AF 1 — score = 0.822569. Computing the AF 1 for the complete and decimated data set with the new parameters reveal that the AF Icomplete = 0.92 and the AF 1 decimated = 0.736818. This shows that the sampling rate of the data set highly influences tracking performance. Interestingly, if we compare the values found for the weights in previous and current experiments (S9 Fig), the weight given to the distance among centroids is reduced and the weight given to other features such as cell perimeter and rotation is increased. Additionally, we compared a 2D simplification of our system to the SeedWaterSegmenter, which utilizes a watershed algorithm to detect cell extents [21] , in the detection of cells, and found that our method performs better in obtaining cell locations. See supplementary material for further experimental details.

Discussion

We have shown that the pipeline detects and tracks cells with a high accuracy when properly parameterized. Clearly the most difficult part of the segmentation process is to locate the A] vertices, as indicated by a lower Fl-score, compared to the detection of edges and the computation of correspondences across frames. However, we described the errors produced by the verteX detection algorithm, so it should be straightforward to manually correct these errors to improve the accuracy of segmentation. We have also provided an example of how to use the generated cell centroid trajectories to study global cellular behaviors during tissue morphogenesis.

Thus, our methods should be generally applicable to the analysis of epithelial morphogenesis in developing embryos in a range of species. We are currently using the system to explore the morphogenesis of the epithelium of the notum and leg imaginal discs to detect patterns of cellular behavior that correlate with the regional subdivision of these structures. We are particularly interested in exploring the behaviors that contribute to the narrowing and elongation of the epithelium of the leg imaginal disc and the relative contribution of presumptive joints and segments to this process [33, 34]. In the notum, we are particularly interested in the relationship between the anterior-posterior subdivision of the notum and the pattern of cell rearrangements in each region [35, 36]. Addressing these questions will help us understand how patterning of epithelial sheets at early stages of development affect epithelial dynamics and mechanics that generate the final form of adult structures at later stages.

The detection of cell outlines and motion can be used to examine the evolution of cell topologies and the relative contribution of subcellular and cellular deformations to tissue morphogenesis. As the system stores the locations of vertices, edges and centroids in the data structure, it could be used to parameterize force models of epithelial remodeling such as cell vertex models [37, 38]. The detection of mitosis and apoptosis can be used to investigate the contribution of lineage relationships and patterns of apoptosis to tissue development.

The activity of contractile actomyosin networks and their coupling to the AIs can alter cell shape and cell contacts, as well as spatial patterns of cell proliferation and apoptosis. In contrast, adhesive forces mediated by cell adhesion proteins promote formation and expansion of cell contacts. Differential regulation of these contractile and adhesive networks in space and time can be used to deform epithelial sheets in predictable ways. Understanding how the activities of these networks are related to the changes in cell behavior can suggest mechanisms that drive tissue morphogenesis. The proposed system might be extended to relate changes in abundance, dynamics and polarity of contractile and adhesive molecular assemblies to the changes in cell and tissue shape. The segmentation and tracking of AIs in 3D and the curation of vertices, edges, cells and trajectories would enable the extension of the pipeline to relate molecular dynamics with cell and tissue remodeling for multiscale analysis of epithelial tissue morphogenesis.

Materials and methods

Tissue Preparation and Image Acquisition

For notum imaging, the pupal case was removed to expose the head and dorsal thoraX at the desired developmental stages. For mounting, a slab of 1.5% agarose gel was placed atop a 30 mm coverslip. 2 intact pupae were mounted in tandem on the coverslip through a slit made in the agarose slab. A silicone gasket (Sylgard 184, Dow Corning) was fitted to surround the agarose slab and a chamber constructed from acrylic was fitted atop the gasket to seal the The precise settings for acquisition of time lapse movies varied for each eXperiment but were dictated by the need to minimize photobleaching of the E-cadzzGFP reporter, tissue damage and image corruption by noise, and to maximize image resolution (supplementary Sl Table). Approximately 10 sections at 1 to 1.5 micrometer intervals with 50% overlap were collected every 10 minutes to provide sufficient temporal resolution to track the cells and capture mitosis and apoptosis events.

We employ an Elastic-Net Prior [42] to perform the deconvolution, as many of the voxels are eXpected to be zero. CLAHE was employed to locally normalize the intensity of the 3D volumes, enhancing the A] structures.

Adherens Junction Segmentation

It is an undirected graph defined by the tuple G A 2 (VA, EA), Where VA is the set of AIs vertices—the points Where three or more cells touch—and the set of edges EA is such that there is an edge e 6 BA for each pair of contiguous vertices in the tissue. We build the A] graph from the output of the pla-teness function 73(x) proposed by Mosaliganti et al. [29] measuring how likely each voxel is part of the AIs and a vertexness function V(x) measuring how likely each voxel is an A] vertex. Both functions are built from spatial second order derivatives of the image intensity. See supplementary material for an accurate description of the functions.

A super-vertex is the 3D region around each verteX v 6 VA such that v is the vertex minimizing a time of travel cost function to be introduced later. We build the formal definition of a superverteX from the notion of Voronoi Region around a vertex. The Voronoi Region [43] around a vertex v 6 VA is the subset of all the points in the 3D space that are closer to 1/ than to any other vertex in VA:

For each point on the Voronoi Region of a vertex v we compute the time needed to travel to 1/, Where the travel speed at each point x E R3 is given by F(x) 2 73(x) + V(x) > O. The time needed to travel from each point x E Voronoi(v) to the verteX v is found as the solution to the partial differential equation: with boundary condition T (v) = O. This equation is a well known partial differential equation known as the Eikonal equation. It is employed to obtain the travel time for a wave to reach a point when traveling through a scalar field. Based on the notion of a Voronoi region around a vertex and the Eikonal equation employed to measure the travel-time cost from a point to its nearest verteX we provide a formal definition for a superverteX. The superverteX of a verteX v 6 VA is defined to be the subset of the points in the Voronoi Region of v such that the solution T: of Eq 2 with initial verteX v is lower than a threshold T5:

An edge (v, w) is added to EA if and only if superverteX(v) and superverteX(w) touch. The computed A] graph is shown in Fig 5C.

Cell Identification

However, the A] graph does not explicitly model cells and their spatial extent and connectivity. To build a tissue model at the level of cells and edges and perform inferences, it is possible to eXploit the properties of the A] graph, and transform it into another graph to represent cells and their neighborhood relationships.

Intuitively, a graph is planar when it might be drawn in a plane with none of its edges crossing one another. Formally, a graph G = (V, E) is planar if there exist an embedding of the vertices in R2, f : V —> R2 such that for all pairs of edges (61, b), (c, d) E E, with a 7E 19 7S 6 7E d E V, the line segment from f (a) to f(b) does not cross the line segment from f(c) to f(d). The edges of any planar graph G divide the space into regions called faces. The main implication for a graph being planar is the eXistence of an equivalent dual graph. For every planar graph G drawn in a plane there exists a graph G* whose vertices correspond to the faces of G and there is an edge between every pair of adjacent faces, i.e. the faces have at least one edge in common.

The Cell graph GC 2 (VC, EC) is defined as the dual graph of G A, removing the verteX corresponding to the outer face. Each one of the vertices v 6 VC represents a cell. Each one of the edges (v, w) 6 EC represents cell adjacency among 1/ and w. Thus, for each cell in the tissue there is a vertex in the Cell graph. A cell is defined by the pair (c,A) Where c 6 VC is the verteX representing the cell in the cell graph and the set A = a1, . . ., a K, ak 6 VA is the clockwise ordered list of the vertices of the A] graph delimiting the cell. The polygonal representation of Aj allows the derivation of moments such as cell perimeter p(c) 6 1K”, cell area a(c) E 11?, cell centroid b(c) E R3, cell Width w (c) 6 1K“, length h(c) 6 RV or rotation r(c) E (—71, 71'). All these features are later employed to match cells among adjacent frames for cell tracking. This representation also provides an eXplicit model of cell neighborhood necessary to compute tissue measurements such as those proposed in [46] or [47].

Cell Tracking

The developed cell tracking methods rely on the methods we have developed to identify cell locations and their spatial properties. The method we employed to solve the cell correspondence problem among frames is a variation of the coupled min cost-max flow framework reported in [10]. The solution to the cell correspondence problem is obtained as the solution to a flow transportation problem in a directed graph. A total of N + M units of flow need to be sent from a source vertex T to a sink vertex T— traversing a network formed by set of vertices and a set of arcs connecting the vertices that encodes the cell association problem. Arcs have a maximum capacity and an associated cost for sending units of flow through them. The set of arcs minimizing the cost for sending the N + M units of flow through the network gives the solution to the cell correspondence problem. Flow has to be preserved among the arcs of the network, i.e., the same amount of flow that gets into a vertex needs to be sent to others, except at source and sink vertices. The cell association cost attached to the arcs is computed from the cell moments obtained in the previous section. The cost of associating a cell to a given cell in the next frame is computed as a weighted sum of the difference among their moments. In case of a mitosis event, the cost is computed in a similar way but from the union of the polygons corresponding to a pair of adjacent hypothesized sibling cells. The weight for cells entering and leaving the scene is proportional to the distance from the centroid to the tissue perimeter. The weight for the apoptosis hypothesis is proportional to the area of the dying cell.

The proposed framework differs from the original in the way correspondence hypotheses are formulated. As we track cells embedded in an epithelial sheet rather than freely moving particles we exploit the neighborhood relationships among cells to drop association hypothesis not corresponding to adjacent cells.

Analysis of Notum Morphogenesis

The cell centroid trajectories extracted by our cell tracking system are not smooth, as it solves temporal correspondences among the observed positions of the cells Without performing any correction to the observed positions. We employ a Kalman filter to model the joint probability distribution P(O, H) of the sequence of Where the state transition distribution P(ht+1|ht) is a Gaussian distribution such hm ~ MAht, 2 H) With covariance 2 H and the observation distribution P(ot|ht) is also a Gaussian distribution such of E MBht, 20) With covariance matrix 20. The state transition matrix A and the observation matrix B are respectively defined as:

In particular, Viterbi algorithm returns the sequence H3k maximizing the posterior probability distribution P(H | 0). Using this method we obtain a smoothed estimation of the position and velocities of the cells that are employed in our next analysis.

We build a strain rate tensor [31] to study the relationship among the position p,- 2 (xi, yi)T and the velocity 1. = 1., i)T of a cell i. The strain rate tensor is a Iacobian matrix I satisfying:

The strain rate tensor might be decomposed into a symmetric part 8 = % ( I —|— I T) and an antisymmetric part Q = i ( I — I T). The symmetric part provides the eXpansion coefficient of the vector field as trace(€), While the antisymmetric part provides the rotation coefficient as acos(Q 12). The temporal evolution of these values provides a quantitative description of the global behavior of the tissue.

Performance Evaluation

There are two type of errors that detection algorithms might produce, known as False Positives (FF) and False Negatives (FN). A FF is produced when the algorithm marks a detection that is not real, while a FN is produced when the algorithm misses a real detection. Thus, to assess the quality of a given set of detections, it is possible to compute three metrics known as Precision, Recall and F1 measure defined as follows:

Detectors commonly have a detection threshold (the proposed here do) that has to be adjusted a priori to some value and conditions Precision and Recall values of the algorithm. Note that Precision and Recall are conflicting measures: a high detection threshold produces high precision and low recall (very few but true detections, but many missed detections), while a low detection threshold produces low precision but a high recall (not a lot of missed detections, but a lot of false detections). Thus, algorithms assessment should include tests for many threshold levels to produce Precision-Recall curves. The sensitivity of parameters was eXplored based on the changes of the F1 measures.

To this end we consider it as a label prediction problem. For each cell in a given frame the task is to predict the correspondence of a cell with the same cell in next frame, if it should disappear following apoptosis, or if it should be associated with two new cells following mitosis. Thus, an average F1 measure might be obtained for each tracking configuration as the harmonic mean of the individual F1 measures:

Supporting Information

The output of the plateness function proposed by Mosaliganti et al. is employed to detect the Als [29]. AC) Slices of the outupt of the plateness function 730(X) computed for different values of a (A a = 0.14, B a = 0.45 and C a = 0.60). The value of 73‘7 (X) depends on the signal intensity variations among Als produced by the imaging process and variations in the structure of AIs, achieving the maxima at real scale a. (TIFF)

A) A slice of the planarity response at a = 0.14 before applying diffusion. Planarity values are noisy and the response at the A13 is uneven. B) The same slice after 50 iterations of membrane enhancement diffusion. Most of the noisy values have been removed, the detected AIS obtain a more uniform intensity profile, and some of the cytoplasmic signal has been cleared. (TIFF)

[29] decreases around AIs vertices. Arrows point to voxels where the effect is well observed. This feature is exploited to detect cell vertices. The image shows a slice of the ouput of the plateness function 730(x) with a = 0.6.

The output of the Vertexness function here proposed has been overimposed in black over the plateness function outputs shown in 81 Fig at the corresponding scales (A a = 0.14, B a = 0.45 and C a = 0.60). Note that at higher scales (B and C) vertices Which are close to each other tend to merge, While at lower scales vertices tend to appear at non-vertex locations along the AIs. A set up such as the one proposed in panel B is desired as it provides an accurate detection of AIs and vertices. In A the scale is too lOW resulting in high noise, While in C the scale is too high resulting in detection of blurred features. (TIFF)

A) TWO cells, 12 divides to produces r2 and r3. B) The graph we built to represent all the correspondence hypotheses. Arcs in red represent cell association, in blue cells entering the scene, in green mitosis, in pink apoptosis and in gray cells leaving the scene. C) The arcs of the graph eXpected to represent the desired solution (TIFF)

Details from Fig 4C. Green vertices represent true detections, blue, missed detections, and red, false detections. A) Common pattern of vertices detected at bristle locations, Where many vertices are not detected but one is falsely detected at the center. B) Indentations appear along edges between vertices as regions With high curvature that are detected as vertices. (TIFF)

Global shows the harmonic mean of the Average Fl-scores obtained for the different datasets. The difference at the optimal between the global measure and the Average Fl-scores of each dataset is not significant, but the global measure drops fast as parameter values deviate from the optimal. A) Centroids. B) Area. C) Perimeter. D) Width. E) Rotation. F) Length.

Similar to S7 Fig, global shows the harmonic mean of the Average Fl-scores obtained for the different datasets. The difference at the optimal between the global measure and the Average Fl-scores of each dataset is not significant, but the global measure drops fast as parameter values deviate from the optimal. A) Cell Association. B) Cell entering the scene. C) Cell mitosis. D) Cell Apoptosis. E) Cell leaving the scene.

wb, wa, wp, ww, we and w, are respectively the weights given to the the distances among cell centroids, areas, perimeters, Widths, rotations and heights to compute cell association costs. The distance between cell centroids (wb) receives less importance in tracking cells in decimated data, While the distance between cells perimeters (WP) and rotation (we) receives more importance.

TTT obtains a F1 score higher than SeedWaterSegmenter. These F1 scores have been obtained after tuning the parameters of SeedWaterSegmenter and TTT to obtain the highest score for each system. For SeedWaterSegmenter only one parameter was tuned, compared to seven parameters that were tuned for TTT. Although SeedWaterSegmenter is less accurate than TTT, it is easier to use. (TIFF)

Vertices and edges of the A] graph can be edited to correct for segmentation errors and obtain accurate data for further analysis. Vertices and edges are selected employing the mouse. A 3D cursor is employed to remove, add or move vertices. Note the red vertex selected by the 3D cursor.

Additional material describing technical details of the proposed system. (PDF)

Description of the datasets employed to evaluate the segmentation and tracking pipeline. (TIFF)

Acknowledgments

We would like to thank Adam Martin (MIT, Biology), Eric Miller and Brian Tracey (Tufts, Medford) for discussions and advice.

Author Contributions

Performed the experiments: RC VM VH SDS BHdM. Analyzed the data: RC. Contributed reagents/materials/analysis tools: BHdM SDS VH. Wrote the paper: RC ID VH SDS.