Stochastic Neural Networks for Automatic Cell Tracking in Microscopy Image Sequences of Bacterial Colonies

Our work targets automated analysis to quantify the growth dynamics of a population of bacilliform bacteria. We propose an innovative approach to frame-sequence tracking of deformable-cell motion by the automated minimization of a new, specific cost functional. This minimization is implemented by dedicated Boltzmann machines (stochastic recurrent neural networks). Automated detection of cell divisions is handled similarly by successive minimizations of two cost functions, alternating the identification of children pairs and parent identification. We validate the proposed automatic cell tracking algorithm using (i) recordings of simulated cell colonies that closely mimic the growth dynamics of E. coli in microfluidic traps and (ii) real data. On a batch of 1100 simulated image frames, cell registration accuracies per frame ranged from 94.5% to 100%, with a high average. Our initial tests using experimental image sequences (i.e., real data) of E. coli colonies also yield convincing results, with a registration accuracy ranging from 90% to 100%.


Introduction
Technology advances have led to increasing magnitudes of data generation with increasing levels of precision [26,90]. However, data generation presently far outpaces data analysis and drives the requirement for analyzing such large-scale data sets with automated tools [15,51,97]. The main goal of the present work is to develop computational methods for an automated analysis of microscopy image sequences of colonies of E. coli growing in a single layer. Such recordings can be obtained from colonies growing in microfluidic devices, and they provide a detailed view of individual cell-growth dynamics as well as population-level, inter-cellular mechanical and chemical interactions [18,30,68].
However, to understand both variability and lineage-based correlations in cellular response to environmental factors and signals from other cells requires the tracking of large numbers of individual cells across many generations. This can be challenging, as large cell numbers tightly packed in microfluidic devices can compromise spatial resolution, and toxicity effects can place limits on the temporal resolution of the recordings [34,42]. One approach to better understand and control the behavior of these bacterial colonies is to develop computational methods that capture the dynamics of gene networks within single cells [18,49,106]. For these methods to have a practical impact, one ultimately has to fit the models to the data, which allows us to infer hidden parameters (i.e., characteristics of the behavior of cells that cannot be measured directly). Image analysis and pattern recognition techniques for biological imaging data [27,45,70], like the methods developed in the present work, can be used to track lineages and thus automatically infer how gene expression varies over time. These methods can serve as an indispensable tool to extract information to fit and validate both coarse and detailed models of bacterial population, thus allowing us to infer model parameters from recordings. We show five frames out of a total of 150 frames of an image sequence showing the dynamics of E. coli in a microfluidic device [6] (real laboratory image data). These cells are are about 1 µm in diameter and on average 3 µm in length, and they divide about every 30 minutes. The original images exported from the microscope are 0.11 µm/pixel. We report results for these real datasets in Sec. 4.
Here we describe an algorithm that provides quantitative information about the population dynamics, including the life cycle and lineage of cells within a population from recordings of cells growing in a mono-layer. A typical sequence of frames of cells growing in a microfluidic trap is shown in Fig. 1. We describe the design and validation of algorithms for tracking individual cells in sequences of such images [5,49,57]. After segmentation of individual image frames to identify each cell, tracking individual cells from frame to frame is a combinatorial problem. To solve this problem we take into account the unknown cell growth, cell motion, and cell divisions that occur between frames. Segmentation and tracking are complicated by imaging noise and artifacts, overlap of bacteria, similarity of important cell characteristics across the population (shape; length; and diameter), tight packing of bacteria, and large interframe durations resulting in significant cell motion, and up to a 30% increase in individual cell volume.
1.1. Related Work. The present work focusses on tracking E. coli in time series of images. A comparison of different cell-tracking algorithms can be found in [37,99]. Multi-object tracking in video sequences and object recognition in time series of images is a challenging task that arises in numerous applications in computer vision [66,109]. In(biomedical) image processing, motion tracking is often referred to as "image registration" [37,59,63,65,64] or "optical flow " [23,40,31,62]. Typical solutions used in the defense industry, for instance, track small numbers of fast moving targets by image sequence analysis at pixel levels and use sophisticated reconstruction of the optical flow, combined with real time segmentation, and quick combinatoric exploration at each image frame. Initially, we did implement several well known algorithms for reconstruction of the optical flow but the results we obtained were not satisfactory due to long interframe times and high noise levels. Moreover, we are not interested in tracking individual pixels but rather cells (i.e., rod-shaped, deformable shapes), while recognizing events of cell division and recording cell lineage. Consequently, we decided to first segment each image frame to isolate each cell, and then to match cells between successive frames.
As for the problem at hand, one approach proposed in prior work to simplify the tracking task is to make the experimental setup more rigid by confining individual cell lineages to small tubes; the associated microfluidic device is called a "mother machine" [17,44,60,76,87,94]. The microfluidic devices we consider here yield more complicated data as cells are allowed to move and multiply freely in two dimensions (constrained to a mono-layer). We refer to Fig. 1 for a typical sequence of experimental images considered in the present work.
Turning to methods that work on more complex biological cell imaging data, we can distinguish different classes of tracking methods. "Model-based evolution methods" operate on the image intensities directly. They rely on particle filters [7,75,93] or active contour models [47,53,105,108,88]. These methods work well if the cells are not tightly packed. However, they may lead to erroneous results if the cells are close together, the inter-cellular boundaries are blurry, or the cells move significantly. Our work belongs to another class-the so called "detection-and-association methods" [16,22,20,46,80,85,98,104,43,110,91], which first detect cells in each frame and then solve the tracking problem/association task across successive frames. (We note that the segmentation and tracking of cells does not necessarily need to be implemented in two distinct steps. In many image sequence analyses, implementing these two steps jointly can be beneficial [7,38,78,80,77,110]. However, for the clarity of exposition and easier implementation of our new tracking technique, we present these steps separately.) Doing so necessitates the segmentation of cells within individual frames. We refer to [100] for an overview of cell segmentation approaches. Deep learning strategies have been widely used for this task [4,35,61,72,77,85,84,96,97,110]. We consider a framework based on convolutional neural networks (CNNs). Others have also used CNNs for cell segmentation [3,61,74,83,84]. We omit a detailed discussion of our segmentation approach within the main body of this paper, as we do not view it as our main contribution (see Sec. 1.2). However, the interested reader is referred to Sec. D for some insights. To solve the tracking problem after the cell detection, many of the methods cited above use hand-crafted association scores based on the proximity of the cells and shape similarity measures [46,22,98,110]. We follow this approach here. We note that we not only consider local association scores between cells but also include measures for the integrity of a cell's neighborhood (i.e., "context information").
Our method is tailored for tracking cells in tightly packed colonies of rod-shaped E. coli bacteria. This problem has been considered previously [16,80,97,104]. However, we are not aware of any large-scale datasets that provide ground truth tracking data for these types of recordings, but note that there are community efforts for providing a framework for testing cell tracking algorithms [67,99] (see, e.g., [7,58]). 1 Works that consider these data are for example [7,61,74,78,77,110]. The cells in this dataset have significantly different characteristics compared to those considered in the present work. As we describe below, our approach is based on distinct characteristics of the bacteria cells and, consequently, does not directly apply to these data. Therefore, we have developed our own validation and calibration framework (see Sec. 2.1).
Standard graph matching algorithms (see, e.g., [103]) do not directly apply to our problem. Indeed, a fundamental complication is that cells can divide between successive images. Hence, each assignment from one frame to its successor is not a one-to-one mapping but a one-to-many correspondence. More advanced graph matching strategies are described in [79,36]. Graph-based matching strategies for cell-tracking that are somewhat related to our approach are described in [28,56,54,58,55]. Like the methods mentioned above, they consider various association scores for tracking. Individual cells are represented as nodes, and neighbors are connected through edges. Our approach also introduces cost terms for structural matching of local neighborhoods by specific scoring for single nodes, pairs of nodes, and triplets of nodes, after a (modified) Delaunay triangulation. By using a graph-like structure, cell divisions can be identified by detecting changes in the topology of the graph [56,54]. We tested a similar strategy, but came to the conclusion that we cannot reliably construct neighborhood networks between frames for which topology changes only occur due to cell division; the main issue we observed is that the significant motion of cells between frames can introduce additional topology changes in our neighborhood structure. Consequently, we decided to relax these assumptions. [103,102,81,48] implement multi-target tracking in videos by stochastic models based on random finite set densities and variants thereof. The fit to the data is based on Gibbs sampling to maximize the posterior likelihood. A key challenge of these approaches is the estimation of an adequate finite number of Gibbs sampling iterations when one computes posterior distributions. Most Gibbs samplers are ergodic Markov chains on a finite but huge state spaces, so that their natural exponential rate of convergence is not a practically reassuring feature.
As mentioned above, some recent works jointly solve the tracking and segmentation problem [7,38,78,80,77,110]. Contrary to observations we have made in our data, these approaches rely (with the exception of [80]) on the fact that the tracking problem is inherent to the segmentation problem ("tracking-by-detection methods" [110]; see also [97]). That is, the key assumption made by many of these algorithms is that cells belonging to the same lineage overlap across frames (see also [20]). In this case, cell-overlap can serve as a good proxy for cell-tracking [110]. We note that in our data we cannot guarantee that the frame rate is sufficiently high for this assumption to hold. [78,74,38] exploited machine learning techniques for segmentation and motion tracking. One key challenge here is to provide adequate training data for these methods to be successful. Here, we describe simulation-based techniques that can be extended to produce training data, which we use for parameter tuning [106,107].
The works that are most similar to ours are [16,80,104]. They perform a local search to identify the best cell-tracking candidates across frames. One key difference across these works are the matching criteria. Moreover, [16,80] employ a local greedy-search, whereas we consider stochastic neural network dynamics for optimization. [104] construct score matrices within a score based neighborhood tracking method; an integer programming method is used to generate frameto-frame correspondences between cells and the lineage map. Other approaches that consider linear programming to maximize an association score function for cell tracking can be found in [21,20,110].
As we have mentioned in the abstract we obtain a tracking accuracy that ranges from 90% to 100%, respectively. Overall, our method is competitive with existing approaches: For example, [16] report a tracking accuracy of up to 97% for data that is similar to ours, while [28] report a tracking accuracy (spatial, temporal, and cell division detection) at the order of 95% (between about 93% and 98%, respectively). The second group also reports results for their prior approach [56], with an accuracy at the order of 90% (ranging from about 87% to 92%, respectively). Accuracies reported in [55] range from about 92% to 97%, respectively. This work also includes a comparison to one of their earlier approaches [54] with an accuracy of up to 85% and 89% if the datasets are prealigned. We note that the data considered in [28,55,54,56] is quite different from ours. [7,58,77,78,110] consider the data from the cell tracking challenge [67,99] to evaluate the performance of their methods. As in the previously mentioned work, this data is again quite different from ours. To evaluate the performance of the methodology the so called acyclic oriented graph matching measure [69] is considered. We refer to the webpage of the cell tracking challenge for details on the evaluation metrics (see http://celltrackingchallenge.net/evaluation-methodology). Based on these, the reported tracking scores are between 0.873 and 0.902 [7], 0.901 and 1.00 [58], 0.950 and 0.987 [110], 0.788 and 0.982 [77] and 0.765 and 0.915 [78] depending on the considered data set, respectively.

1.2.
Contributions. For image segmentation, we first apply two well-known, powerful variational segmentation algorithms to generate a large training set of correctly delineated single cells. We can then train a CNN dedicated to segmenting out each single cell. Using a CNN significantly reduces the runtime of our computational framework for cell identification. The frame-to-frame tracking of individual cells in tightly packed colonies is a significantly more challenging task, and is hence the main topic discussed in the present work. We develop a set of innovative automatic cell tracking algorithms based on the successive minimization of three dedicated cost functionals. For each pair of successive image frames, minimizing these cost functionals over all potential cell registration mappings poses significant computational and mathematical challenges. Standard gradient descent algorithms are inefficient for these discrete and highly combinatorial minimization problems. Instead, we implement the stochastic neural network dynamics of BMs, with architectures and energy functions tailored to effectively solve our combinatorial tracking problem. Our major contributions are: i) The design of a multi-stage cell tracking algorithm, that starts with a parentchildren pairing step, followed by removal of identified parent-children triplets, and concludes with a cell to cell registration step. ii) The design of dedicated BM architectures, with several energy functions, respectively, minimized by true parent-children pairing and by true cell-to-cell registration. Energy minimizations are then implemented by simulation of BM stochastic dynamics.
iii) The development of automatic algorithms for the estimation of unknown weight parameters of our BM energy functions, using convex-concave programming tools [2,32,89]. iv) The evaluation of our methodology on synthetic and real image sequences of cell colonies. The massive effort involved in human expert annotation of cell colony recordings limits the availability of "ground truth tracking" data for dense bacterial colonies. We therefore first validated the accuracy of our cell tracking algorithms on recordings of simulated cell colonies, generated by the dedicated cell colony simulation software [106,107]. This provided us with ground truth frame-by-frame registration for cell lineages, enabling us to validate our methodology.

