Robust cell tracking in epithelial tissues through identification of maximum common subgraphs

Tracking of cells in live-imaging microscopy videos of epithelial sheets is a powerful tool for investigating fundamental processes in embryonic development. Characterizing cell growth, proliferation, intercalation and apoptosis in epithelia helps us to understand how morphogenetic processes such as tissue invagination and extension are locally regulated and controlled. Accurate cell tracking requires correctly resolving cells entering or leaving the field of view between frames, cell neighbour exchanges, cell removals and cell divisions. However, current tracking methods for epithelial sheets are not robust to large morphogenetic deformations and require significant manual interventions. Here, we present a novel algorithm for epithelial cell tracking, exploiting the graph-theoretic concept of a ‘maximum common subgraph’ to track cells between frames of a video. Our algorithm does not require the adjustment of tissue-specific parameters, and scales in sub-quadratic time with tissue size. It does not rely on precise positional information, permitting large cell movements between frames and enabling tracking in datasets acquired at low temporal resolution due to experimental constraints such as phototoxicity. To demonstrate the method, we perform tracking on the Drosophila embryonic epidermis and compare cell–cell rearrangements to previous studies in other tissues. Our implementation is open source and generally applicable to epithelial tissues.


Introduction
Live-imaging microscopy is a powerful, and increasingly quantitative, tool for gaining insight into fundamental processes during embryonic development [1][2][3]. Quantitative information on cell growth, proliferation, death, shape changes and movement extracted from live-imaging reveals how such processes are regulated to give correct tissue-level behaviour. This approach has been particularly successful in characterizing the growth and patterning of embryonic epithelial tissues in a number of model organisms [4][5][6][7][8][9].
A common experimental technique for visualizing cell shapes in an epithelial sheet is to fluorescently tag a molecule marking cell boundaries, such as E-cadherin (figure 1a). The analysis of time-lapse microscopy data obtained from such tissues is extremely challenging [2,3], especially in cases of imaging data of rapidly evolving tissues, and when limitations of, for example, microscope speed, imaging resolution or phototoxicity prohibit the creation of datasets with high temporal and spatial resolution.
The analysis of time-lapse microscopy data comprises two major steps: segmentation and tracking (registration). Segmentation must be performed for each frame of a video and involves the identification of objects and landmarks, such as cell shapes (figure 1b). Automated segmentation is hindered by various factors such as noise in fluorescent signals, uneven illumination of the sample or overlapping cells in a two-dimensional projection. Often, manual correction is necessary to address over-segmentation, where too many cells are detected, or under-segmentation, where too few cells are detected [10 -12]. Tracking involves the association of segmented cells across video frames (figure 1c) and requires resolving cellular movement, cell division, cell death and cells entering and leaving the field of view [12].
Numerous algorithms are available for the segmentation and tracking of cellular-resolution microscopy data [10,11,13]. Common methods for cell tracking use optimization techniques to minimize differences in cellular properties between two frames [11,[14][15][16][17]. The min-cost max-flow algorithm [14] uses linear integer programming to minimize differences in cell areas, perimeters, orientations and locations between frames, whereas multiple-parameter tracking [15] employs global optimization to minimize differences in cell shapes as well as locations. By contrast, multi-temporal association tracking [16,17] minimizes differences in cell locations and sizes by using a probabilistic approach that finds the most probable extension to existing cell trajectories. Chain-graph models [18] minimize differences in cell velocity while overcoming mis-segmentation by verifying that each segmented object continues or begins a cell trajectory in successive frames. Optical flow (warping) between successive frames can be used to guide cell tracking as well as segmentation [19]. It is also possible to combine segmentation and tracking of two-dimensional microscopy videos by interpreting time as a third spatial dimension and employing three-dimensional segmentation techniques [20]. The nearest-neighbour method associates two cells in consecutive frames with each other if their respective centroids have minimal distance within the field of view [10], or if their overlap in pixels within the field of view is maximal [21,22]. Particle image velocimetry, a technique originally developed to analyse fluid flow [23], has also been employed to track cells in epithelial tissues [24].
Software implementations and computational tools for cell tracking include FARSIGHT [25] (segmentation only), SeedWaterSegmenter [10] (nearest-neighbour tracking), ilastik [18] (chain-graph models), Tufts Tissue Tracker [11] (min-cost max-flow algorithm), Tracking with Gaussian Mixture Models [26] (nearest-neighbour tracking), Packing Analyzer [27] ( particle image velocimetry) and EpiTools [13] (nearest-neighbour tracking). These algorithms and software tools primarily rely on there being small differences in cell positions and shapes across consecutive images. Their performance is therefore hindered when analysing data from in vivo studies where phototoxicity provides a barrier to high-temporal resolution imaging [28 -30]. To address this limitation, we propose a novel algorithm for cell tracking that uses only the connectivity of cell apical surfaces (figure 1). By representing the cell sheet as a physical network in which each pair of adjacent cells shares an edge, we show that cells can be tracked between successive frames by finding the maximum common subgraph (MCS) of the two networks: the largest network of connected cells that is contained in these two consecutive frames. It is then possible to track any remaining cells based on their adjacency to cells tracked using the MCS. Our algorithm does not require the tuning of parameters to a specific application, and scales in sub-quadratic time with the number of cells in the sheet, making it amenable to the analysis of large tissues.
We demonstrate here that our algorithm resolves tissue movements, cell neighbour exchanges, cell division and cell removal (for example, by delamination, extrusion or death) in a large number of in silico datasets, and successfully tracks cells across sample segmented frames from in vivo microscopy data of a stage-11 Drosophila embryo. We further show how our algorithm may be used to gain insight into tissue homeostasis by measuring, for example, the rate of cell rearrangement in the tissue. In particular, we find a large amount of cell rearrangement within the observed dataset despite the absence of gross morphogenetic movement. The remainder of the paper is structured as follows. In §2, we describe the algorithm for cell tracking. In §3, we analyse the performance of the algorithm on in silico and in vivo datasets. Finally, in §4, we discuss future extensions and potential applications.

Material and methods
In this section, we provide a conceptual overview of the core principles underlying our cell tracking algorithm. We focus on providing an accessible, non-technical description rather than including all details required to implement the algorithm from scratch. A comprehensive mathematical description of the algorithm is provided in the electronic supplementary material.
The input to the algorithm is a set of segmented images obtained from a live-imaging microscopy dataset of the apical surface of an epithelial cell sheet. For each image, the segmentation is assumed to have correctly identified which cells are adjacent and the locations of junctions where three or more cells meet. Various publicly available segmentation tools can be used for this segmentation step, for example, SeedWaterSegmenter [10] or ilastik [18]. The segmentation is used to generate a polygonal approximation to the cell tessellation (figure 1b,c). Such approximations are an adequate assumption for many epithelia [11,31 -34].
Our algorithm tracks cells by interpreting the polygonal representations arising from the segmentation as networks (graphs) of cells. Examples of such networks are shown in figure 2a. In this representation, each cell corresponds to a vertex of the network, and two vertices are connected by an edge if the corresponding cells are adjacent. Our algorithm tracks cells across consecutive images by aligning the networks of cells that correspond to these images. This network alignment is achieved in three steps. First, we generate an initial tracking for subsets of the cells in each pair of consecutive images by finding the MCS between the two corresponding networks (figure 2b). Second, this MCS is reduced to avoid tracking errors (figure 2c). Third, remaining untracked cells are tracked based on their adjacency to cells within the MCS (figure 2d ). In the final output of the algorithm, each tracked cell of a frame is paired with exactly one cell in the subsequent frame.
The key step in this network alignment approach is the identification of an MCS [35,36]. An MCS comprises the largest sub-network that is contained in two larger networks; thus finding an MCS can be understood as recognizing patterns of connections that are preserved between two networks. In this work, the structure of the MCS roughly corresponds to cells that do not rearrange between consecutive images, except for a few cells at its boundaries.
In figure 2b, we visualize the MCS generated by our algorithm as a collection of green (light) and purple (dark) cells. Most of the highlighted cells in figure 2b are tracked correctly by the MCS. Three cells in each frame are marked by a yellow (bright) dot. Within the two cell networks, these cells are members of the MCS. However, these cells are not tracked correctly by the MCS. This mismatch arises as the MCS is found based on the connectivity of cells within the network alone. The fewer connections a cell has to other cells in the MCS, the less information about the cell's position and shape is encoded by these network connections, and so the greater the possibility of mismatches. To avoid such tracking errors, we remove any cells that have only a few connections within the MCS, as well as small isolated clusters of cells. All cells that are removed from the tracking in the second step of our algorithm are shown in purple (dark) in figure 2b. In figure 2c, we highlight the cells that are tracked after applying this second step of our algorithm.
Cells that are untracked after reducing the MCS are then tracked based on their connections to previously tracked cells. This last step of our algorithm comprises starting from the MCS and iteratively 'growing' the set of tracked cells by adding cell-to-cell matches to the tracking that maximize the number of preserved connections to other tracked cells. In this step, the algorithm also resolves cell neighbour exchanges, cell removals and cell divisions by identifying changes in the network structure that are characteristic of these events. For example, a cell neighbour exchange corresponds to the deletion of a network connection while a new connection is added.

The maximum common subgraph is identified through repeated seeding and iterative extension
In computer science, MCS finding has been known to be an NP-hard problem [35,36]: the time to find an exact MCS of two networks increases exponentially with the size of the networks, which poses a computational barrier to the use of MCS-finding algorithms in applications. We overcome this computational barrier in this work by constructing the MCS iteratively from the MCSs of smaller subgraphs, exploiting the planar structure of our cell networks to reduce the complexity of the problem.
To start the construction of the MCS, the algorithm identifies a match between two cells in the consecutive images for which the structure of the network of their surrounding cells is identical. Here, the network of surrounding cells is restricted to the network formed by a cell's neighbours and its second nearest neighbours (figure 3b). If no such initial match can be found, the algorithm instead searches for an initial match where only the first-order neighbourhood is preserved, under the condition that this neighbourhood does not touch the boundary of the tissue. This latter condition avoids tracking errors that can occur on the tissue boundary, where cells have few neighbours.
Once the initial match (a 'seed') is found, the algorithm iteratively adds further cells to the MCS. At each step of this iteration, a cell in the first network is picked that is adjacent to the existing MCS and has a minimal number of potential matches. This number of potential matches is determined based on how many cells in the second network have the same number of neighbours as the considered cell while preserving connections to already tracked cells. Among the choice of potential matches, the algorithm identifies an optimal match based on the local network structure of these cells' neighbours. A cell in the second network is identified as an optimal match if the network structure of its neighbourhood is most similar to the cell in the first match. This choice is made based on local MCSs between the neighbourhood-networks of the cell in the first image and each potential match. Note that the optimal choice may exclude the cell from the tracking entirely. In this case, most neighbours are included in the local MCS when the considered cell is not tracked, indicating, for example, a cell removal event. In this case, the cell is not mapped and the algorithm proceeds by inspecting another cell in the first match. Cells in the first frame for which no match in the second frame has been found may be re-inspected at later stages of the algorithm as the size of the identified MCS increases. Once no more adjacent cells can be added to the MCS through this iterative extension, the iteration continues the search among untracked cells in the first network that are not adjacent to the existing MCS. As soon as at least one cell has been added to the MCS in this way, the algorithm again restricts its search to adjacent cells. The algorithm halts once no further cell-to-cell matches can be found. During the construction of the MCS, the algorithm ignores any potential cell-to-cell matches where the corresponding cell centroids are more than a cutoff distance d max apart within the field of view. Throughout the paper, we choose d max to be 10 average cell lengths.
Once the MCS is complete, any cells that have less than three isolated connections to other cells in the MCS are removed from the tracking. Any clusters of 10 or fewer cells are also removed from the tracking result. Both of these steps help to minimize tracking errors (figure 2b,c).

Cells are added to the tracking result by inspecting connections to previously tracked neighbours
Through the identification of the MCS, the algorithm tracks most of the cells that do not rearrange between consecutive frames. When adding cells to the tracking, the algorithm ensures that a cell cannot gain more tracked neighbours between consecutive frames than the number of tracked neighbours preserved between these frames. The algorithm also requires a cell to have at least two tracked neighbours in order to be added to the tracking in this way. Once all possible cells have been tracked, the algorithm resolves division events. Division locations can be identified as regions in the second frame that contain more cells than the corresponding region in the first frame. As the algorithm will have found exactly one match in the second network for each tracked cell in the first network and vice versa, there are thus untracked cells in the second frame wherever a cell divides between two consecutive frames. The algorithm attempts to resolve division events by identifying changes in cell-to-cell connectivity that are characteristic to dividing cells (figure 4). For example, two cells adjacent to each division must gain a neighbour (grey cells in figure 4a), and in many cases the mother and daughter cells are easily identified as the cells that are shared neighbours of these cells adjacent to the division event. However, one of the daughter cells may be four-or three-sided (figure 4b,c). In these cases, the algorithm is not able to determine the mother and daughter cells based on their network properties alone. Instead, the algorithm takes the geometric shape of the cells into account. The mother and daughter cells are chosen by identifying which pair of potential daughter cells has the closest position to their potential mother cell.
Cell deaths are identified as cells in the first frame that do not have a tracked match in the second frame and that are not on the boundary of the region of tracked cells.

Code availability
The code used in this article is publicly available under the 3-clause BSD licence as the MCSTracker project (https:// github.com/kursawe/MCSTracker). The project is implemented in pure Python, employs unit testing [37] and is fully documented. Graphs in our code are represented using the Python package NetworkX [38].

Generation of in silico datasets
To test the algorithm, we generate in silico datasets that include examples of cell divisions, removals and neighbour exchanges, as well as tissue movement. These datasets are generated using Voronoi tessellations modified using Lloyd's relaxation, which resemble cell packings in a variety of epithelial tissues [33,39].
To generate polygonal patterns of size m Â n, where m and n are natural numbers, (m þ g) Â (n þ g) Voronoi seeds are distributed uniformly at random in a two-dimensional domain V of width m þ g and height n þ g (figure 5a). Here, g denotes the size of a boundary region that is introduced to reduce the impact of the Voronoi boundary on the patterns. The domain V is surrounded by two additional rows of evenly spaced seeds on each side. The inner row is a distance of 0.5 length units to V, and the seed spacing is 1.0. The outer row has a distance of 1.5 to V, and the seeds are shifted parallel to the first row by a distance of 0.5. The Voronoi tessellation of all these seeds is then constructed.
In each Lloyd's relaxation step, the polygons (or infinitely large areas) corresponding to the regularly spaced seeds outside V are removed from the tessellation. Next, the centroid of each remaining polygon is calculated and registered as a new seed. Further seeds are added that again correspond to two rows of evenly spaced seeds outside V. A new Voronoi tessellation is then constructed (figure 5b). This procedure is repeated for L relaxation steps, after which all generated polygons are discarded except those whose centroids lie within a rectangular domain of size n Â m area units whose centroid coincides with that of V (figure 5b).
The polygonal tessellations have approximately m Â n polygons of average area 1.0. During the generation of the tessellations, evenly spaced seeds outside V are added to prevent the occurrence of infinitely large polygons inside V. The boundary of size g is added in between the generated tessellation and the evenly spaced seeds to reduce the effect of the evenly spaced boundary seeds on the tessellation. Throughout this study, we use g ¼ 8 and n L ¼ 4, resulting in cell packings similar to those observed, for example, in the Drosophila wing imaginal disc [33]. We provide further details of how tissue rearrangements are implemented in the Results section.

Experimental methods
Live imaging of cell proliferation was performed in stage-11 Drosophila embryos expressing a tagged version of DE-Cadherin (DE-Cadherin::GFP) using a spinning disc confocal microscope, as described in [40]. For the embryo set-up, a modified version of the standard live-imaging protocol was used [41].

Data segmentation
Microscopy images were segmented using pixel classification in ilastik [18]. The classifier was trained to recognize cell outlines and the segmentation of each frame was manually corrected. A watershed algorithm was used to identify the precise shape of the cell outlines. Each segmented frame was converted to a 16-bit greyscale image where pixels belonging to different cells had different integer values. Polygonal tessellations for the tracking algorithm were generated from the segmented image in two steps. First, all junctions between three or more cells were identified as points where pixels of three or more different cells met; second, vertices were assigned to cells. Then, edges shorter than two pixels (0.5 mm) were removed and replaced by a single vertex at the midpoint of the edge. Finally, polygons at the boundary of the tissue were removed from the simulation. This removal was necessary since cell shapes at the tissue boundary are poorly approximated by polygons due to missing vertices. Note that our algorithm can interpret segmentations saved using either ilastik [18] or SeedWaterSegmenter [10].

In silico testing of the algorithm
To assess the performance of the algorithm, we begin by applying it to in silico datasets that include cell neighbour exchanges, tissue movement, cell removal and cell division. In each case, we compare the outcome of the tracking algorithm to the ground truth.
We begin by assessing the ability of the algorithm to resolve permutations in otherwise identical tissues (figure 6a). In this test, a random tessellation of size nine by nine cells is created as described in the Material and methods section, and integer identifiers c i are assigned to each cell. Next, an identical copy of the tissue is created in which the integer identifiers are randomly shuffled. A ground truth mapping from the first to the second integer identifiers is Neither the MCS detection algorithm nor the post-processing algorithm is able to resolve such mappings, which involve fewer than four cells in each dataset (fewer than five per cent of the tissue). We design four further tests of tissue rearrangements (figure 6b-e). The first test comprises tissue movements between images (figure 6b). In this test, a tissue of size fifteen by eight cells is generated as described in the Material and methods section. Two smaller tissues of width seven units are cut out of this tissue, which each cover the full height of the tissue, and which are horizontally translated relative to each other by a distance of two cell lengths. The position of each three-cell junction in both tissues is shifted such that the x-coordinate of the left-most junction in each tissue is 0.
The second test (figure 6c) generates cell neighbour exchanges, also called T1 transitions [42,43]. In our implementation of T1 transitions, an edge shared by two cells is replaced by a new perpendicular edge (of length l T1 ¼ 0.2 units) such that the local cell connectivity changes (figure 2b). We create two identical copies of a tissue of size nine by nine cells. In the second copy, a T1 transition is performed on an edge in the centre of the tissue.
The third test involves cell removal ( figure 6d ). In this test, we first generate two identical copies of a tissue of size nine by nine cells. In the second copy, we replace the central cell by a vertex shared by its neighbouring cells, a rearrangement similar to the so-called T2 transitions [42]. The final test involves cell divisions (figure 6e). Here, we once again create two identical copies of size nine by nine cells. In the second copy, a cell in the centre of the tissue is bisected by introducing a straight line in a random direction through the centroid of that cell.
For all tests generated in this way, integer cell identifiers in the second tissue are randomly shuffled, and a ground truth is generated. We run 100 realizations of each test case, and compare the tracking outcome to the ground truth. In all cases, we find that cells are tracked correctly, with at most three unmatched cells at the boundary of the sheet.
In figure 6, all cells identified after the cleaning step, in which weakly connected cells are removed from the MCS, are coloured green, whereas cells that are identified by the postprocessing algorithm are coloured red. Note that the exact number of cells that are identified by the post-processing algorithm varies between individual realizations of the tests. In many cases, the cells identified by the post-processing algorithm include cells that are adjacent to those undergoing division, removal or neighbour exchange.
We next analyse the extent to which the success of our tracking algorithm depends on the number of Lloyd's relaxation steps, n L , used to generate the in silico datasets. To investigate this, we iteratively increase n L , thus generating tissues with increasingly homogeneous graph structures, and repeat all tests. We find that the algorithm successfully passes all tests for all values of n L from 4 up to 14.

Algorithm performance for large numbers of cell neighbour exchanges
To assess the performance of the algorithm when applied to tissues exhibiting large numbers of cell neighbour exchanges, we next apply the algorithm to in silico datasets with increasing numbers of cell neighbour exchanges between frames (figure 7). The number of correctly tracked cells decreases as the number of cell neighbour exchanges increases. However, the number of incorrectly tracked cells remains below 20% throughout the analysed range of neighbour exchanges, and decreases to zero as the number of edge swaps exceeds 10%. The number of untracked cells increases rapidly as the percentage of cell-cell interfaces that are swapped between successive images increases from 5 to 10%. Note that the percentage of cells involved in these neighbour exchanges is larger than the percentage of cell-cell interfaces that are swapped, as an individual T1 transition changes the cell neighbour relations of four cells, and each cell shares multiple inner edges. For example, rearranging 5% of the inner edges of the tissue affects roughly 40% of the cells in the tissue, while rearranging 10% of the tissue edges affects up to 70% of the cells. The number of (correctly or incorrectly) tracked cells drops to zero if the tissue rearranges so much that the neighbourhood of each cell changes; in this case, a first match cannot be found to initialize the MCS construction algorithm. Figure 8 shows the first three of 21 segmented image frames of the lateral epidermis of a stage-11 Drosophila embryo to which the algorithm was applied. During stage 11, gross morphogenetic movements do not occur but the tissue is very active with a large number of proliferations occurring within a short duration, making this a much more challenging tissue on which to perform cell segmentation and tracking than the wing imaginal disc, where many previous efforts have been made [4,13,31,44]. Cell delamination is also more common than in the wing imaginal disc during normal development. This stage of development thus offers a true test of the present algorithm.

Application of the algorithm to in vivo data
The images were taken 5 min apart over a time span of 100 min. These first three images comprise 271, 263 and 263 cells, respectively. Our algorithm tracks 247 cells between the first and second images, 245 cells between the second and third images and 234 cells across all three images. The centroids of cells of previous images are superimposed on the tracking results in figure 8, illustrating that the algorithm successfully tracks cells in situations where it is difficult to match cells between images based on the centroid positions alone. Cells that include only their corresponding centroid from the previous image are coloured in green, while cells that do not include their corresponding centroid from the previous image, and cells that include multiple centroids from the previous image, are coloured in purple. In the first frame, we highlight 'higher-order' junctions (shared by four or more cells) by yellow asterisks. Such junctions occur frequently throughout the dataset.   The data in figure 8 are the first three out of 21 frames. In figure 9, we show the results of the analysis of the full dataset, including all 21 frames. During the period of measurement, the total number of cells increases from 280 to 330 cells, whereas the total number of tracked cells increases from 270 to roughly 310 cells. As the number of cells in the tissue rises, the total number of cell rearrangements increases, whereas the average cell area decreases. Here, the number of cell rearrangements is measured by counting how many cells change their cell neighbour number between consecutive frames. For all frames, the number of tracked cells is lower than the number of cells in the tissue. A visual inspection of the tracked data reveals that the difference between the total number of cells and the number of tracked cells is largely due to cells entering or leaving the field of view. The percentage of cells that our algorithm tracks is lowest (84%) when the rates of cell division and cell rearrangement are highest, which occurs at 70 min. Here, the number of tracked cells decreases because the algorithm is not yet able to resolve division events immediately adjacent to rearrangements as well as multiple adjacent divisions.
As cell rearrangements are one of the most difficult aspects of cell tracking, and our in vivo data exhibit a high frequency of such events, it is natural to ask what percentage of cells are correctly tracked. To estimate this percentage, we compare the results in figure 9, where up to 30% of cells are involved in neighbour exchange between frames in a population of up to 340 cells, with the results shown in figure 7, where in the case of 400 cells and 4% of edges undergoing T1 transitions between frames (corresponding to 30% of cells involved in neighbour exchanges) we find the percentage of correctly tracked cells to be 85%. This provides a lower bound for the success rate of the algorithm on the in vivo frames. When up to 3% of edges undergo T1 transitions (corresponding to 25% of cells in the tissue involved in neighbour exchanges), the success rate of the algorithm is 98%.
The tracking of epithelial in vivo data enables quantitative assessment of dynamic changes in cellular morphology. The tracking results in figure 9 reveals that the analysed section of the epidermis undergoes 60 cell rearrangements per 5 min initially and around 100 cell rearrangements per 5 min at the end of the observed time interval. The average ratio of the maximal area and the minimal area observed for individual cells during the period of measurement is 4.2, indicating that on average cells increase their apical area by a factor of four during mitotic rounding. A total of 18 cell deaths are tracked in the dataset. A striking feature is the level of T1 transitions occurring during this stage of development, even in the absence of gross morphogenetic movements found in earlier or later stages.

Calculation times
To analyse the scaling of the calculation times with tissue size, we repeat the permutation test with tissues of square dimension of varying size on a desktop computer with an Intel i5-6500 T CPU (2.5 GHz) and 8 GB RAM. We find that the calculation times scale subquadratically with cell number ( figure 10). The calculation times for the experimental images analysed in figure 8 vary more widely than for the in silico datasets. For the tracking between the first and second frames in figure 8, the algorithm required 43 s to run, whereas between the second and the third frames the algorithm required 9 s. This is due to differences in the time required to find the first correct mapping; in the first example, 154 cells were searched before the first correct mapping was found, whereas in the second example only 12 cells were searched. This means that the number of cells considered when finding the initial mappings depends on the graph structure of the analysed frames and impacts on the calculation time of the algorithm. In total, analysing all 21 frames of the in vivo data presented in figures 8 and 9 requires 19 min of calculation time.

Discussion
Cell tracking in epithelial sheets has the potential to generate a vast amount of quantitative data to inform our understanding of the contributions of different cellular processes to tissue morphogenesis. However, cell tracking is notoriously difficult, especially for the complex morphogenetic processes that occur as embryogenesis proceeds. Here, we present an algorithm based on MCS detection for the tracking of cells in segmented images of epithelial sheets. Our algorithm successfully tracks cells in in vivo images of the Drosophila embryonic epidermis, a challenging dataset compared to other tissues, as well as in randomly generated in silico datasets, without the need for the adjustment of tissue-specific parameters such as weights for individual terms in a global minimization scheme [14]. The use of in silico data to test our algorithm allows us to analyse the performance of our algorithm for a large range of experimentally observed cell rearrangements and tessellations.
The tracking of cells in in vivo datasets such as presented in figures 8 and 9 provides quantitative insight into tissue homeostasis. Using our algorithm, we measure example quantities that would not be accessible without a robust cell tracking method. The amount of cell rearrangement, the extent of mitotic rounding and the occurrence of cell death in the observed frames each can be used to learn about tissue homeostasis in developing epithelia. Within the analysed dataset, we find a significant number of T1 transitions despite the absence of gross morphogenetic movements. This may be driven by the large number of proliferation events that occur. Further, using in vivo imaging together with our tracking algorithm allows the observation of cell death or cell delamination without the need for fluorescent markers of apoptosis. Future applications of the algorithm to such processes may, for example, provide novel insight to tissue size control in the Drosophila embryonic epidermis [9,45] and can also be adapted to study the dynamics of epithelial wound closure. In this and other systems, cell tracking may enable the observation of cell death due to delamination as opposed to apoptosis [46]. Access to quantification of cell rearrangement and area changes has recently provided insight to wing morphogenesis in Drosophila [43].
Our algorithm is able to track cells that undergo significant movement and neighbour exchanges between frames. For example, we can correctly track cells in tissues where more than 40% of the cells rearrange between successive movie frames (figure 7). In addition, even comparably large gaps in the initial MCS can be filled in during the post-processing step (figures 2 and 8). For example, in the first tracking step in figure 8, only 182 of the 246 tracked cells were identified by the MCS algorithm, and it was possible to track the 64 remaining cells during the post-processing step. For comparison, Heller et al. [13] report 15 cell rearrangements per 1000 cells per hour at an imaging interval of 6 min for their timelapse microscopy data of Drosophila wing imaginal discs. In addition, the experimental data shown in figures 2, 8 and 9 include junctions shared by four or more cells (yellow asterisks in figure 8), while our in silico data include multiple instances of such junctions (figure 6d). Therefore, higher-order junctions, such as multicellular rosettes [47,48], do not pose a challenge to our algorithm.
Our algorithm is able to correctly track cells in all considered test cases. However, on rare occasions, a few cells at the tissue boundary cannot be tracked. It may be possible to adapt the algorithm to track these cells, if this is considered necessary for the application at hand. In the current version of the algorithm, two connections to already tracked cells that are preserved between two time frames are a condition to add a cell-to-cell mapping in the post-processing algorithm. Further analysis of cases where this condition is not fulfilled may reveal ways to relax it.
When generating in silico data to test the algorithm, we used Voronoi tessellations in combination with Lloyd's relaxation to generate data that resemble tissues in the Drosophila wing imaginal disc [33]. We expect the algorithm to perform less well on tissues whose network structure is nearly homogeneous. For example, in an epithelial sheet where cells are arranged in a hexagonal fashion, such as the early Drosophila embryonic epidermis [49] or the late pupal Drosophila wing [50], the local adjacency network of each cell is identical, and hence a network-based tracking algorithm may not be able to distinguish cells. When generating in silico tissues, we use four Lloyd's relaxation steps after Voronoi tessellation. With each Lloyd's relaxation step, the homogeneity of the tissue increases. We were able to successfully repeat all in silico tests on virtual tissues that were generated using up to 14 Lloyd's relaxation steps. Hence, we expect the algorithm to be suitable for tissues that can be well described with 14 or fewer Lloyd's relaxation steps, such as the chick neural tube embryonic epithelium, or the Drosophila eye disc [33]. The algorithm relies on being able to generate polygonal tessellations from segmented video microscopy data. In particular, all in silico tests we conducted consider tissues where each cell has at least three neighbours. Conceptually, it would be possible to apply the algorithm to tissues in which individual cells may have only two neighbours, although such examples have not been included in the present analysis.
In microscopy videos including division events, we expect the algorithm to perform well in tissues in which no adjacent divisions occur between successive movie frames, and in which cells adjacent to the dividing cell do not undergo rearrangements before the next frame is captured. Our algorithm is designed to identify mother and daughter cells of a division event by establishing the bordering cells that gain an edge during the division event. In the case of two adjacent divisions, and if cells adjacent to a division event gain edges due to cell rearrangements, the dividing cell cannot be correctly identified. An example of a typical tracking error for two adjacent divisions is shown in figure 11. In cases where the division resolution step fails, our Python implementation returns all tracked cells of the post-processing step, and gives a warning that the division has not been resolved. In these cases, manual correction methods could be used for incorrectly tracked cells in the vicinity of division events.
The parameters of the algorithm are chosen to maximize its robustness and avoid the necessity to adjust the parameters to individual applications. For example, the cutoff length d max that determines the distance below which two cells in consecutive movie frames are considered mappable to each other was chosen to be 10 times the average cell length in the tissue, which is significantly larger than the movement that is to be expected between consecutive frames of a live-imaging microscopy video. However, parameter adjustments may be possible for individual applications in order to decrease the algorithm calculation times. For example, the size of the extended neighbourhood considered in the initial step or the iterative extension could be reduced to include only nearest neighbours instead of nearest neighbours and second nearest neighbours in case the tissue is sufficiently heterogeneous. Similarly, one might decrease the d max for possible cell pairings if the cell positions are not expected to vary significantly between time frames.
Adjustments may be possible to extend the applicability of the algorithm to a wider range of tissues. For example, instead of automatic detection of the initial seeds for the MCS detection algorithm, a small set of seeds could be manually supplied to guide the tracking. This should improve the performance of the algorithm on homogeneous tissues. In such cases, irregular boundaries may also help to aid the initial seeding. During the construction of the MCS, non-adjacent cells are considered for addition to the MCS whenever the extension of the MCS by adjacent cells is not possible. An alternative option to extend an intermediate MCS may be to repeat the initial seeding algorithm.
In this work, we have sought to keep geometrical input to the algorithm to a minimum. Cases where geometric data are taken into account comprise division events where one of the daughter cells is four-or three-sided, because in these cases we are not able to make a decision on which cell is the second daughter cell based on network adjacency alone. If future applications reveal cases where the algorithm performs poorly due to a large number of cell neighbour exchanges or high degree of tissue homogeneity, then it may be possible to construct algorithms that combine information on the network topology with data on cell shapes, cell positions and cell movements to improve performance. For example, information on network topology could be integrated into previous algorithms that minimize differences between geometric properties of cells, such as cell size and location [14], with information about network connectivity.
In cell tracking applications, the scaling of the algorithm with tissue size is crucial. Potential applications range from systems of 30 cells (Drosophila embryonic epidermal P compartments [9]), to 10 000 cells (Drosophila imaginal wing disc [31]; the wing pouch has about 3000 cells [40]). Calculation times in the presented algorithm scale subquadratically with cell number, making it suitable for applications of varying sizes. For example, extrapolating the data in figure 10, a tissue of 10 000 cells could be tracked across two frames within 10 min. The scaling of the algorithm is polynomial despite the fact that it is based on MCS detection, which is known to scale exponentially in the general case. MCS detection has a wide range of research applications, including protein interaction networks [51,52] and finding the binding sites of chemical structures [53]. Our approach of reducing the MCS search to a localized search may have applications in other areas where the networks are inherently planar.
Our algorithm is designed to track cells in segmented microscopy videos of epithelial sheets in two dimensions. However, it may be possible to apply the algorithm to datasets of epithelial sheets that are embedded in a threedimensional environment, such as the Drosophila imaginal wing disc [4], or the Drosophila embryonic epidermis [6,9], including tissues that can be mapped onto a cylinder or ellipsoid, such as the mouse visceral endoderm [48].
A large number of cell tracking algorithms have been developed for various applications [10 -27]. Further efforts are required to compare these algorithms with our own, and to identify the algorithm best suited for an individual dataset. In the cell-tracking challenge [54], the authors provide microscopy videos from a variety of in vitro cell cultures, including, for example, mouse embryonic stem cells and human squamous lung carcinoma cells, together with ground truth segmentation and tracking data as benchmarks for cell tracking and segmentation algorithms. However, many of the published algorithms above have not yet been applied to the challenge, and benchmark datasets for epithelial sheets are not currently available. The fully segmented dataset published within the MCSTracker project can provide a benchmark for future epithelial cell tracking applications. In [55], in silico datasets are used as benchmarking datasets for particle tracking algorithms.
The proposed algorithm provides a two-dimensional tracking solution specialized for epithelial sheets that attempts to maximize the information that can be gained from the packing that is typical to such tissues. It may, however, be possible to extend this algorithm to applications of cell tracking, where cells are not physically connected by constructing adjacency networks from Voronoi tessellations that use the cell locations as seeds. We hope that, as segmentation tools are developed further, the combination of our algorithm with these tools will lead to further insights into cellular behaviour in epithelial tissues.
Data accessibility. The datasets and code supporting this article are publicly available at https://github.com/kursawe/MCSTracker. Authors' contributions. J.K., R.B., J.J.Z., R.E.B. and A.G.F. conceived the study and contributed to the manuscript. J.K. developed the algorithm and carried out the data segmentation and performance analysis. J.K., R.E.B. and A.G.F. drafted the manuscript. J.J.Z. provided microscopy data. All authors gave final approval for publication.
Competing interests. We declare we have no competing interests. Funding. J.K. acknowledges funding from the Engineering and Physical Sciences Research Council through a studentship. R.B. acknowledges funding from ANR grant ANR-16-CE23-0003-01. J.Z. acknowledges funding support from the National Science Foundation (awards CBET-1553826 and CBET-1403887).