Datasets
Below we introduce the datasets used to evaluate the performance of the proposed methodology. The synthetic data is described in Sec. 2.1. The experimental data (real imaging data) is described in Sec. 2.2.
2.1. Synthetic Videos of Simulated Cell Colonies. To validate our cell tracking algorithms, we consider simulated image sequences of dense cell populations. We refer to [106,107] for a detailed description of this mathematical model and its implementation. 2 The simulated cell colony dynamics are driven by an agent based model [106,107], which emulates live colonies of growing, moving, and dividing rod-like E. coli cells in a 2D microfluidic trap environment. Between two successive frames J, J + , cells are allowed to move until they nearly bump into each other, and to grow at multiplicative rate denoted g.rate with an average value of 1.05 per minute (plus/minus a small random perturbation).
The cells are modeled as 2D spherocylinders of constant 1 µm width. Each cell grew exponentially in length with a doubling time of 20 minutes. To prevent division synchronization across the population when a mother cell of length L div divides, the two daughter cells are assigned random birth lengths L 0 (b 1 ) = L 1 = δL div and L 0 (b 2 ) =: L 2 = (1 − δ)L div , where δ > 0 is a random number sampled independently at each division from a uniform distribution on [0.45, 0.55]. Consequently, a bacterial cell b of length L div divides into two cells b 1 and b 2 , their lengths L 1 , L 2 satisfy L 1 + L 2 = L div and L i /L div , i = 1, 2, is a random number. The cells have a length of approximately 2µm after division and 4µm right before division. We refer to [107] for additional details. The simulation keeps track of cell lineage, cell size, and cell location (among other parameters). The main output of each such simulation considered here is a binary image sequence of the cell colony with a fixed interframe duration. Each such synthetic image sequence is used as the sole input to our cell tracking Table 1. Benchmark datasets. To test the tracking software we consider simulated data. We have generated data of varying complexity with different interframe durations. We note that we also consider these data to train our algorithms for tracking cells. We report the label for each dataset, the interframe duration, as well as the number of frames generated. We set the cell growth factor to g.rate = 1.05 per minute. We refer to the text for details about how these data have been generated.

Label
Interframe Duration Number of Frames algorithm. The remaining meta-data generated by the simulations are only used as ground truth to evaluate the performance of our tracking algorithms. We consider several several benchmark datasets of synthetic image sequences of simulated cell colonies of different complexity. We refer to these benchmarks as BENCH1 (500 frames), BENCH2 (300 frames), BENCH3 (300 frames), and BENCH6 (100 frames), with an interframe duration of 1, 2, 3 and 6 minutes, respectively. Notice that there is no explicit noise on the growth rate. However, due to crowding of cells the growth rate will vary from cell-to-cell. The generated binary images are of size 600 × 600 pixels. We summarize these benchmarks in Tab. 1. The associated image sequences involve between 100 up to 500 frames, respectively. In Fig. 2 we display an example of two simulated consecutive frames separated by 1 minute. To simplify our presentation and validation tests, we control our simulations to make sure that cells will not exit the region of interest from one frame to the next, and we exclude cells that are only partially visible in the current frames.
2.2. Laboratory Image Sequences (Real Biological Data). We also verify the performance of our approach on real datasets of E. coli bacteria. These bacteria are about 1 µm in diameter and on average 3 µm in length, and they divide about every 30 minutes. The original images exported from the microscope are 0.11µm/pixel. The microscopy experimental data was obtained using JS006 [95] (BW25113 ∆araC ∆lacI) E. coli strains containing a plasmid constitutively expressing yellow or cyan fluorescent protein (sfyfp or sfcfp) for identification. The plasmid also contains an ampicillin resistance gene and p15A origin. These cells were grown overnight in LB medium with 100 µg/mL ampicillin for 18 hours. These cultures were diluted in the morning into 1/1000 into 50 mL fresh LB with 100 µg/ml ampicillin and grown for 3 hours until they reached an OD600 of about 0.3. The cells were then concentrated by centrifuging 30 ml of culture at 2000 g for 5 minutes and then resuspending in 10 ml of fresh LB. The concentrated culture was loaded into a hallway microfluidic device prewarmed and flushed with 0.1% (v/v) Tween-20 [29]. In the microfluidic device, the cells were provided with continuous fresh LB with 100 µg/ml ampicillin and 0.075% (v/v) Tween-20. The microfluidic device was placed onto an 60× oil objective and imaged every 6 min at phase contrast, YFP, and CFP filter settings using an inverted fluorescence microscope. We show a representative dataset in Fig. 1. 2.3. Cell Characteristics. Next, we discuss characteristics of the E. coli bacteria important for our tracking algorithm.
Cell Geometry. In accordance with the dynamics of bacterial colonies in microfluidic traps, the dynamic simulation software generates colonies of rod-shaped bacteria. Cell shapes can be approximated by long and thin ellipsoids, which are geometrically well identified by their center, their long axis, and the two endpoints of this long axis. The center c(b) is the centroid of all pixels belonging to cell b.
For the synthetic images of size 600 × 600 that we considered (see Sec. 2.1), we take ρ = 80 pixels. We write b 1 ∼ b 2 for short, whenever b 1 , b 2 are neighbors (i.e, satisfy the three conditions identified above).
Cell Motion. Let J, J + denote two successive images (i.e., frames). Denote B = B(J), B + = B(J + ) the associated sets of cells. Superpose temporarily the images J and J + so that they then have the same center pixel. Any cell b ∈ B, which does not divide in the interframe J → J + , becomes a cell b + in image J + . The "motion vector " of cell b from frame J to J + is then defined by v(b) = c(b + ) − c(b). If the cell b does divide between J and J + , denote b div the last position reached by cell b at the time of cell division, and define similarly the motion v(b) = c(b div ) − c(b). In our experimental recordings of real bacterial colonies with interframe duration 6 min, there is a fixed number w > 0 such that v(b) ≤ w/2 for all cells b ∈ B(J) for all pairs J, J + . In particular, we observed that for real image sequences, w = 100 pixels is an adequate choice. Consequently, we select w = 100 pixels for all simulated image sequences of BENCH6. For BENCH1 we select w = 45 pixels, again based on a comparison with real experimental recordings. Overall, the metaparameter w is assumed to be a fixed number and to be known, since w/2 is an observable upper bound for the cell motion norm for a particular image sequence of a lab experiment.
Target Window. Recall that J, J+ are temporarily superposed. Let U (b) ⊂ J + be a square window of width w, with the same center as cell b. The target window W (b) is the set of all cells in B + having their centers in U (b). Since v(b) ≤ w/2, the cell b + must belong to the target window W (b) ⊂ B + .

Methodology
3.1. Registration Mappings. Next we discuss our assumptions on a valid registration mapping that establishes cell-to-cell correspondences between two frames. Let J, J + denote two successive images, with cell sets B and B + , respectively. As above, we let N = card(B), and N + = card(B + ).
Our goal is to track each cell from J to J + . For each cell b ∈ B, there exist three possible evolutions between J and J + : Case 1:: Cell b ∈ B did not divide in the interframe J → J + , and has become a cell f (b) ∈ B + ; that is, f (b) has grown and moved during the interframe time interval. Case 2:: Cell b ∈ B divided between J and J + , and generated two children To simplify our exposition, we ignore Case 3. We discuss Case 3 in greater detail in the conclusions in Sec. 5. Consequently, a valid (true) registration mapping f will take values in the set

3.2.
Calibration of Cost Function Weights. With the notation we introduced, fix any two finite sets A, A + . Let G := {g : A → A + } be the set of all mappings g : A → A + . Fix m penalty functions pen k (g) ≥ 0, k = 1, . . . , m. Let g * ∈ G be the ground truth mapping we want to discover through minimization in g of some given cost function COST(g) defined by the linear combination of the penalty functions pen k (g), the contributions of which are controlled by the cost function weights λ k > 0. In this section, we present a generic weight calibration algorithm, extending a technique introduced and applied in [10,11] for Markov random fields based image analysis.
The cost function must perform well (with the same weights) for hundreds of pairs of (synthetic) images J, J + . We consider one such synthetic pair for which the ground truth registration mapping f ∈ G is known, and use it to compute an adequate set of weights, which will then be used on all other synthetic pairs J, J + . Notice, that for experimental recordings of real cell colonies, no ground truth registration mappings f is available. In this case, f should be replaced by a set of user constructed, correct partial mappings defined on small subsets of A. The proposed weight calibration algorithm will also work in those situations.
We now show how knowing one ground truth mapping f can be used to derive the best feasible weights ensuring that f should be a plausible minimizer of the cost functional COST(g) over g ∈ G.
Let PEN(g) = [pen 1 (g), . . . , pen m (g)] be the vector of m penalties for any mapping g ∈ G. Let Λ = [λ 1 , . . . , λ m ] be the weights vectors. Then, COST(g) = Λ, PEN(g) . Replacing g by another mapping h = g induces the penalty changes ∆ PEN g,h = PEN(h) − PEN(g) and the cost change ∆ COST(g, h) = Λ, ∆ PEN g,h . Now, fix any known ground truth mapping f ∈ G. We want f to be a minimizer of COST, so we should have ∆ COST(f, f ) ≥ 0 for all modifications f → f ∈ G. For each a ∈ A, select an arbitrary s(a) ∈ W (a) (where W (a) is the target window for cell a; see Sec. 2.3), to define a new mapping f = f a from A to A + by f a (a) = s(a), and f a (x) ≡ f (x) for all x = a. Since f is a minimizer of COST, this single point modification f → f a must generate the following cost increase Then, the positive vector Λ ∈ R m , Λ 0, should verify the set of linear constraints Λ, V a ≥ 0 for all a ∈ A. There may be too many such linear constraints. Consequently, we relax these constraints by introducing a vector y = [y(a)] ∈ R card(A) , y 0, of slack variables y(a) ≥ 0 indexed by all the a ∈ A. (In optimization slack variables are introduced as additional unknowns to transform inequality constraints to an equality constraint and a non-negativity constraint on the slack variables.) We require the unknown positive vector Λ and the slack variables vector y to verify the system of linear constraints: The normalizing constant 1000 can be arbitrarily changed by rescaling. We seek high positive values for ∆ COST(f, f a ) and small L 1 -norm for the slack variable vector y. So, we will seek two vectors Λ ∈ R m and y ∈ R card(A) solving the following convex-concave minimization problem, where γ > 0 is a user selected (large) meta parameter: To numerically solve the constrained minimization problem (2), we use the libraries CVXPY and DCCP (disciplined convexconcave programming) [2,32,89]. DCCP is a package for convex-concave programming designed to solve non-convex problems. 3 It can handle objective functions and constraints with any known curvature as defined by the rules of disciplined convex programming [24]. We give examples of numerically computed weight vectors Λ below. The computing time was less than 30 seconds for the data that we have prepared. For simplicity, we just considered one step changes in our computations, which make the overlap penalty weak. To increase the accuracy of the model it is possible to consider a larger number of samples (i.e., multi-step changes). Note that the solutions Λ of (2) are of course not unique, even after normalization by rescaling.

3.3.
Cell Divisions and Parent-Children Short Lineages. Next we discuss how we tackle the assignment problem when cells divide.
3.3.1. Cell Divisions. We now outline a cost function based methodology to detect cell divisions.
The first step will be to seek the most likely parent for each potential pair of children cells. For our video recordings of actual cell populations, during any interframe, we may have n out cells exiting the field of view and n in cells entering it, so that | card(trueCH ) − DIV | may be of the order of n in + n out . To take this into account, we relax the constraint in (3) as follows where REL is a fixed bound estimated from our experiments. For simplicity, we have restricted our methodology to the situation where n in and n out are always 0. But even in that case, there was a computational advantage to using the slightly relaxed constraint (4) with REL = 1.

Most Likely Parent Cell for a Given Children Pair.
For successive images J, J + with 1 minute interframe, define the set PCH of plausible children pairs by where the threshold τ > 0 is user selected and fixed for the whole benchmark set BENCH1 of synthetic image sequences.
To evaluate if a pair of cells (b 1 , b 2 ) ∈ PCH can qualify as a pair of children generated by division of a parent cell b ∈ B, we now quantify the geometric distortion between b and (b 1 , b 2 ). Cell division of b into b 1 , b 2 ∈ B + occurs with small motions of b 1 , b 2 . During the short interframe duration the initial centers c 1 , c 2 of b 1 , b 2 in image J move by at most w/2 pixels each (see Sec. 2.3), and their initial distance to the center c of b is roughly at Define the set SHLIN of potential short lineages as the set all triplets verifying the preceding constraint (6). For each potential lineage (b, b 1 , b 2 ) ∈ SHLIN , define three terms penalizing the geometric distortions between a parent b ∈ B and a pair of children (b 1 , b 2 ) ∈ PCH by the following formulas, where we denote c, c 1 , c 2 , the centers of cells b, b 1 , b 2 and A, A 1 , A 2 their long axes, respectively: (i) center distortion cen Here, angle denotes "angles between non-oriented straight lines," with a range from 0 to π/2. Introduce three positive weights λ cen , λ siz , λ ang (to be estimated), and for every short lineage For each plausible pair of children (b 1 , b 2 ) ∈ PCH , we will compute the most likely parent cell (7) over all b ∈ B, as summarized by the formula To force this minimization to yield a reliable estimate of b * = parent(b 1 , b 2 ) for most true pairs of children (b 1 , b 2 ), we calibrate the weights λ j , j ∈ {cen, siz, ang} by the algorithm outlined in Sec. 3.2, using as "ground truth set" a fairly small set of visually identified true short lineages (parent, children). For fixed (b 1 , b 2 ), the set of potential parent cells b ∈ B has very small size due to the constraint (6). Hence, brute force minimization of the functional dist is still a greedy minimization in the sense that other soft constraints introduced further on are not taken into consideration during this preliminary fast computation of b * .

3.3.3.
Penalties to Enforce Adequate Parent-Children Links. Any true pair of children cells pch = (b 1 , b 2 ) should belong to PCH , but must also verify lineage and geometric constraints which we now enforce via several penalties. Note that the new penalties introduced here are fully distinct from the three penalties specified above to define dist(b, b 1 , b 2 ).
"Lineage" Penalty. Valid children pairs (b 1 , b 2 ) ∈ PCH should be correctly matchable with their most likely parent cell b * = parent(b 1 , b 2 ) (see (8)). So, we define the "lineage" penalty Notice that the computation of lin(b 1 , b 2 ) is quite fast. "Gap" Penalty. Denote tips(b) the set of two endpoints of any cell b. For any pair pch = (b 1 , b 2 ) ∈ PCH , define endpoints x 1 ∈ tips(b 1 ), x 2 ∈ tips(b 2 ) and the gap penalty gap(b 1 , b 2 ) by (9) gap "Dev " Penalty. For rod-shaped bacteria, a true pair (b 1 , b 2 ) ∈ PCH of just born children must have a small gap(b 1 , b 2 ) = x 1 − x 2 and roughly aligned cells b 1 and b 2 . For (b 1 , b 2 ) ∈ PCH , we quantify the deviation from alignment dev(b 1 , b 2 ) as follows. Let x 1 , x 2 be the closest endpoints of b 1 , b 2 (see (9)). Let str 12 be the straight line linking the centers c 1 , "Ratio" Penalty. True children pairs must have nearly equal lengths. So, for (b 1 , b 2 ) ∈ PCH with lengths L 1 , L 2 we define the length ratio penalty by "Rank " penalty. Let L min be the minimum cell length over all cells in B + . In B + , children Given two successive images J, J + , we seek the set X = trueCH of true children pairs in B + ×B + which is an unknown subset of PCH . In Sec. 3.5 below we replace X by its indicator function z and we build a cost function E(z) which should be nearly minimized when z is close to the indicator of trueCH . A key term of E(z) will be a weighted linear combination of the penalty functions {lin, gap, dev, ratio, rank}. Since these penalties are different from those introduced in Sec. 3.3.2, we estimate their weights in the cost function E(z) by the algorithm outlined in Sec. 3.2. The minimization of E(z) will be implemented by simulations of a BM with energy function E(z). We present these stochastic neural networks in the next section.

Generic Boltzmann Machines (BMs).
Minimization of our main cost functionals is a heavily combinatorial task, since the unknown variable is a mapping between two finite sets of sizes ranging from 80 to 120. To handle these minimizations, we use BMs originally introduced by Hinton et al. (see [1,39]). Indeed, these recurrent stochastic neural networks can efficiently emulate some forms of simulated annealing.
Each BM implemented here is a network BM = {U 1 , . . . , U N } of N stochastic neurons U j . In the BM context, the time t = 0, 1, 2, . . . is discretized and represents the number of steps in a Markov chain, where the successive updates Z(t) → Z(t + 1) of the BM configuration Z(t) are analogous to the steps of a Gibbs sampler. The configuration Z(t) = {Z 1 (t), . . . , Z N (t)} of the whole network BM at time t is defined by the random states Z j (t) of all neuron U j . Each Z j (t) belongs to a fixed finite set W (j). Hence Z(t) belongs to the configurations set CONF = W (1) × · · · × W (N ).
Neurons interactivity is specified by a finite set CLQ of cliques. Each clique K is a subset of S = {1, . . . , N }. During configuration updates Z(t) → Z(t + 1), neurons may interact only if they are in the same clique. Here, all cliques K are of small sizes 1, or 2, or 3.
For each clique K, one specifies an energy function J K (z) defined for all z ∈ CONF , with J K (z) depending only on the z j such that j ∈ K. The full energy E(z) of configuration z is then defined by The BM stochastic dynamics Z(t) → Z(t + 1) is driven by the energy function E(z), and by a fixed decreasing sequence of virtual temperatures Temp(t) > 0, tending slowly to 0 as t → ∞. Here we use standard temperature schemes of the form Temp(t) ≡ cη t with fixed c > 0 and slow decay rate 0.99 < η < 1.
We have implemented the classical "asynchronous" BM dynamics. At each time t, only one random neuron U j may modify its state, after reading the states of all neurons belonging to cliques containing U j . A much faster alternative, implementable on GPUs, is the "synchronous" BM dynamics, where at each time t roughly 50% of all neurons may simultaneously modify their states (see [9,8,13]). The detailed BM dynamics is presented in the appendix (see Sec. A).
When the virtual temperatures Temp(t) decrease slowly enough to 0, the energy E(Z(t)) converges in probability to a local minimum of the BM energy E(z) over all configurations z ∈ CONF .
3.5. Optimized Set of Parent-Children Triplets. Next, we formulate the search for bona fide parent-children triplets as an optimization problem. For brevity, this outline is restricted to situations where (3) holds, as is the case for our synthetic image data. Simple modifications extend this approach to the relaxed constraint (4) We now define a binary BM constituted by m binary stochastic neurons U j , j = 1 . . . m. At time t = 0, 1, 2, . . ., each U j has a random binary valued state Z j (t) = 1 or 0. The random configuration Let SUB be the set of all subsets of PCH . Each configuration z ∈ CONF is the indicator function of a subset sub(z) of PCH . We view each sub(z) ∈ SUB as a possible estimate for the unknown set trueCH ⊂ B + × B + of true children pairs (b 1 , b 2 ). For each potential estimate sub(z) of trueCH , the "lack of quality" of the estimate sub(z) will be penalized by the energy function E(z) ≥ 0 of our binary BM. We now specify the energy E(z) for all z ∈ CONF by combining the penalty terms introduced above. Note that the penalty terms introduced in Sec. 3.3.2 are quite different from those introduced in Sec. 3.3.3. No cell in B + can be assigned to more than one parent in b. To enforce this constraint, define the symmetric m × m binary matrix [Q j,k ] by (i) Q j,k = 1 if j = k and the two pairs pch j , pch k have one cell in common, (ii) Q j,k = 0 if j = k and the two pairs pch j , pch k have no cell in common, (iii) Q j,j = 0 for all j.
The quadratic penalty z → z, Qz is non-negative for z ∈ CONF , and must be zero if sub(z) = trueCH . Introduce six positive weight parameters to be selected further on λ j , j ∈ {lin, gap, dev, rat, rank, Q}. Define the vector V ∈ R m as a weighted linear combination of the penalty vectors LIN , GAP , DEV , RAT , RANK For any configuration z ∈ CONF , the BM energy E(z) is defined by the quadratic function We already know that the unknown set trueCH of true children pairs must have cardinal DIV = N + − N . So we seek a configuration z * ∈ CONF minimizing the energy E(z) under the rigid constraint card{sub(z)} = DIV . Let ONE ∈ R m be the vector with all its coordinates equal to 1. The constraint on z can be reformulated as ONE , z = DIV . We want the unknown trueCH to be close to the solution z * of the constrained minimization problem To force this minimization to yield a reliable estimate of trueCH , we calibrate the six weights λ j , j ∈ {lin, gap, dev, rat, rank, Q} by the algorithm in Sec. 3.2. Denote CONF 1 the set of all z ∈ CONF such that ONE , z = DIV .
To minimize E(z) under the constraint z ∈ CONF 1 , fix a slowly decreasing temperature scheme Temp(t) as in Sec. 3.4. We need to force the BM stochastic configurations Z(t) to remain in CONF 1 . Then, for large time step t, the Z(t) will converge in probability to a configuration z * ∈ CONF 1 approximately minimizing E(z) under the constraint z ∈ CONF 1 . Start with any Z(0) ∈ CONF 1 . Assume that for 0 ≤ s ≤ t, one has already dynamically generated BM configurations Z(s) ∈ CONF 1 . Then, randomly select two sites j, k such that  A priori reduction of m significantly reduces the computing times, and can be implemented by trimming away the pairs pch j ∈ PCH for which the penalties LIN j , GAP j , DEV j , RAT j , and RANK j are all larger than predetermined empirical thresholds. We performed a study on 100 successive (synthetic) images. We show scatter plots for the most informative penalty terms in Fig. 3. These plots allow us to determine adequate thresholds for the penalty terms. We observed that for the synthetic and real data we considered the trimming of DEV , GAP , and RANK reduced the percentage of invalid children pairs by 95%, therefore drastically reducing the combinatorial complexity of the problem. The quadratic energy function E(z) is the sum of clique energies J K (z) involving only cliques of cardinality 1 and 2. For any clique K = {j} of cardinality 1, with 1 ≤ j ≤ m, one has J K (z) = V j z j . For any clique K = {j, k} of cardinality 2, with 1 ≤ j < k ≤ m, one has J K (z) = 2Q j,k z j z k . A key computational step when generating Z(t + 1) is to evaluate the energy change ∆E when one flips the binary values Z j (t) = 1 and Z k (t) = 0 by the new value (1 − z i ) for a fixed single site i. This step is quite fast since it uses only the numbers V j , V k , and q(j), Z(t) , q(k), Z(t) , where q(i) is the i th row of the matrix Q.
3.6.2. Children Pairing: Implementation on Synthetic Videos. We have implemented our children pairing algorithms on synthetic image sequences having 100 to 500 image frames with 1 minutes interframe (benchmark set BENCH1; see Sec. 2.1). The cell motion bound w/2 per interframe was defined by w = 20 pixels. The parameter τ that defines the sets P CH of plausible children pairs (see (5)) was set at τ = 45 pixels.
The known true cell registrations indicated that in our typical BENCH1 image sequence, the successive sets P CH had average cardinals of 120, while the number of true children pairs per Table 2. Accuracies of parent-children pairing algorithm. We applied our parent-children pairing algorithm to three long synthetic image sequences BENCH1 (500 frames), BENCH2 (300 frames), and BENCH3 (300 frames), with interframe intervals of 1, 2, 3 minutes, respectively. The table summarizes the resulting pcpaccuracies. Note that pcp-accuracies are practically always at 100%. For BENCH2 pcp-accuracies are 100% for 298 frames out of 300, and for the remaining two frames, accuracies were still high at 93% and 96%. For BENCH3 the average pcp-accuracy for the 3 minute interframe is 99%.
sequence pcp-accuracy frames BENCH1 acc = 100% 500 out of 500 BENCH2 acc = 100% 298 out of 300 BENCH2 99% ≥ acc ≥ 93% 2 out of 300 BENCH3 acc = 100% 271 out of 300 BENCH3 99% ≥ acc ≥ 95% 17 out of 300 BENCH3 94% ≥ acc ≥ 90% 12 out of 300 PCH roughly ranged from 2 to 6 with a median of 4. The size of the reduced configuration space CON F 1 per image frame thus ranged from 10 4 to 120 6 /6! = 4.2 · 10 9 with a median of 9 · 10 6 Our weights estimation technique introduced in Sec  1 , b 2 ) ∈ SHL is obtained by direct comparison to the known ground truth registration J → J + . For each frame J, we define the pcp-accuracy of our Parent-Children Pairing algorithm as the ratio VAL/DIV . We have tested our parent-children matching algorithm on three long synthetic image sequences BENCH1 (500 frames), BENCH2 (300 frames), and BENCH3 (300 frames), with respective interframes of 1, 2, and 3 minutes. For each frame J k , we computed the pcp-accuracy between J k and J k+1 .
We report the accuracies of our parent-children pairing algorithms in Tab. 2. For BENCH1, all 500 pcp-accuracies reached 100%. For BENCH2, pcp-accuracies reach 100% for 298 frames out of 300, and for the remaining two frames, accuracies were still high at 93% and 96%. For BENCH3, where interframe duration was longest (3 minutes), the 300 pcp-accuracies decreased slightly but still averaged 99%, and never fell below 90%.  For each (b, b 1 , b 2 ) ∈ SHL, eliminate from B the parent cell, b, and eliminate from B + the two children cells b 1 , b 2 . We are left with two residual sets, resB ⊂ B and resB + ⊂ B + , having the same cardinality, N − DIV = N + − 2DIV . Assuming that our set SHC of short lineages is correctly determined, the cells b ∈ redB should not divide in the interframe J → J + , and hence have a single (still unknown) registration f (b) ∈ redB + . Thus, the still unknown part of the registration f is a bijection from redB to redB + .
Let divB = B − redB and divB + = B + − redB + . For each b ∈ divB, the cell b divides into the unique pair of cells, (b 1 , b 2 ) ∈ divB + × divB + , such that (b, b 1 , b 2 ) ∈ SHL. Hence, we can set f (b) = (b 1 , b 2 ) for all b ∈ divB . Thus, the remaining problem to solve is to compute the bijective registration f : redB → redB + . We have reduced the registration discovery to a new problem, where no cell divisions occur in the interframe duration. In what follows, we present our algorithm to solve this registration problem.

Automatic Cell Registration after Reduction to Cases with No Cell Division.
As indicated above, we can explicitly reduce the generic cell tracking problem to a problem where there is no cell division. We consider images J, J + with associated cell sets B, B + such that N = card(B) = card(B + ). Hence, there are no cell divisions in the interframe J → J + and the map f of this reduced problem is (in principle) a bijection f : B → B + with card(B) = card(B + ). In Fig. 4 we show two typical successive images we use for testing with no cell division generated by the simulation software [106,107] (see Sec. 2.1).
3.8.1. The Set MAP of Many-to-One Cell Registrations. We have reduced the registration search to a situation where during the interframe J → J + , no cell has divided, no cell has disappeared, and no cell has suddenly emerged in B + without originating from B. The unknown registration f : B → B + should then in principle be injective and onto. However, for computational efficiency, we will temporarily relax the bijectivity constraint on f . We will seek f in the set MAP of all many-to-one mappings f : B → B + such that for each b ∈ B, the cell f (b) is in the target window W (b) ⊂ B + (see Sec. 2.3).

Registration Cost Functional.
To design a cost functional cost(f ), which should be roughly minimized when f ∈ MAP is very close to the true registration from B to B + , we linearly combine penalties match(f ), over(f ), stab(f ), flip(f ) weighted by unknown positive weights λ match , λ over , λ stab , λ flip , to write, for all registrations f ∈ MAP , (11) cost(f ) = λ match match(f ) + λ over over(f ) + λ stab stab(f ) + λ flip flip(f ).
We specify the individual terms that appear in (11) below. Ideally, the minimizer of cost(f ) over all f ∈ MAP is close to the unknown true registration mapping f : B → B + . To enforce a good approximation of this situation, we first estimate efficient positive weights by applying our calibration algorithm (see Sec. 3.2). The actual minimization of cost(f ) over all f ∈ MAP is then implemented by a BM described in Sec. 3.9.
Cell Matching Likelihood: match(f ). Here, we extend a pseudo likelihood approach used to estimate parameters in Markov random fields modeling by Gibbs distributions (see [52]). Recall that g.rate is the known average cell growth rate. For any cells b ∈ B, b + ∈ B + , the geometric quality of the matching b → b + relies on three main characteristics: (i) motion c(b + ) − c(b) of the cell center c(b), (ii) angle between the long axes A(b) and A(b + ), (iii) cell length ratio is the geometric angle between the straight lines carrying A(b) and A(b + ).
Fix b ∈ B, and let b run through the whole target window W (b). The finite set of values thus reached by the kinetic penalties kin(b, b ) has two smallest values kin 1 (b), kin 2 (b). Define list.kin = b∈B {kin 1 (b), kin 2 (b)}, which is a list of 2N "low " kinetic penalty values. Repeat this procedure for the penalties dis(b, b ) and rot(b, b ) to similarly define a list.dis of 2N "low " distortion penalty values, and a list.rot of 2N "low " rotation penalty values.
The three sets list.kin, list.dis, list.rot can be viewed as three random samples of size 2N , respectively, generated by three unknown probability distributions P kin , P dis , P rot . We approximate these three probabilities by their empirical cumulative distribution functions CDF kin , CDF dis , CDF rot , which can be readily computed. We now use the right tails of these three CDFs to compute separate probabilistic evaluations of how likely the matching of cell b ∈ B with cell b + ∈ W (b) is. For any fixed mapping f ∈ MAP , and any b ∈ B, set b + = f (b). Compute the three penalties vkin = kin(b, b + ), vdis = dis(b, b + ), vrot = rot(b, b + ), and define three associated "likelihoods" for the matching High values of the penalties vkin, vdis, vrot thus will yield three small likelihoods for the matching b → b + = f (b). With this, we can define a "joint likelihood " 0 ≤ LIK(b, b + ) ≤ 1 evaluating how likely is the matching b → b + = f (b): Note that higher values of LIK(b, b + ) correspond to a better geometric quality for the matching of b with b + = f (b). To avoid vanishingly small likelihoods, whenever LIK(b, b + ) < 10 −6 , we replace , as should be expected for bona fide cells registrations. But for this mapping f , we have z3 above z2 above z1, whereas for the original cells we had b 1 above b 2 above b3. Our cost function penalizes flips of this nature.
it by 10 −6 . Then, for any mapping f ∈ MAP , we define its likelihood lik(f ) by the finite product The product of these N likelihoods is typically very small, since N = card(B) can be large. So, we evaluate the geometric matching quality match(f ) of the mapping f via the averaged log-likelihood of f , namely, .
Good registrations f ∈ MAP should yield small values for the criterion match(f ).
Overlap: over(f ). We expect bona fide cell registrations f ∈ MAP to be bijections. Consequently, we want to penalize mappings f which are many-to-one. We say that two distinct cells . The total number of overlapping pairs (b, b ) for f defines the overlap penalty: Neighbor Stability: stab(f ). Let B = {b 1 , . . . , b N }. Denote G i the set of all neighbors for cell b i in B (i.e., b j ∼ b i ⇐⇒ b j ∈ G i ; see Sec. 2.3). For bona fide registrations f ∈ MAP , and for most pairs of neighbors b i ∼ b j in B, we expect f (b i ) and f (b j ) to remain neighbors in B + . Consequently, we penalize the lack of "neighbors stability" for f by For any registration f ∈ MAP , define the flip penalty for f by where G(b) is the neighborhood of cell b in B. In Fig. 5 we illustrate an example of an unwanted cell flip.

BM Minimization of Registration Cost Function.
In what follows, we define the optimization problem for the registration of cells from one frame to another (i.e., cell tracking), as well as associated methodology and parameter estimates. To any configuration z = {z 1 , . . . , z N } ∈ CONF , we associate a unique cell registration f ∈ MAP defined by f (b j ) = z j for all j, denoted by f = map(z). This determines a bijection z → f = map(z) from CONF onto MAP . The inverse of map : CONF → MAP will be called range : MAP → CONF , and is defined by z = range(f ), when z j = f (b j ) for all j.

BM Energy Function E(z)
. We now define the energy function E(z) ≥ 0 of our BM for all z ∈ CONF . Denote E * = minimize z∈CONF E(z). Since f → z = range(f ) is a bijection from MAP to CONF , we must have Our goal is to minimize cost(f ), and we know that BM simulations should roughly minimize E(z) over all z ∈ CONF . So, we define the BM energy function E(z) by forcing (13) cost(f ) = E(range(f )) for any registration mapping f ∈ MAP , which-due to the preceding subsection-is equivalent to (14) E(z) = cost(map(z)) for all configurations z ∈ CONF . The next subsection will explicitly express the energy E(z) in terms of cliques of neurons. Due to (13) and (14) we have For large time t, the BM stochastic configuration Z(t) tends with high probability to concentrate on configurations z ∈ CONF , which roughly minimize E(z). The random registration F t = map(Z(t)) will belong to MAP and verify Z(t) = range(F t ), so that E(Z(t)) = E(range(F t )) = cost(F t )). Consequently, for large t-with high probability-the random mapping F t = map(Z(t)) will have a value of the cost functional cost(F t ) close to minimize f ∈MAP cost(f ). Cliques in CL 1 . For each clique K = {i} in CL 1 , and each z ∈ CONF , define its energy J match,K (z) = J match,K (z i ) by where LIK is given by (12). Set J match,K ≡ 0 for K in CL 2 ∪ CL 3 . For all z ∈ CONF , define the energy E match (z) by 3.9.4. Test Set of 100 Synthetic Image Pairs. As shown above, the minimization of cost(f ) over all registrations f ∈ MAP is equivalent to seeking BM configurations z ∈ CONF with minimal energy E(z). We have implemented this minimization of E(z) by the long term asynchronous dynamics of the BM just defined. This algorithm was designed for the registration of image pairs exhibiting no cell division, and was, therefore, implemented after the automatic reduction of the generic registration problem, as indicated earlier. 3.9.5. Implementation of BM minimization for cost(f ). The numbers clq 1 , clq 2 , clq 3 of cliques in CL 1 , CL 2 , CL 3 have the following rough ranges 80 ≤ clq1 ≤ 100, 160 ≤ clq 2 ≤ 250, and 450 ≤ clq 3 ≤ 600. For k = 1, 2, 3, denote val (k) the numbers of non-zero values for J K (z) when z runs through CONF and K runs through all cliques of cardinality k. One easily checks the rough upper bounds val (1) < 4 000; val (2) < 200, 000; val (3) < 300, 000. Hence, to automatically register B to B + , one could pre-compute and store all the possible values of J K (z) for all cliques K ∈ CL 1 ∪ CL 2 ∪ CL 3 and all the configurations z ∈ CONF . This accelerates the key computing steps of the asynchronous BM dynamics, namely, for the evaluation of energy change ∆E = E(z ) − E(z), when configurations z and z differs at only one site j ∈ S. Indeed, the single site modification z j → z j affects only the energy values J K (z) for the very small number r(j) of cliques K, which contain the site j. In our benchmark sets of synthetic images, one had r(j) < 24 for all j ∈ S. Hence, the computation of ∆E was fast since it requires retrieving at most 24 pairs of pre-computed J K (z), J K (z ), and evaluating the 24 differences J K (z ) − J K (z). Another practical acceleration step is to replace the ubiquitous computations of probabilities p(t) = exp(−D/Temp(t)) by simply testing the value −D/Temp(t) against 100 precomputed logarithmic thresholds. In our implementation of ABM dynamics, we used virtual temperature schemes such as Temp(t) = 50 · ρ t with 0.995 ≤ ρ ≤ 0.999. The BM simulation was stopped when the stochastic energy E(Z(t)) had remained roughly stable during the last N steps. Since all target windows W (j) had cardinality smaller 40, the initial configuration Z(0) = x was computed via where the likelihoods LIK were defined by (12).
3.9.6. Weight Calibration. For the pair of successive synthetic images J, J + displayed in Fig. 4, we have N = card(B) = card(B + ) = 513 cells. The ground truth registration f is known by construction; we used it to apply the weight calibration described in Sec. 3.2. We set the metaparameter γ to 10 10 and obtained the vector of weights (17) Λ * = [λ * match , λ * over , λ * stab , λ * flip ] = [110,300,300,290]. These weights are kept fixed for all the 100 pairs of images taken from the set BENCH6. The determined weights are used in the cost function cost(f ) defined above. This correctly parametrized the BM energy function E(z). We then simulated the BM stochastic dynamics to minimize the BM energy E(Z(t)).
3.9.7. BM Simulations. We launched 100 simulations of the asynchronous BM dynamics, one for each pair of successive images in our test set of 100 images taken from BENCH6. For each such pair, the ground truth mapping f : B → B + was known by construction and the stochastic minimization of the BM energy generated an estimated cells registration f : B → B + . For each pair B, B + Table 3. Registration accuracy for synthetic image sequence BENCH 100 . We consider 100 pairs of consecutive synthetic images taken from the benchmark dataset BENCH6. Automatic registration was implemented by BM minimization of the cost function cost(f ), which was parametrized by the vector of optimized weights Λ * in (17). The average registration accuracy was 99%.  Each image contains about 100 to 150 cells. Consequently, the runtime for the algorithm is approximately 20 seconds per cell for our prototype implementation. We note that this is only a rough estimate. The runtime depends on several factors, such as the number of cells in an image; the number of mother and daughter cells (i.e., how many cells divide); the size of the neighborhood of each individual cell (window size); the weights used in the cost function (which affects the number of epochs); etc. We note that the temperature scheme had not been optimized yet, so that these computing times are upper bounds. Earlier SBM studies [12,14] indicate that the same energy minimizations on GPUs could provide a computational speedup by a factor ranging between 30 and 50. We report registration accuracies in Tab. 3. For each pair of images in the considered set of 100 images, the accuracy of automatic registration was larger than 94.5%. The overall average registration accuracy was quite high at 99%.

Results
In this section, we report results for the registration for cell dynamics involving growth, motion, and cell divisions. of one minute. We still impose the mild constraint that no cell is lost between two successive images. The main difference with the earlier benchmark of 100 images from BENCH6 is that cells are allowed to freely divide during interframes, as well as to grow and to move. For the full implementation on 100 pairs of successive images, we first execute the parent-children pairing, and remove the identified parent-children triplets; we can then apply our cell registration algorithmic on the reduced sets cells. Our image sequence contained 760 true parent-children triplets, which we automatically identified with an accuracy of 100%. As outlined earlier, we removed all these identified cell triplets and then applied our tracking algorithm. This left us with a total of 12,631 cells (spread over 100 frames). Full automatic registration was then implemented with an accuracy higher than 99.5%.

4.2.
Tests of Cell Registration Algorithms on Laboratory Image Sequences. To test our cell tracking algorithm on pairs of consecutive images extracted from recorded image sequences of bacterial colonies (real data), we had to automatically delineate all individual cells in each image. Representative frames of this data are shown in Fig. 1. We describe these data in more detail in Sec. 2.2. We will only briefly outline the overall segmentation approach to not distract from our main contribution-the cell tracking algorithm. We use the watershed algorithm [19] (also used, e.g., in [54]) to segment each frame into individual image segments containing one single cell each. Consequently, these regions represent over segmentations of the individual cells; we only know that each region will contain a bacteria cell b. To segment individual cells, an additional step is necessary. We then apply adhoc nonlinear filters to remove minor segmentation artifacts. In a second step, we then identified the contour of each single cell b by applying the Mumford-Shah algorithm [73] within the image segment containing a cell b. Since this procedure is quite time consuming for large images, we have implemented it to produce a training set of delineated individual cells to train a CNN for image segmentation. After automatic training, this CNN substantially reduces the runtime of the cell segmentation/delineation procedure. We show the resulting segmentations in Fig. 6. We provide additional information regarding our approach for the segmentation of individual bacteria cells in the appendix (see Sec. D).
After each cell has been identified (i.e., segmented out) in each pair J, J + of successive images, we transform J, J + into binary images, where cells appear in white on a black background. For each resulting pair B, B + of successive sets of cells, we apply the parent-children pairing algorithm outlined in Sec. 3.3 to identify all the short lineages. For the two successive images in COL1, the Table 4. Cost function weights for parent-children pairing in the COL1 images displayed in Fig. 6.
Weights λ cen λ siz λ ang λ gap λ dev λ rat λ rank λ over  Figure 8. Cell tracking results for the short image sequence COL2 in Fig. 6. The interframe duration for COL2 is six minutes. COL2 involves four successive images J(t i ), i = 0, 1, 2, 3. In our figure, each one of the three rows displays the automatic cell registration results between images J(t i ) and J(t i+1 ) for i = 0, 1, 2. We report the accuracies of parent-children pairing and of the registration in Tab. 5. Left column: Results for parent-children pairing. Each parent-children triplet is identified by the same color for each parent cell an its two children. Middle column: Display of the automatically computed registration after removing the parent-children triplets already identified in order to generate two reduced sets redB and redB + of cells. Again, the same color is used for each pair of automatically registered cells. The white cells in redB + are cells which could not be registered to some cell in redB . Right column: To differentiate between errors induced during automatic identification of and errors generated by automatic registration between redB and redB + , we manually removed all "true" parent-children triplets and then applied our registration algorithm to this "cleaned " (reduced) cell sets redB * and redB * + .
discovered short lineages are shown in Fig. 7 (left pair of images). Here, color designates the cell triplet algorithmically identified: parent cell in image J and its two children in image J + . We then remove each identified "parent" from B and its two children from B + . This yields the reduced cell sets redB and redB + . We can then apply our tracking algorithm (see 3.7) dedicated to situations where cells do not divide during the interframe. For image sequences of live cell colonies we had to re-calibrate most of our weight parameters. The weight parameters used for these image sequences are summarized in Tab. 4.
The BM temperature scheme was Temp(t) = 2000 (0.995) t , with the number of epochs capped at 5000. We illustrate our COL1 automatic registration results in Fig. 7 (right pair of images). Here, if cell b ∈ redB has been automatically registered onto cell b + ∈ redB + , b, b + share the same color. The cells colored in white in redB + are cells which the registration algorithm did not succeed in matching to some cell in redB . These errors can essentially be attributed to errors in the parent-children pairing step. By visual inspection we have determined that there are 14 true Table 5. Cell tracking accuracy for the short image sequence COL2 in Fig. 6 with an interframe of six minutes. We report the ratio of correctly predicted cell matches over the total number of true cell matches and the associated percentages. The accuracy results quantify four distinct percentages of correct detections (i) for parent cells in image J, (ii) for children cells in image J + , (iii) for parent-children triplets, and (iv) for registered pairs of cells (b, b + ) ∈ redB × redB + .
task accuracy parent-children triplets in the successive images of COL1. Our parent-children pairing algorithm did correctly identify 11 of these 14 triplets.To check further the performance of our registration algorithm on live images, we also report automatic registration results for "manually prepared " true versions of redB and redB + , obtained by removing "manually" the true parent-children triplets determined by visual inspection. For the short image sequence COL2, results are displayed in Fig. 8. The display setup is the same: The left column shows the results of automatic parent-children pairing. The middle column illustrates the computed registration after automatic removal of the computer identified parent-children triplets. The third column displays the computed registration after removing "manually" the true parent-children triplets determined by visual inspection. Note that the overall matching accuracy can be improved if we reduce errors in the parent-children pairing. We report quantitative accuracies in Tab. 5. For parent-children pairing, accuracy ranges between 70% and 78%. For pure registration after correct parent-children pairing, accuracy ranges between 90% and 100%.

Conclusions and Future Work
We have developed a methodology for automatic cell tracking in recordings of dense bacterial colonies growing in a mono-layer. We have also validated our approach using synthetic data from agent based simulations, as well as experimental recordings of E. coli colonies growing in microfluidic traps. Our next goal is to streamline our implementation for systematic cell registration on experimentally acquired recordings of such cell colonies, to enable automated quantitative analysis and modeling of cell population dynamics and lineages.
There are a number of challenges for our cell tracking algorithm: Inherent imaging artifacts such as noise or intensity drifts, cells overlaps, similarity of cell shape characteristics across the population, tight packing of cells, somewhat large interframe times, cell growth combined with cell motion and cell divisions, represent just a few of these challenges. Overall, the cell tracking problem has combinatorial complexity, and for large frames is beyond the concrete patience of human experts. We tackle these challenges by developing a two-stage algorithm that first identifies parent-children triplets and subsequently computes cell registration from one frame to the next, after reducing the two original cell sets by automatic removal of the identified parent-children triplets. Our algorithms specify innovative cost functions dedicated to these registration challenges. These cost functions have combinatorial complexity. To discover good registrations we minimize these cost functions numerically by intensive stochastic simulations of specifically structured BMs. We have validated the potential of our approach by reporting promising results obtained on long synthetic image sequences of simulated cell colonies (which naturally provide a ground truth for cell registration from one frame to the next). We have also successfully tested our algorithms on experimental recordings of live bacterial colonies.
The choice of adequate cost functions to drive each major cost optimization step in our multi-step cell tracking algorithms is essential to obtaining good tracking. Selecting the proper formulation had a strong impact on actual tracking accuracy. Our cost functions are fundamentally nonlinear, which entails additional complications. We introduced a set of meta-parameters for each cost function, and proposed an original learning algorithm to automatically identify good ranges for these meta-parameters.
Our BMs are focused on stochastic minimization of dedicated cost functions. An interesting feature of BMs we will explore in future work is the simplicity of their natural massive parallelization for fast stochastic minimization [10]. This allows us to mitigate the slow convergence typically observed for Gibbs samplers on discrete state spaces with high cardinality. Parallelized BMs implement a form of massively parallel simulated annealing. Sequential simulated annealing has been explored by physicists [71,50,86,25] seeking to minimize spin-glasses energies. For these cliquebased energies, reaching global minima requires unfeasible CPU times, and much faster parallel simulated annealing yields only good local minima, via a sophisticated but still greedy stochastic search. Parallel stochasticity favors ending in rather stable local minima, which in turn enforces low sensitivity to small changes in energy parameters. Robustness to small changes in the coefficients of our cost functions is a desirable feature, since our algorithmic calibration of cost coefficients focuses on computing good ranges for these meta-parameters. We do not aim to seek global minima, generally a very elusive search, because computing speed and scalability are important features in our problem. Recall the established results of Huber [41] showing that optimal estimators of the mean for a Gaussian distribution lose efficiency very quickly when the Gaussian data are slightly perturbed.
In future work, we will further improve the stability and accuracy of our cell registration algorithms by exploring natural modifications of our cost functions. In the present work, we have not yet explicitly considered the case of cells vanishing between successive frames. This is a critical issue that can occur due to cells exiting or entering the field of view as well as due to errors in cell segmentation. The problem is somewhat controlled and/or mitigated in our experimental setup, where we expect cells to enter or vanish close to a precisely positioned trap edge and/or near frame boundaries. Since we intend to track lineages, each frame-to-frame error of this type may be problematic, and it will be instrumental for our future work to address these issues.
Linking parents to children involves an optimization distinct from the final optimization of frame-to-frame registrations. This did reduce computing time without reducing the quality for our benchmark results. However, in future work one could attempt to iterate this sequence of two optimizations in order to reach a better minimum.
We note that our algorithm does work for experimental setups in which the frame rate of the video recordings is not fixed. This will require an adaptive parameter selection that depends on the frame rate. This can be implemented based on a trivial rescaling procedure. However, note that for larger interframe times, more errors will impact tracking results. Indeed, large interframe durations intensify fluctuations in key parameters of cell dynamics, and increase the range of cell displacement, imposing searches in larger cell neighborhoods for cell pairing, as well as increased combinatorial complexity.
We have considered synthetic data to evaluate the performance of our method. One clear practical issue is that some of the parameters of our tracking algorithms may change when applied to laboratory image sequences acquired from colonies of different cells, with various image acquisition setups. One can design a computational framework to automatically fit the parameters of the simulation model to the imaging data acquired on specific live cell colonies, using specific camera hardware and setup. In future work, we will attempt to implement this type of fitting for our simulation model, before launching intensive model simulations to calibrate the parameters of our new tracking algorithms. We have not yet removed physical scales in the implementation of our tracking algorithm. Implementing such a non-dimensionalization will allow us to reduce the sensitivity of our methodology with respect to new datasets.
Identification of full lineages is an interesting concrete goal for cell tracking. Evaluating the accuracy of lineage identification on real cell colonies is quite challenging since it requires inheritable biological tagging of cells. This is probably feasible for populations mixing two or three cell types, but not for individualized tagging in populations of moderate size. However even partial tagging of sub-populations would provide some control on lineage identification accuracies.
Acknowledgements. This work was partly supported by the National Science Foundation There are two main options to implement the Markov chain dynamics Z(t) → Z(t + 1) (see [9]). (v) For all j = M , the new state Z j (t + 1) of neuron U j remains equal to its current state U j (t).
A.2. Synchronous BM Dynamics. Fix a synchrony parameter 0 < α < 1, usually around 50%. At each time t, all neurons U j synchronously, but independently compute their own random binary tag tag j (t), equal to 1 with probability α, and to 0 with probability (1 − α). Let SYN (t) be the set of all neurons. All the neurons U j such that tag j (t) = 1 then synchronously and independently compute their new random states Z j (t + 1) ∈ W (j) by applying the updating procedure given above. And for all j such that tag j (t) = 0, the new state Z j (t + 1) of U j remains equal to Z j (t).
A.3. Comparing Asynchronous and Synchronous BM Dynamics. As t becomes large, and for temperatures Temp(t) slowly decreasing to 0, both BM dynamics generate with high probability configurations Z(t) which provide deep local minima E(Z(t)) of the BM energy function. The asynchronous dynamics can be fairly slow. But the synchronous dynamics is much faster since it emulates efficient forms of parallelel simulated annealing (see [10,82]) and is directly implementable on GPUs.
Our tracking module is a collection of python functions and has been released to the public at https://github.com/scopagroup/BacTrak. We refer to [106,107] for a detailed description of this mathematical model and its implementation. The code for generating the synthetic data has been released at https://github.com/jwinkle/eQ. We note that detailed installation instructions for the software can be found on this page. The parameters for this agent-based simulation software are as follows: Cells were modeled as 2D spherocylinders of constant, 1 µm width. The computational framework takes into account mechanical constraints that can impact cell growth and influence other aspects of cell behavior. The growth rate of the cells is exponential and is controlled by the doubling time. The time until cells double is set to 20 min (default setting; resulting in a growth rate of g.rate = 1.05). The cells have a length of approximately 2µm after division and 4µm right before division (minimum division length of 4 µm; subject to some random perturbation). In our data set of simulated videos, there is no "trap wall" (as opposed to the simulations carried out in [106,107]). The "trap" encompassing all cells on a given frame has a size of 30 µm × 30 µm subdivided into 400 × 400 pixels of size 0.075 µm × 0.075 µm. The size of the resulting binary image used in our tracking algorithm is 600 × 600 pixels. (We add a boundary of 100 pixels on each side). Bacteria are moving, growing, dividing within the trap. However, at this stage of our study, we consider only video segments where no cell disappears and where cells do not enter the trap from outside so that the trap is a confined environment. Cells move only due to soft shocks interactions with other neighboring cells. The time interval between any two successive image frames ranges from one minute to six minutes (see Tab. 1). All other simulation parameters remain unchanged; i.e., we use the default parameters specified in the simulation software.

Appendix D. Cell Segmentation
In the next couple of section we outline the framework we have developed to segment individual cells from real world laboratory imaging data. In a first step, we consider traditional segmentation algorithms-a watershed algorithm [33,19,101] in combination with a variational contour based model-to generate a sufficiently large dataset to train a neuronal network. The actual segmentations on real data can subsequently be carried out efficiently using segmentation predictions generated by the trained neuronal network. Note that the proposed segmentation algorithm is only included for completeness. We do not view this as a major contribution of the present work. D.1. Watershed Algorithm. We consider a watershed algorithm based on immersion that compares high intensity values to local intensity minima for cell segmentation [33,19,101].
We consider Matlab's implementation of the watershed algorithm in the present work. This version of the watershed algorithm is unseeded and yields n regions R = {R 1 , R 2 , . . . , R n }. To identify these regions, we perform a statistical analysis of each image histogram to compute adaptive rough thresholds for interiors and exterior of cells. This leads to watershed results which identify each cell by a segment slightly larger than the cell itself. The very small percentage of oversegmented cells is automatically detected by cell length and width computations through PCA analysis of each cell shape viewed as a cloud of planar points. Since our segments are slightly too wide we reduce each segment to the exact outer cell contour by applying a Mumford-Shah algorithm to each segment computed by the watershed algorithm. In an ideal case, after applying the watershed algorithm, each individual bacteria cell b i , i = 1, . . . , n, will be located in a single region R i ⊂ R 2 . Table 6. Statistics of some quantities of interest related to the intensity of boundary segments and regions. These quantities allow us to define heuristics to identify erroneous segmentations computed by the watershed algorithm. We state the characteristic and report the minimum, maximum 5% quantile, mean, and standard deviation for the reported quantities of interest.

D.2. Segmentation Errors: Correction
Steps. We define the boundary segment B i,j as a nonempty intersection of two region's boundaries, i.e., B i,j = ∂S i ∩ ∂S j . Moreover, we denote the area of a region R i as area(R i ). We know that the interior of a bacteria cell b i has a lower intensity than the exterior region of a cell. More precisely, the interior of a cell tends to have intensity values of zero whereas the exterior of a cell (i.e., the background) tends to have an intensity that is close or equal to one. For this reason, we define a function for the intensity of the boundary. To remove outliers, we consider the average intensity value of the pixels located along a boundary segment. We denote this mean intensity value along a boundary B i,j by mint(B i,j ) and the average intensity of a region R i by mint(R i ). One difficulty is that we cannot assume that the intensity of the pixels on the interior of each cell corresponds to the same value (i.e., there exist intensity and contrast drifts depending on location). We hypothesize that if mint(B i,j ) of a boundary segment is close to the average intensity of the regions on both sides of the boundary segment B i,j this boundary segment does not separate two bacteria cells; it is erroneous. Conversely, if the difference between the mean intensity along a boundary segment and the mean intensity of the interior regions it separates is high, we consider that the boundary segment represents a good segmentation (i.e., represents a segment that does separate two cells). To quantify this notion, we define the height of a boundary segment as H i,j = mint(B i,j ) − (mint(R i ) + mint(R j ))/2. In Tab. 6, we report some statistics associated with the quantities of interest introduced above. There are several key observations we can draw from this table which confirm our qualitative (i.e., visual) assessment of the segmentation results. Most notably, we can observe that there seem to exist outliers in terms of cell size. Moreover, we can observe that in some cases we obtain a height of the boundary segment that is negative, and by that nonsensical. These observations allow us to develop some heuristic rules to remove erroneous segmentations.
We introduce the following post-processing steps: i) We connect small regions to their neighbors (i.e., regions that are too small in area to realistically contain any cells). We select the threshold for the area to be 65. This threshold is selected in accordance to the scale of the image and the expected size of bacteria cells observed in the image data. We merge each small region with one of its neighboring regions by removing the segment that separates the two. To select an appropriate region for merging, we choose the region that gave the lowest height H i,j from all available candidate regions that share the same boundary segment. ii) We remove all boundary segments B ij with a height H ij that is below the 5% quantile of all heights. iii) We remove all incomplete regions from our segmentation. We define a region as incomplete, if the region or the associated boundary segments touch an edge of the image. This step is necessary since we cannot guarantee that the regions close to the boundary contain an entire cell or only parts of a cell. Consequently, we decided to remove them to prevent any issues with our post-analysis. D.3. Cell Boundary Detection. The next step is to identify the boundaries of individual cells contained within a subregion defined by the watershed algorithm. To identify the boundaries of the cells (and by that segment the individual cells) we use the Mumford-Shah algorithm [73]. Notice that we can execute the Mumford-Shah algorithm for each region R i separately making this an embarrassingly parallelizable problem. Denote the cell in each R i region by b i . We divide each of these regions into three different zones. The first zone is the interior of the cell b i denoted by in(b i ). The second zone is exterior of the cell (i.e., the background) contained in the region and denoted by out(b i ). The third zone is the boundary of the cell b i , denoted by ∂b i . The Mumford-Shah algorithm represents a variational approach that allows us to segment cartoon like images.
Mathematically speaking, we model information contained in each region R i as piecewise-smooth functions. In our model, the associated regions we seek to identify are given by the zones defined above-the interior and the exterior of the cell b i . Let u int (b i ) denote the mean intensity for the interior of the cell b i and u ext (b i ) denote the mean intensity for the exterior of the cell b i . With this definition, we obtain the cost functional where the first two terms measure the discrepancy between the piecewise smooth function u ext and u int and the image intensities u and the third term is a penalty that measures the length of the boundary of a particular cell b i with parameter ν > 0. Notice, that our formulation slightly deviates from the traditional definition of the Mumford-Shah cost functional; we drop the penalty for the smoothness of the function u. The minimizer of the cost function cost MS defined above provides the sought after segmentation: the boundary, interior, and exterior of a cell. We have implemented the minimization of the cost function formula for each cell separately. D.4. Convolutional Neural Networks (CNNs). Next, we introduce our actual method for cell segmentation that can be efficiently applied to a large dataset (as opposed to the prototype method described above to generate the underlying training data). The biggest issue with the methodology outlined above is that our prototype implementation is computationally costly. While we envision that an improved implementation as well as the use of parallel computing can significantly reduce the time to solution, we decided to not further pursue a reduction in runtime but extend our methodology by taking advantage of existing machine learning algorithms. Replacing the approach outlined above by CNNs allowed us to reduce the runtime by factor of 60 to less than 3 minutes, without any significant loss in accuracy.
Training and Testing Data. In the absence of any ground truth data set for the classification of rod-shape bacteria cells from movies of cell populations, we consider the output of the Mumford-Shah algorithm introduced above as ground truth classification for training and testing our machine learning methodology. Above, we introduced three different zones: The interior in(b i ), the exterior ext(b i ), and the boundary ∂b i of a cell b i . We reduce these three regions to two zones-the interior and exterior of a cell b i . We assign pixels that belong to int(b i ) the label 0 and pixels that belong to ext(b i ) and ∂b i the label of 1. For an image of size 200 × 200 we obtain 40,000 binary labels. We limit the training of the CNN to a subregion of size 200 × 200 in the center of each preprocessed image to avoid issues associated with mislabeled training data of cells located at the boundary of our data. We consider X as the set of features and Y as the set of labels. We want to assign to each pixel a label of either 0 or 1. For pixel p, we define X p to be a 7 × 7 square window with center p located in the original image. The corresponding label Y p is denoted by C(p), which corresponds to the class of the pixel p in the binarized image.  CNN Algorithm. The considered CNN algorithm consists of two parts, i) the convolutional auto-encoder and ii) a fully connected multiLayer perceptron (MLP). The input for the autoencoder is a window of 7 × 7 pixels. In the first layer of the encoder, we have a 5 × 5 × 4 convolution layer Conv1 with 3 × 3 kernel. We feed Conv1 to a max-pooling layer MPool2 with one stride and pooling window 2 × 2. The output of MPool2 is the input of a 3 × 3 × 8 convolution layer Conv3. For decoding, we have almost the same structure in reverse order: We feed Conv3 to a 5 × 5 × 4 deconvolution with 3 × 3 kernel. Subsequently, we feed the output of this layer to a 7 × 7 × 1 deconvolution with 3 × 3 kernel. The decoder's output is a window of 7 × 7 pixels. We compare this output with the input window (since it is an auto-encoder, features and labels are the same) by using the mean square error as a cost function. We train the auto-encoder for all training sets using a mini-batch gradient descent. When the training is finished, we freeze the weights for Conv1 and Conv3.
After training the auto-encoder and freezing the weights, we feed X as the input to Conv1 and get the output of Conv3 denoted byX. In the next step, we train an MLP with featuresX and labels Y . We flattenX, which is a 3 × 3 × 8 matrix to a vector of size 72 × 1, called FCL4. FCL4 is fully connected to the hidden layer HID5 with 10 nodes. We use ReLu as a nonlinear function for HID5. We connect HID5 to the output layer OUT6, which possess two nodes for the two classes 0 and 1. We use a softmax function to find two probabilistic outputs p 0 and p 1 = 1 − p 0 for related classes. We use maximum-entropy as a cost function. We train the MLP for training set of (X, Y ) with mini-batch gradient descent.
We have trained the model with two images of size 200 × 200 pixels; the training set is 80,000 7 × 7 images. We train the model for 100 epochs. The accuracy of the model for the image is 93%. The confusion matrix is shown in Tab. 7. Based on this confusion matrix we can observe that the proposed methodology can predict the pixels located in the interior of a cell quite well. However, we can also observe that there is a slightly lower accuracy for the pixels outside the cells. This can be probably explained by the fact that the data sets are tightly packed with cells so that we have available more observations of foreground pixels (interior of cells) than pixels that belong to the background.