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

Abstract

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%.

1. Introduction

Technology advances have led to increasing magnitudes of data generation with increasing levels of precision [1,2]. However, data generation presently far outpaces data analysis and drives the requirement for analyzing such large-scale data sets with automated tools [3,4,5]. 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 [6,7,8].

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 [9,10]. 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 [6,11,12]. 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 [13,14,15], 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.

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 Figure 1. We describe the design and validation of algorithms for tracking individual cells in sequences of such images [11,16,17]. 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 focuses on tracking E. coli in a time series of images. A comparison of different cell-tracking algorithms can be found in [19,20]. 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 [21,22]. In (biomedical) image processing, motion tracking is often referred to as “image registration” [19,23,24,25,26] or “optical flow” [27,28,29,30]. 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” [31,32,33,34,35,36]. 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 Figure 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 [37,38,39] or active contour models [40,41,42,43,44]. 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” [45,46,47,48,49,50,51,52,53,54,55], 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 [37,49,54,56,57,58]. 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 [59] for an overview of cell segmentation approaches. Deep learning strategies have been widely used for this task [5,50,54,58,60,61,62,63,64,65]. We consider a framework based on convolutional neural networks (CNNs). Others have also used CNNs for cell segmentation [62,64,66,67,68]. 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 Section 1.2). However, the interested reader is referred to Appendix 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,48,51,54]. 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 [5,45,49,52]. 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 [20,69] (see, e.g., [37,70]). (Cell tracking challenge: http://celltrackingchallenge.net (accessed on 15 December 2021).) Works that consider these data are for example [37,54,57,58,62,67]. 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 Section 2.1).

Standard graph matching algorithms (see, e.g., [71]) 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 [72,73]. Graph-based matching strategies for cell-tracking that are somewhat related to our approach are described in [70,74,75,76,77]. 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 [75,76]. 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.

Refs. [71,78,79,80] implement multi-target tracking in videos by stochastic models based on random finite set densities and variants thereof. The fit to the data are 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 [37,49,54,56,57,58]. Contrary to observations we have made in our data, these approaches rely (with the exception of [49]) on the fact that the tracking problem is inherent to the segmentation problem (“tracking-by-detection methods” [54]; see also [5]). That is, the key assumption made by many of these algorithms is that cells belonging to the same lineage overlap across frames (see also [47]). In this case, cell-overlap can serve as a good proxy for cell-tracking [54]. We note that in our data we cannot guarantee that the frame rate is sufficiently high for this assumption to hold. Refs. [56,57,67] 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 [12,81].

The works that are most similar to ours are [45,49,52]. 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, Refs. [45,49] employ a local greedy-search, whereas we consider stochastic neural network dynamics for optimization. Ref. [52] constructs score matrices within a score based neighborhood tracking method; an integer programming method is used to generate frame-to-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 [47,54,82].

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, Ref. [45] reports a tracking accuracy of up to 97% for data that are similar to ours, while Ref. [74] reports 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 [75], with an accuracy at the order of 90% (ranging from about 87% to 92%, respectively). Accuracies reported in [77] range from about 92% to 97%, respectively. This work also includes a comparison to one of their earlier approaches [76] with an accuracy of up to 85% and 89% if the datasets are pre-aligned. We note that the data considered in [74,75,76,77] are quite different from ours. Refs. [37,54,57,58,70] consider the data from the cell tracking challenge [20,69] to evaluate the performance of their methods. As in the previously mentioned work, these data are again quite different from ours. To evaluate the performance of the methodology, the so-called acyclic oriented graph matching measure [83] is considered. We refer to the webpage of the cell tracking challenge for details on the evaluation metrics (see http://celltrackingchallenge.net/evaluation-methodology, accessed on 15 December 2021). Based on these, the reported tracking scores are between 0.873 and 0.902 [37], 0.901 and 1.00 [70], 0.950 and 0.987 [54], 0.788 and 0.982 [58] and 0.765 and 0.915 [57] 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 parent–children 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 [84,85,86]. (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. Therefore, we first validated the accuracy of our cell tracking algorithms on recordings of simulated cell colonies, generated by the dedicated cell colony simulation software [12,81]. This provided us with ground truth frame-by-frame registration for cell lineages, enabling us to validate our methodology.

1.3. Outline

In Section 2, we describe the synthetic image sequence (see Section 2.1) and experimental data (see Section 2.2) of cell colonies considered as benchmarks for our cell tracking algorithms. In Section 2.3, we describe key cell characteristics considered in our tracking methodology to define metrics that enter our cost functionals. Our tracking approach is developed in greater detail in Section 3. We define valid cell registration mappings between successive image frames in Section 3.1. We outline how to automatically calibrate the weights of our various penalty terms in Section 3.2. Our algorithms for pairing parent cells with their children and for cell-to-cell registration are developed in Section 3.3, Section 3.4, Section 3.5, Section 3.6, Section 3.7, Section 3.8, Section 3.9. We present our main validation results on long image sequences (time series of images) in Section 4 and conclude with Section 5.

2. Datasets

Below, we introduce the datasets used to evaluate the performance of the proposed methodology. The synthetic data are described in Section 2.1. The experimental data (real imaging data) are described in Section 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 [12,81] 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, accessed on 15 December 2021) The simulated cell colony dynamics are driven by an agent based model [12,81], 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.

The cells are modeled as 2D spherocylinders of constant 1 $\mathsf{\mu }$m width. Each cell grew exponentially in length with a doubling time of 20 min. To prevent division synchronization across the population when a mother cell of length $L_{\mathrm{div}}$ divides, the two daughter cells are assigned random birth lengths $L_0\left(b_1\right)=L_1=\delta L_{\mathrm{div}}$ and $L_0\left(b_2\right)=:L_2=(1-\delta )L_{\mathrm{div}}$, where $\delta >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_{\mathrm{div}}$ divides into two cells $b_1$ and $b_2$, their lengths $L_1$, $L_2$ satisfy $L_1+L_2=L_{\mathrm{div}}$ and $L_i/L_{\mathrm{div}}$, $i=1,2$, is a random number. The cells have a length of approximately 2 $\mathsf{\mu }$m after division and 4 $\mathsf{\mu }$m right before division. We refer to [81] 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 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 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 min, respectively. Notice that there is no explicit noise on the growth rate. However, due to the crowding of cells, the growth rate will vary from cell-to-cell. The generated binary images are of size $600\times 600$ pixels. We summarize these benchmarks in

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 min. We refer to the text for details about how these data have been generated.

Label	Interframe Duration	Number of Frames
BENCH1	1 min	500
BENCH2	2 min	300
BENCH3	3 min	300
BENCH6	6 min	100

. The associated image sequences involve between 100 up to 500 frames, respectively. In Figure 2, we display an example of two simulated consecutive frames separated by 1 min. 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 $\mathsf{\mu }$m in diameter and on average 3 $\mathsf{\mu }$m in length, and they divide about every 30 min. The original images exported from the microscope are 0.11 $\mathsf{\mu }$m/pixel. The microscopy experimental data were obtained using JS006 [87] (BW25113 $\mathsf{\Delta }$araC $\mathsf{\Delta }$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 $\mathsf{\mu }$g/mL ampicillin for 18 h. These cultures were diluted in the morning into 1/1000 into 50 mL fresh LB with 100 $\mathsf{\mu }$g/mL ampicillin and grown for 3 h 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 min 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 [88]. In the microfluidic device, the cells were provided with continuous fresh LB with 100 $\mathsf{\mu }$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 Figure 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\left(b\right)$ is the centroid of all pixels belonging to cell b. The long axis $A\left(b\right)$ of cell b is computed by principal component analysis (PCA). The endpoints $e\left(b\right)$ and $h\left(b\right)$ of cell b are the first and last cell pixels closest to $A\left(b\right)$; see Figure 2 (right) for a schematic illustration.

Cell Neighbors. For each image frame J, denote $B=B\left(J\right)$ as the set of fully visible cells in J, and by $N=N\left(J\right)=card\left(B\right)$ the number of these cells. Let V be the set of all cell centers $c\left(b\right)$ with $b\in B$. Denote $\mathit{delV}$ the Delaunay triangulation [89] of the finite planar set V with N vertices. We say that two cells $b_1$, $b_2$ in B are neighbors if they verify the following three conditions: (i) $(b_1,b_2)$ are connected by the edge $\mathit{edg}$ of one triangle in $\mathit{delV}$. (ii) The edge $\mathit{edg}$ does not intersect any other cell in B. (iii) Their centers verify $\parallel c\left(b_1\right)-c\left(b_2\right)\parallel \le \rho$, where $\rho >0$ is a user defined parameter.

For the synthetic images of size $600\times 600$ that we considered (see Section 2.1), we take $\rho =80$ pixels. We write $b_1\sim 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\left(J\right)$, $B_{+}=B\left(J_{+}\right)$ as the associated sets of cells. Superpose temporarily the images J and $J_{+}$ so that they then have the same center pixel. Any cell $b\in B$, which does not divide in the interframe $J\to J_{+}$, becomes a cell $b_{+}$ in image $J_{+}$. The “motion vector” of cell b from frame J to $J_{+}$ is then defined by $v\left(b\right)=c\left(b_{+}\right)-c\left(b\right)$. If the cell b does divide between J and $J_{+}$, denote $b_{\mathrm{div}}$ as the last position reached by cell b at the time of cell division, and define similarly the motion $v\left(b\right)=c\left(b_{\mathrm{div}}\right)-c\left(b\right)$. In our experimental recordings of real bacterial colonies with interframe duration 6 min, there is a fixed number$w>0$ such that $\parallel v\left(b\right)\parallel \le w/2$ for all cells $b\in B\left(J\right)$ 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 meta-parameter 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\left(b\right)\subset J_{+}$ be a square window of width w, with the same center as cell b. The target window$W\left(b\right)$ is the set of all cells in $B_{+}$ having their centers in $U\left(b\right)$. Since $\parallel v\left(b\right)\parallel \le w/2$, the cell $b_{+}$ must belong to the target window $W\left(b\right)\subset B_{+}$.

3. 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\left(B\right)$, and $N_{+}=card\left(B_{+}\right)$. Our goal is to track each cell from J to $J_{+}$. For each cell $b\in B$, there exist three possible evolutions between J and $J_{+}$:

Case 1: 

Cell $b\in B$ did not divide in the interframe $J\to J_{+}$, and has become a cell $f\left(b\right)\in B_{+}$; that is, $f\left(b\right)$ has grown and moved during the interframe time interval.

Case 2: 

Cell $b\in B$ divided between J and $J_{+}$, and generated two children cells $b_1,b_2\in B_{+}$; we then denote $f\left(b\right)=(b_1,b_2)\in B_{+}\times B_{+}$.

Case 3: 

Cell $b\in B$ disappeared in the interframe $J\to J_{+}$, so that $f\left(b\right)$ is not defined.

To simplify our exposition, we ignore Case 3. We discuss Case 3 in greater detail in the conclusions in Section 5. Consequently, a valid (true) registration mapping f will take values in the set $\left\{B_{+}\right\}\cup \{B_{+}\times B_{+}\}$.

3.2. Calibration of Cost Function Weights

With the notation we introduced, fix any two finite sets A, $A_{+}$. Let $G:=\{g:A\to A_{+}\}$ be the set of all mappings g: $A\to A_{+}$. Fix m penalty functions ${pen}_k\left(g\right)\ge 0$, $k=1,\dots ,m$. Let $g^{*}\in G$ be the ground truth mapping we want to discover through minimization in g of some given cost function $COST\left(g\right)$ defined by the linear combination of the penalty functions ${pen}_k\left(g\right)$, the contributions of which are controlled by the cost function weights ${\lambda }_k>0$. In this section, we present a generic weight calibration algorithm, extending a technique introduced and applied in [90,91] 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\in 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 are 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\left(g\right)$ over $g\in G$. Let $PEN\left(g\right)=[{pen}_1\left(g\right),\dots ,{pen}_m\left(g\right)]$ be the vector of m penalties for any mapping $g\in G$. Let $\mathsf{\Lambda }=[{\lambda }_1,\dots ,{\lambda }_m]$ be the weight vectors. Then, $COST\left(g\right)=\langle \mathsf{\Lambda },PEN\left(g\right)\rangle$. Replacing g by another mapping $h\ne g$ induces the penalty changes $\mathsf{\Delta }{PEN}_{g,h}=PEN\left(h\right)-PEN\left(g\right)$ and the cost change $\mathsf{\Delta }COST(g,h)=\langle \mathsf{\Lambda },\mathsf{\Delta }{PEN}_{g,h}\rangle$. Now, fix any known ground truth mapping $f\in G$. We want f to be a minimizer of COST, so we should have $\mathsf{\Delta }COST(f,f^{\prime })\ge 0$ for all modifications $f\to f^{\prime }\in G$.

For each $a\in A$, select an arbitrary $s\left(a\right)\in W\left(a\right)$ (where $W\left(a\right)$ is the target window for cell a; see Section 2.3), to define a new mapping $f^{\prime }=f_a^{\prime }$ from A to $A_{+}$ by $f_a^{\prime }\left(a\right)=s\left(a\right)$, and $f_a^{\prime }\left(x\right)\equiv f\left(x\right)$ for all $x\ne a$. Since f is a minimizer of COST, this single point modification $f\to f_a^{\prime }$ must generate the following cost increase

$$\langle \mathsf{\Lambda },\mathsf{\Delta }PEN(f,f_a^{\prime })\rangle =\mathsf{\Delta }COST(f,f_a^{\prime })\ge 0.$$

Denote $V_a\in {\mathbb{R}}^m$ the vector $V_a=\mathsf{\Delta }PEN(f,f_a^{\prime })$. Then, the positive vector $\mathsf{\Lambda }\in {\mathbb{R}}^m$, $\mathsf{\Lambda }⪰0$, should verify the set of linear constraints $\langle \mathsf{\Lambda },V_a\rangle \ge 0$ for all $a\in A$. There may be too many such linear constraints. Consequently, we relax these constraints by introducing a vector $y=\left[y\left(a\right)\right]\in {\mathbb{R}}^{card\left(A\right)}$, $y⪰0$, of slack variables $y\left(a\right)\ge 0$ indexed by all the $a\in 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 $\mathsf{\Lambda }$ and the slack variables vector y to verify the system of linear constraints:

$$\begin{array}{cccc}\hfill \langle \mathsf{\Lambda },V_a\rangle +y\left(a\right)& =0\hfill & \hfill & \mathrm{for}\mathrm{all}a\in A\hfill \\ \hfill \mathsf{\Lambda }⪰0,y& ⪰0\hfill & \hfill & \\ \hfill \langle \mathsf{\Lambda },Z\rangle & \le 1000\hfill & \hfill \end{array}$$

(1)

where $Z=[1,\dots ,1]\in {\mathbb{R}}^m$. The normalizing constant 1000 can be arbitrarily changed by rescaling. We seek high positive values for $\mathsf{\Delta }COST(f,f_a^{\prime })$ and small $L_1$-norm for the slack variable vector y. Thus, we will seek two vectors $\mathsf{\Lambda }\in {\mathbb{R}}^m$ and $y\in {\mathbb{R}}^{card\left(A\right)}$ solving the following convex-concave minimization problem, where $\gamma >0$ is a user selected (large) meta parameter:

$$\underset{\mathsf{\Lambda },y}{minimize}\gamma {\parallel y\parallel }_{L1}-\sum _{a\in A}{\left[\langle \mathsf{\Lambda },V_a\rangle \right]}^{+}$$

(2)

subject to (1), where we denote ${\left[x\right]}^{+}:=max(x,0)$ for arbitrary x. To numerically solve the constrained minimization problem (2), we use the libraries CVXPY and DCCP (disciplined convex-concave programming) [84,85,86]. DCCP is a package for convex-concave programming designed to solve non-convex problems. (DCCP can be downloaded at https://github.com/cvxgrp/dccp (last accessed on 20 January 2022)) It can handle objective functions and constraints with any known curvature as defined by the rules of disciplined convex programming [92]. We give examples of numerically computed weight vectors $\mathsf{\Lambda }$ below. The computing time was less than 30 s 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 $\mathsf{\Lambda }$ 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. Fix two successive synthetic image frames J, $J_{+}$ with short interframe time equal to 1 minute. Their cell sets B, $B_{+}$ have cardinality N and $N_{+}$, respectively. For our synthetic image sequences, all cells $b\in B$ still exist in $B_{+}$—either as whole cells or after dividing into two children cells, and no new cell enters the field of view during the interframe $J\to J_{+}$. This forces $N_{+}\ge N$, and implies that the number $\mathit{DIV}$ of cell divisions occurring in this interframe verifies $\mathit{DIV}=N_{+}-N$. Each children pair $(b_1,b_2)\in B_{+}\times B_{+}$ is born from a single parent $b\in B$. Thus, the set $\mathit{trueCH}$ of all such true children pairs must then verify

$$card\left(\mathit{trueCH}\right)=\mathit{DIV}=N_{+}-N.$$

(3)

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(\mathit{trueCH})-\mathit{DIV}|$ may be of the order of $n_{in}+n_{out}$. To take this into account, we relax the constraint in (3) as follows:

$$|card(\mathit{trueCH})-\mathit{DIV}|\le \mathit{REL},$$

(4)

where $\mathit{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. However, even in that case, there was a computational advantage to using the slightly relaxed constraint (4) with $\mathit{REL}=1$.

3.3.2. Most Likely Parent Cell for a Given Children Pair

For successive images J, $J_{+}$ with 1 min interframe, define the set $\mathit{PCH}$ of plausible children pairs by

$$\mathit{PCH}=\{(b_1,b_2)\in B_{+}\times B_{+}\mathrm{with}\mathrm{centers}{c}_1,c_2\mathrm{verifying}\parallel c_1-c_2\parallel <\tau \},$$

(5)

where the threshold $\tau >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)\in \mathit{PCH}$ can qualify as a pair of children generated by division of a parent cell $b\in B$, we now quantify the “geometric distortion” between b and $(b_1,b_2)$. Cell division of b into $b_1,b_2\in 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 Section 2.3), and their initial distance to the center c of b is roughly at most $\parallel A\left(b\right)\parallel /4$, where $A\left(b\right)$ is the long axis of cell b. Hence, the centers c, $c_1$, $c_2$ of b, $b_1$, $b_2$ should verify the constraint

$$max\{\parallel c_1-c\parallel ,\parallel c_2-c\parallel \}\le w+\parallel A\parallel /4.$$

(6)

Define the set $\mathit{SHLIN}$ of potential short lineages as the set all triplets $(b,b_1,b_2)$ with $b\in B$, $(b_1,b_2)\in \mathit{PCH}$, verifying the preceding constraint (6). For each potential lineage $(b,b_1,b_2)\in \mathit{SHLIN}$, define three terms penalizing the geometric distortions between a parent $b\in B$ and a pair of children $(b_1,b_2)\in \mathit{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(b,b_1,b_2)=\parallel c-(c_1+c_2)/2\parallel$, ($ii$) size distortion $siz(b,b_1,b_2)=|\parallel A\parallel -(\parallel A_1\parallel +\parallel A_2\parallel \left)\right|$, and ($iii$) angle distortion

$$ang(b,b_1,b_2)=angle(A,A_1)+angle(A,A_2)+angle(A,c_2-c_1).$$

Here, angle denotes “angles between non-oriented straight lines,” with a range from 0 to $\pi /2$. Introduce three positive weights ${\lambda }_{\mathrm{cen}}$, ${\lambda }_{\mathrm{siz}}$, ${\lambda }_{\mathrm{ang}}$ (to be estimated), and for every short lineage $(b,b_1,b_2)\in \mathit{SHLIN}$ define its distortion cost by

$$dist(b,b_1,b_2)={\lambda }_{\mathrm{cen}}cen(b,b_1,b_2)+{\lambda }_{\mathrm{siz}}siz(b,b_1,b_2)+{\lambda }_{\mathrm{ang}}ang(b,b_1,b_2).$$

(7)

For each plausible pair of children $(b_1,b_2)\in \mathit{PCH}$, we will compute the most likely parent cell$b^{*}=parent(b_1,b_2)$ as the cell $b^{*}\in B$ minimizing $dist(b,b_1,b_2)$ in (7) over all $b\in B$, as summarized by the formula

$$b^{*}=parent(b_1,b_2)=\underset{\{b\in B\mid (b,b_1,b_2)\in \mathit{SHLIN}\}}{argmin}dist(b,b_1,b_2).$$

(8)

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 ${\lambda }_j$, $j\in \{\mathrm{cen},\mathrm{siz},\mathrm{ang}\}$ by the algorithm outlined in Section 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\in B$ has very small size due to the constraint (6). Hence, brute force minimization of the functional $dist(b,b_1,b_2)$ in (7) over all $b\in B$ such that $(b,b_1,b_2)\in \mathit{SHLIN}$, is a fast computation for each $(b_1,b_2)$ in $\mathit{PCH}$. The distortion minimizing $b=b^{*}$ yields the most likely parent cell $parent(b_1,b_2)=b^{*}$. The brute force minimization in b of $dist(b,b_1,b_2)$ 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 $\mathit{pch}=(b_1,b_2)$ should belong to $\mathit{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)\in \mathit{PCH}$ should be correctly matchable with their most likely parent cell $b^{*}=parent(b_1,b_2)$ (see (8)). Thus, we define the “lineage” penalty $lin(b_1,b_2)=dist(b^{*},b_1,b_2)$ by

$$lin(b_1,b_2)=\underset{\{b\in b\mid (b,b_1,b_2)\in \mathit{shlin}\}}{argmin}dist(b,b_1,b_2)=dist(parent(b_1,b_2),b_1,b_2).$$

Notice that the computation of $lin(b_1,b_2)$ is quite fast.

“Gap” Penalty. Denote $\mathit{tips}\left(b\right)$ as the set of two endpoints of any cell b. For any pair $\mathit{pch}=(b_1,b_2)\in \mathit{PCH}$, define endpoints $x_1\in \mathit{tips}\left(b_1\right),x_2\in \mathit{tips}\left(b_2\right)$ and the gap penalty $gap(b_1,b_2)$ by

$$gap(b_1,b_2)=\parallel x_1-x_2\parallel =min\{\parallel x-y\parallel \mathrm{for}(x,y)\in \mathit{TIPS}\}$$

(9)

with $\mathit{TIPS}=\mathit{tips}\left(b_1\right)\times \mathit{tips}\left(b_2\right)$.

“Dev” Penalty. For rod-shaped bacteria, a true pair $(b_1,b_2)\in \mathit{PCH}$ of just born children must have a small $gap(b_1,b_2)=\parallel x_1-x_2\parallel$ and roughly aligned cells $b_1$ and $b_2$. For $(b_1,b_2)\in \mathit{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 ${\mathit{str}}_{12}$ be the straight line linking the centers $c_1$, $c_2$ of $b_1$, $b_2$. Let $d_1$, $d_2$ be the distances from $x_1$, $x_2$ to the line ${\mathit{str}}_{12}$. Then, set

$$dev(b_1,b_2)=\frac{d_1+d_2}{\parallel c_2-c_1\parallel }.$$

“Ratio” Penalty. True children pairs must have nearly equal lengths. Thus, for $(b_1,b_2)\in \mathit{PCH}$ with lengths $L_1$, $L_2$, we define the length ratio penalty by

$$ratio(b_1,b_2)=|(L_1/L_2)+(L_2/L_1)-2|.$$

“Rank” Penalty. Let $L_{\min }$ be the minimum cell length over all cells in $B_{+}$. In $B_{+}$, children pairs $(b_1,b_2)$ just born during interframe $J\to J_{+}$ must have lengths $L_1$, $L_2$ close to $L_{min}$. Thus, for $(b_1,b_2)\in \mathit{PCH}$, we define the rank penalty by

$$rank(b_1,b_2)=|(L_1/L_{\min })-1|+|(L_2/L_{\min })-1|.$$

Given two successive images J, $J_{+}$, we seek the set $X=\mathit{trueCH}$ of true children pairs in $B_{+}\times B_{+}$, which is an unknown subset of $\mathit{PCH}$. In Section 3.5 below, we replace X by its indicator function z and we build a cost function $E\left(z\right)$ which should be nearly minimized when z is close to the indicator of $\mathit{trueCH}$. A key term of $E\left(z\right)$ will be a weighted linear combination of the penalty functions $\{lin,gap,dev,ratio,rank\}$. Since these penalties are different from those introduced in Section 3.3.2, we estimate their weights in the cost function $E\left(z\right)$ by the algorithm outlined in Section 3.2. The minimization of $E\left(z\right)$ will be implemented by simulations of a BM with energy function $E\left(z\right)$. We present these stochastic neural networks in the next section.

3.4. 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 [93,94]). Indeed, these recurrent stochastic neural networks can efficiently emulate some forms of simulated annealing.

Each BM implemented here is a network $\mathit{BM}=\{U_1,\dots ,U_N\}$ of N stochastic neurons $U_j$. In the BM context, the time $t=0,1,2,\dots$ is discretized and represents the number of steps in a Markov chain, where the successive updates $Z\left(t\right)\to Z(t+1)$ of the BM configuration $Z\left(t\right)$ are analogous to the steps of a Gibbs sampler. The configuration$Z\left(t\right)=\{Z_1\left(t\right),\dots ,Z_N\left(t\right)\}$ of the whole network $\mathit{BM}$ at time t is defined by the random states $Z_j\left(t\right)$ of all neuron $U_j$. Each $Z_j\left(t\right)$ belongs to a fixed finite set $W\left(j\right)$. Hence, $Z\left(t\right)$ belongs to the configurations set $\mathit{CONF}=W\left(1\right)\times \cdots \times W\left(N\right)$.

Neuron interactivity is specified by a finite set $\mathit{CLQ}$ of cliques. Each clique K is a subset of $S=\{1,\dots ,N\}$. During configuration updates $Z\left(t\right)\to 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\left(z\right)$ defined for all $z\in \mathit{CONF}$, with $J_K\left(z\right)$ depending only on the $z_j$ such that $j\in K$. The full energy $E\left(z\right)$ of configuration z is then defined by

$$E\left(z\right)={\sum }_{K\in \mathit{CLQ}}{J}_K\left(z\right).$$

The BM stochastic dynamics $Z\left(t\right)\to Z(t+1)$ is driven by the energy function $E\left(z\right)$, and by a fixed decreasing sequence of virtual temperatures $\mathit{Temp}\left(t\right)>0$, tending slowly to 0 as $t\to \infty$. Here, we use standard temperature schemes of the form $\mathit{Temp}\left(t\right)\equiv c{\eta }^t$ with fixed $c>0$ and slow decay rate $0.99<\eta <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 [95,96,97]). The detailed BM dynamics are presented in the appendix (see Appendix A).

When the virtual temperatures $\mathit{Temp}\left(t\right)$ decrease slowly enough to 0, the energy $E\left(Z\right(t\left)\right)$ converges in probability to a local minimum of the BM energy $E\left(z\right)$ over all configurations $z\in \mathit{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), which we used for lab videos of live cell populations. Fix successive images J, $J_{+}$ with a positive number of cell divisions $\mathit{DIV}=N_{+}-N$. Denote $\mathit{PCH}=\{{\mathit{pch}}_1,{\mathit{pch}}_2,\dots ,{\mathit{pch}}_m\}$ the set of m plausible children pairs $(b_1,b_2)$ in $B_{+}$. The penalties lin, gap, dev, ratio, and rank defined above for all pairs $(b_1,b_2)\in \mathit{PCH}$ determine five numerical vectors $\mathit{LIN}$, $\mathit{GAP}$, $\mathit{DEV}$, $\mathit{RAT}$, $\mathit{RANK}$ in ${\mathbb{R}}^m$ with coordinates ${\mathit{LIN}}_j=lin\left({\mathit{pch}}_j\right)$, ${\mathit{GAP}}_j=gap\left({\mathit{pch}}_j\right)$, ${\mathit{DEV}}_j=dev\left({\mathit{pch}}_j\right)$, ${\mathit{RAT}}_j=ratio\left({\mathit{pch}}_j\right)$, ${\mathit{RANK}}_j=rank\left({\mathit{pch}}_j\right)$.

We now define a binary BM constituted by mbinary stochastic neurons $U_j$, $j=1\cdots m$. At time $t=0,1,2,\dots$, each $U_j$ has a random binary valued state$Z_j\left(t\right)=1$ or 0. The random configuration $Z\left(t\right)=[Z_1\left(t\right),\dots ,Z_m\left(t\right)]$ of this BM belongs to the configuration space $\mathit{CONF}={\{0,1\}}^m$ of all binary vectors $z=[z_1,\dots ,z_m]$. Let $\mathit{SUB}$ be the set of all subsets of $\mathit{PCH}$. Each configuration $z\in \mathit{CONF}$ is the indicator function of a subset $\mathit{sub}\left(z\right)$ of $\mathit{PCH}$. We view each $\mathit{sub}\left(z\right)\in \mathit{SUB}$ as a possible estimate for the unknown set $\mathit{trueCH}\subset B_{+}\times B_{+}$ of true children pairs $(b_1,b_2)$. For each potential estimate $\mathit{sub}\left(z\right)$ of $\mathit{trueCH}$, the “lack of quality” of the estimate $\mathit{sub}\left(z\right)$ will be penalized by the energy function$E\left(z\right)\ge 0$ of our binary BM. We now specify the energy $E\left(z\right)$ for all $z\in \mathit{CONF}$ by combining the penalty terms introduced above. Note that the penalty terms introduced in Section 3.3.2 are quite different from those introduced in Section 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\times m$ binary matrix $\left[Q_{j,k}\right]$ by (i) $Q_{j,k}=1$ if $j\ne k$ and the two pairs ${\mathit{pch}}_j$, ${\mathit{pch}}_k$ have one cell in common, (ii) $Q_{j,k}=0$ if $j\ne k$ and the two pairs ${\mathit{pch}}_j$, ${\mathit{pch}}_k$ have no cell in common, (iii) $Q_{j,j}=0$ for all j.

The quadratic penalty $z\mapsto \langle z,Qz\rangle$ is non-negative for $z\in \mathit{CONF}$, and must be zero if $\mathit{sub}\left(z\right)=\mathit{trueCH}$. Introduce six positive weight parameters to be selected further on ${\lambda }_j$, $j\in \{\mathrm{lin},\mathrm{gap},\mathrm{dev},\mathrm{rat},\mathrm{rank},Q\}$. Define the vector $V\in {\mathbb{R}}^m$ as a weighted linear combination of the penalty vectors $\mathit{LIN}$, $\mathit{GAP}$, $\mathit{DEV}$, $\mathit{RAT}$, $\mathit{RANK}$

$$V={\lambda }_{\mathrm{lin}}\mathit{LIN}+{\lambda }_{\mathrm{gap}}\mathit{GAP}+{\lambda }_{\mathrm{dev}}\mathit{DEV}+{\lambda }_{\mathrm{rat}}\mathit{RAT}+{\lambda }_{\mathrm{rank}}\mathit{RANK}.$$

For any configuration $z\in \mathit{CONF}$, the BM energy $E\left(z\right)$ is defined by the quadratic function

$$E\left(z\right)=\langle V,z\rangle +{\lambda }_Q\langle z,Qz\rangle .$$

(10)

We already know that the unknown set $\mathit{trueCH}$ of true children pairs must have cardinal $\mathit{DIV}=N_{+}-N$. Thus, we seek a configuration $z^{*}\in \mathit{CONF}$ minimizing the energy $E\left(z\right)$ under the rigid constraint $card\left\{\mathit{sub}\right(z\left)\right\}=\mathit{DIV}$. Let $\mathit{ONE}\in {\mathbb{R}}^m$ be the vector with all its coordinates equal to 1. The constraint on z can be reformulated as $\langle \mathit{ONE},z\rangle =\mathit{DIV}$. We want the unknown $\mathit{trueCH}$ to be close to the solution $z^{*}$ of the constrained minimization problem

$$z^{*}=\underset{z\in \mathit{CONF}}{argmin}E\left(z\right)\mathrm{subject}\mathrm{to}\langle \mathit{ONE},z\rangle =\mathit{DIV}.$$

To force this minimization to yield a reliable estimate of $\mathit{trueCH}$, we calibrate the six weights

$${\lambda }_j,j\in \{\mathrm{lin},\mathrm{gap},\mathrm{dev},\mathrm{rat},\mathrm{rank},Q\}$$

by the algorithm in Section 3.2. Denote ${\mathit{CONF}}_1$ the set of all $z\in \mathit{CONF}$ such that $\langle \mathit{ONE},z\rangle =\mathit{DIV}$. To minimize $E\left(z\right)$ under the constraint $z\in {\mathit{CONF}}_1$, fix a slowly decreasing temperature scheme $\mathit{Temp}\left(t\right)$ as in Section 3.4. We need to force the BM stochastic configurations $Z\left(t\right)$ to remain in ${\mathit{CONF}}_1$. Then, for large time step t, the $Z\left(t\right)$ will converge in probability to a configuration $z^{*}\in {\mathit{CONF}}_1$ approximately minimizing $E\left(z\right)$ under the constraint $z\in {\mathit{CONF}}_1$.

Start with any $Z\left(0\right)\in {\mathit{CONF}}_1$. Assume that, for $0\le s\le t$, one has already dynamically generated BM configurations $Z\left(s\right)\in {\mathit{CONF}}_1$. Then, randomly select two sites j, k such that $Z_j\left(t\right)=1$ and $Z_k\left(t\right)=0$. Compute a virtual configuration Y by setting $Y_j=0$, $Y_k=1$, and $Y_i\equiv Z_i$ for all sites i different from j and k. Compute the energy change $\mathsf{\Delta }E=E\left(Y\right)-E\left(Z\right(t\left)\right)$, and the probability $p\left(t\right)=exp(-D/\mathit{Temp}(t\left)\right)$, where $D=max\{0,\mathsf{\Delta }E\}$. Then, randomly select $Z(t+1)=Y$ or $Z(t+1)=Z\left(t\right)$ with respective probabilities $p\left(t\right)$ and $(1-p(t\left)\right)$. Clearly, this forces $Z(t+1)\in {\mathit{CONF}}_1$.

3.6. Performance of Automatic Children Pairing on Synthetic Videos

In the following subsections, we provide experimental results for pairing children and parent cells.

3.6.1. Children Pairing: Fast BM Simulations

For $m=card\left(\mathit{PCH}\right)\le 1000$, one can reduce the computational cost for BM dynamics simulations by pre-computing and storing the $m\times m$ symmetric binary matrix Q, as well as the m-dimensional vectors $\mathit{LIN}$, $\mathit{GAP}$, $\mathit{DEV}$, $\mathit{RAT}$, $\mathit{RANK}$ and their linear combination V. A priori reduction of m significantly reduces the computing times, and can be implemented by trimming away the pairs ${\mathit{pch}}_j\in \mathit{PCH}$ for which the penalties ${\mathit{LIN}}_j$, ${\mathit{GAP}}_j$, ${\mathit{DEV}}_j$, ${\mathit{RAT}}_j$, and ${\mathit{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 Figure 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 $\mathit{DEV}$, $\mathit{GAP}$, and $\mathit{RANK}$ reduced the percentage of invalid children pairs by 95%, therefore drastically reducing the combinatorial complexity of the problem.

The quadratic energy function $E\left(z\right)$ is the sum of clique energies $J_K\left(z\right)$ involving only cliques of cardinality 1 and 2. For any clique $K=\left\{j\right\}$ of cardinality 1, with $1\le j\le m$, one has $J_K\left(z\right)=V_j{z}_j$. For any clique $K=\{j,k\}$ of cardinality 2, with $1\le j<k\le m$, one has $J_K\left(z\right)=2Q_{j,k}{z}_j{z}_k$. A key computational step when generating $Z(t+1)$ is to evaluate the energy change $\mathsf{\Delta }E$ when one flips the binary values $Z_j\left(t\right)=1$ and $Z_k\left(t\right)=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 $\langle q\left(j\right),Z\left(t\right)\rangle$, $\langle q\left(k\right),Z\left(t\right)\rangle$, where $q\left(i\right)$ 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 min interframe (benchmark set BENCH1; see Section 2.1). The cell motion bound $w/2$ per interframe was defined by $w=20$ pixels. The parameter $\tau$ that defines the sets $PCH$ of plausible children pairs (see (5)) was set at $\tau =45$ pixels.

The known true cell registrations indicated that, in our typical BENCH1 image sequence, the successive sets $PCH$ had average cardinalities of 120, while the number of true children pairs per $\mathit{PCH}$ roughly ranged from 2 to 6 with a median of 4. The size of the reduced configuration space $CONF1$ per image frame thus ranged from ${10}^4$ to ${120}^6/6!=4.2\times {10}^9$ with a median of $9\times {10}^6$.

Our weights estimation technique introduced in Section 3.2 yields the weights

$$[{\lambda }_{\mathrm{cen}},{\lambda }_{\mathrm{siz}},{\lambda }_{\mathrm{ang}}]=[0.255,0.05,0.05]$$

and

$$[{\lambda }_{\mathrm{gap}},{\lambda }_{\mathrm{dev}},{\lambda }_{\mathrm{rat}},{\lambda }_{\mathrm{rank}}]=[0.01,1,0.0001,0.05]$$

for the penalties introduced in Section 4. To reduce the computing time for hundreds of BM energy minimizations on the BENCH1 image sequences, we excluded obviously invalid children pairs in each $\mathit{PCH}$ set, by simultaneously thresholding of the penalty terms. The BM temperature scheme was $\mathit{Temp}\left(t\right)=1000{\left(0.995\right)}^t$, with the number of epochs capped at 5000. The average CPU time for BM energy minimization dedicated to optimized children pairing was about 30 seconds per frame. (We provide hardware specifications in Appendix B).

3.6.3. Parent–Children Matching: Accuracy on Synthetic Videos

For each successive image pair J, $J_{+}$, with cells B, $B_{+}$ of cardinality $N<N_{+}$, our parent–children matching algorithm computes a set $\mathit{SHL}$ of short lineages $(b,b_1,b_2)$, where the cell $b\in B$ is expected to be the parent of cells $b_1,b_2\in B_{+}$. Recall that $\mathit{DIV}=N_{+}-N$ provides the number of cell divisions during the interframe $J\to J_{+}$. The number $\mathit{VAL}$ of correctly reconstructed short lineages $(b,b_1,b_2)\in \mathit{SHL}$ is obtained by direct comparison to the known ground truth registration $J\to J_{+}$. For each frame J, we define the pcp-accuracy of our parent–children pairing algorithm as the ratio $\mathit{VAL}/\mathit{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 min. 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

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 min, respectively. The table summarizes the resulting pcp-accuracies. 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 min interframe is 99%.

Sequence	Pcp-Accuracy	Frames
BENCH1	$\mathit{acc}=100\%$	500 out of 500
BENCH2	$\mathit{acc}=100\%$	298 out of 300
BENCH2	$99\%\ge \mathit{acc}\ge 93\%$	2 out of 300
BENCH3	$\mathit{acc}=100\%$	271 out of 300
BENCH3	$99\%\ge \mathit{acc}\ge 95\%$	17 out of 300
BENCH3	$94\%\ge \mathit{acc}\ge 90\%$	12 out of 300

. 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 min), the 300 pcp-accuracies decreased slightly but still averaged 99%, and never fell below 90%.

3.7. Reduction to Registrations with No Cell Division

Fix successive frames $J,J_{+}$ and their cell sets B, $B_{+}$. We seek the unknown registration mapping $f:B\to \{B_{+}\cup (B_{+}\times B_{+})\}$, where $f\left(b\right)\in B_{+}$ iff cell b did not divide during the interframe $J\to J_{+}$ and $f\left(b\right)=(b_1,b_2)\in B_{+}\times B_{+}$ iff cell b divided into $(b_1,b_2)$ during the interframe.

If $card\left(B\right)=N<N_{+}=card\left(B_{+}\right)$, we know that the number of cell divisions during the interframe $J\to J_{+}$ should be $\mathit{DIV}=\mathit{DIV}(B,B_{+})=N_{+}-N>0$. We then apply the parent–children matching algorithm outlined above to compute a set $\mathit{SHL}=\mathit{SHL}(B,B_{+})$ of short lineages $(b,b_1,b_2)$ with $b\in B$, $b_1,b_2\in B_{+}$ and $card\left(\mathit{SHL}\right)=\mathit{DIV}$. For each $(b,b_1,b_2)\in SHL$, the cell b is computed by $b=parent(b_1,b_2)$ as the parent cell of the two children cells $b_1,b_2\in B_{+}$.

For each $(b,b_1,b_2)\in \mathit{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, $\mathit{resB}\subset B$ and ${\mathit{resB}}_{+}\subset B_{+},$ having the same cardinality, $N-\mathit{DIV}=N_{+}-2\mathit{DIV}$. Assuming that our set $\mathit{SHC}$ of short lineages is correctly determined, the cells $b\in \mathit{redB}$ should not divide in the interframe $J\to J_{+}$, and hence have a single (still unknown) registration $f\left(b\right)\in {\mathit{redB}}_{+}$. Thus, the still unknown part of the registration f is a bijection from $\mathit{redB}$ to ${\mathit{redB}}_{+}$.

Let $\mathit{divB}=B-\mathit{redB}$ and ${\mathit{divB}}_{+}=B_{+}-{\mathit{redB}}_{+}$. For each $b\in divB$, the cell b divides into the unique pair of cells, $(b_1,b_2)\in {\mathit{divB}}_{+}\times {\mathit{divB}}_{+},$ such that $(b,b_1,b_2)\in \mathit{SHL}$. Hence, we can set $f\left(b\right)=(b_1,b_2)$ for all $b\in \mathit{divB}$. Thus, the remaining problem to solve is to compute the bijective registration $f:\mathit{redB}\to {\mathit{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.

3.8. 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\left(B\right)=card\left(B_{+}\right)$. Hence, there are no cell divisions in the interframe $J\to J_{+}$ and the map f of this reduced problem is (in principle) a bijection $f:B\to B_{+}$ with $card\left(B\right)=card\left(B_{+}\right)$. In Figure 4, we show two typical successive images we use for testing with no cell division generated by the simulation software [12,81] (see Section 2.1).

3.8.1. The Set $\mathit{MAP}$ of Many-to-One Cell Registrations

We have reduced the registration search to a situation where, during the interframe $J\to 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\to 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 $\mathit{MAP}$ of all many-to-one mappings$f:B\to B_{+}$ such that for each $b\in B$, the cell $f\left(b\right)$ is in the target window $W\left(b\right)\subset B_{+}$ (see Section 2.3).

3.8.2. Registration Cost Functional

To design a cost functional $cost\left(f\right)$, which should be roughly minimized when $f\in \mathit{MAP}$ is very close to the true registration from B to $B_{+}$, we linearly combine penalties $match\left(f\right)$, $over\left(f\right)$, $stab\left(f\right)$, $flip\left(f\right)$ weighted by unknown positive weights ${\lambda }_{\mathrm{match}}$, ${\lambda }_{\mathrm{over}}$, ${\lambda }_{\mathrm{stab}}$, ${\lambda }_{\mathrm{flip}}$, to write, for all registrations $f\in \mathit{MAP}$,

$$cost\left(f\right)={\lambda }_{\mathrm{match}}match\left(f\right)+{\lambda }_{\mathrm{over}}over\left(f\right)+{\lambda }_{\mathrm{stab}}stab\left(f\right)+{\lambda }_{\mathrm{flip}}flip\left(f\right).$$

(11)

We specify the individual terms that appear in (11) below. Ideally, the minimizer of $cost\left(f\right)$ over all $f\in \mathit{MAP}$ is close to the unknown true registration mapping $f:B\to B_{+}$. To enforce a good approximation of this situation, we first estimate efficient positive weights by applying our calibration algorithm (see Section 3.2). The actual minimization of $cost\left(f\right)$ over all $f\in \mathit{MAP}$ is then implemented by a BM described in Section 3.9.

Cell Matching Likelihood: $match\left(f\right)$. Here, we extend a pseudo likelihood approach used to estimate parameters in Markov random fields modeling by Gibbs distributions (see [98]). Recall that $g.\mathit{rate}$ is the known average cell growth rate. For any cells $b\in B$, $b_{+}\in B_{+}$, the geometric quality of the matching $b\mapsto b_{+}$ relies on three main characteristics: (i) motion $c\left(b_{+}\right)-c\left(b\right)$ of the cell center $c\left(b\right)$, ($ii$) angle between the long axes $A\left(b\right)$ and $A\left(b_{+}\right)$, ($iii$) cell length ratio $\parallel A\left(b_{+}\right)\parallel /\parallel A\left(b\right)\parallel$. Thus, for all $b\in B$ and $b_{+}$ in the target window $W\left(b\right)$, define (i) Kinetic energy: $kin(b,b_{+})={\parallel c\left(b\right)-c\left(b_{+}\right)\parallel }^2$. (ii) Distortion of cell length: $dis(b,b_{+})=|log(\parallel A\left(b_{+}\right)\parallel /\parallel A\left(b\right){\parallel )-logg.\mathit{rate}|}^2$. (iii) Rotation angle: $0\le rot(b,b_{+})\le \pi /2$ is the geometric angle between the straight lines carrying $A\left(b\right)$ and $A\left(b_{+}\right)$.

Fix $b\in B$, and let $b^{\prime }$ run through the whole target window $W\left(b\right)$. The finite set of values thus reached by the kinetic penalties $kin(b,b^{\prime })$ has two smallest values ${\mathit{kin}}_1\left(b\right)$, ${\mathit{kin}}_2\left(b\right)$. Define $\mathit{list}.\mathit{kin}={\bigcup }_{b\in B}\{{kin}_1\left(b\right),{kin}_2\left(b\right)\}$, which is a list of $2N$ “low” kinetic penalty values. Repeat this procedure for the penalties $dis(b,b^{\prime })$ and $rot(b,b^{\prime })$ to similarly define a $\mathit{list}.\mathit{dis}$ of $2N$ “low” distortion penalty values, and a $\mathit{list}.\mathit{rot}$ of $2N$ “low” rotation penalty values.

The three sets $\mathit{list}.\mathit{kin}$, $\mathit{list}.\mathit{dis}$, $\mathit{list}.\mathit{rot}$ can be viewed as three random samples of size $2N$, respectively, generated by three unknown probability distributions $P_{\mathrm{kin}}$, $P_{\mathrm{dis}}$, $P_{\mathrm{rot}}$. We approximate these three probabilities by their empirical cumulative distribution functions ${CDF}_{\mathrm{kin}}$, ${CDF}_{\mathrm{dis}}$, ${CDF}_{\mathrm{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\in B$ with cell $b_{+}\in W\left(b\right)$ is. For any fixed mapping $f\in \mathit{MAP}$, and any $b\in B$, set $b_{+}=f\left(b\right)$. Compute the three penalties $\mathit{vkin}=kin(b,b_{+})$, $\mathit{vdis}=dis(b,b_{+})$, $vrot=rot(b,b_{+})$, and define three associated “likelihoods” for the matching $b\to b_{+}=f\left(b\right)$:

$$\begin{array}{cc}\hfill {LIK}_{\mathrm{kin}}(b,b_{+})& =1-{CDF}_{\mathrm{kin}}\left(vkin\right),\hfill \\ \hfill {LIK}_{\mathrm{dis}}(b,b_{+})& =1-{CDF}_{\mathrm{dis}}\left(vdis\right),\hfill \\ \hfill {LIK}_{\mathrm{rot}}(b,b_{+})& =1-{CDF}_{\mathrm{rot}}\left(vrot\right).\hfill \end{array}$$

High values of the penalties $\mathit{vkin}$, $\mathit{vdis}$, $\mathit{vrot}$ thus will yield three small likelihoods for the matching $b\to b_{+}=f\left(b\right)$. With this, we can define a “joint likelihood” $0\le LIK(b,b_{+})\le 1$ evaluating how likely is the matching $b\to b_{+}=f\left(b\right)$:

$$LIK(b,b_{+})=\prod _{j\in \{\mathrm{kin},\mathrm{dis},\mathrm{rot}\}}{LIK}_j(b,b_{+}).$$

(12)

Note that higher values of $LIK(b,b_{+})$ correspond to a better geometric quality for the matching of b with $b_{+}=f\left(b\right)$. To avoid vanishingly small likelihoods, whenever $LIK(b,b_{+})<{10}^{-6}$, we replace it by ${10}^{-6}$. Then, for any mapping $f\in \mathit{MAP}$, we define its likelihood $lik\left(f\right)$ by the finite product

$$lik\left(f\right)=\prod _{b\in B}LIK(b,f\left(b\right)).$$

The product of these N likelihoods is typically very small, since $N=card\left(B\right)$ can be large. Thus, we evaluate the geometric matching quality $match\left(f\right)$ of the mapping f via the averaged log-likelihood of f, namely,

$$match\left(f\right)=-\frac{1}{N}loglik\left(f\right)=-\frac{1}{N}\sum _{b\in B}logLIK(b,f\left(b\right)).$$

Good registrations $f\in \mathit{MAP}$ should yield small values for the criterion $match\left(f\right)$.

Overlap: $over\left(f\right)$. We expect bona fide cell registrations $f\in \mathit{MAP}$ to be bijections. Consequently, we want to penalize mappings f which are many-to-one. We say that two distinct cells $(b,b^{\prime })\in B\times B$ do overlap for the mapping $f\in \mathit{MAP}$ if $f\left(b\right)=f\left(b^{\prime }\right)$. The total number of overlapping pairs $(b,b^{\prime })$ for f defines the overlap penalty:

$$over\left(f\right)=\frac{1}{card\left(B\right)}\sum _{b\in B}\sum _{b^{\prime }\in B}{1}_{f\left(b\right)=f\left(b^{\prime }\right)}.$$

Neighbor Stability: $stab\left(f\right)$. Let $B=\{b_1,\dots ,b_N\}$. Denote $G_i$ as the set of all neighbors for cell $b_i$ in B (i.e., $b_j\sim b_i\iff b_j\in G_i$; see Section 2.3). For bona fide registrations $f\in \mathit{MAP}$, and for most pairs of neighbors $b_i\sim b_j$ in B, we expect $f\left(b_i\right)$ and $f\left(b_j\right)$ to remain neighbors in $B_{+}$. Consequently, we penalize the lack of “neighbors stability” for f by

$$stab\left(f\right)=\sum _i\sum _{j\ne i}\frac{1}{N\left|G_i\right|\left|G_j\right|}{1}_{b_i\sim b_j}{1}_{f\left(b_i\right)\neg \sim f\left(b_j\right)}.$$

Neighbor Flip: $flip\left(f\right)$. Fix any mapping $f\in \mathit{MAP}$, any cell $b\in B$ and any two neighbors $b^{\prime }$, $b^{″}$ of b in B. Let $z=f\left(b\right)$, $z^{\prime }=f\left(b^{\prime }\right)$, $z^{″}=f\left(b^{″}\right)$. Let c, $c^{\prime }$, $c^{″}$ and d, $d^{\prime }$, $d^{″}$ be the centers of cells b, $b^{\prime }$, $b^{″}$ and z, $z^{\prime }$, $z^{″}$. Let $\alpha$ be the oriented angle between $c^{\prime }-c$ and $c^{″}-c$, and let ${\alpha }_f$ be the angle between $d^{\prime }-d$ and $d^{″}-d$, respectively. We say that the mapping f has flipped cells $b^{\prime }$, $b^{″}$ around b, and we set $FLIP(f,b,b^{\prime },b^{″})=1$ if $z^{\prime }$, $z^{″}$ are both neighbors of z, and the two angles $\alpha$, ${\alpha }_f$ have opposite signs. In all other cases, we set $FLIP(f,b,b^{\prime },b^{″})=0$.

For any registration $f\in \mathit{MAP}$, define the flip penalty for f by

$$flip\left(f\right)=\sum _{b\in B}\sum _{b^{\prime }\in B}\sum _{b^{″}\in B}\frac{1}{N{\left|G\left(b\right)\right|}^2}FLIP(f,b,b^{\prime },b^{″}),$$

where $G\left(b\right)$ is the neighborhood of cell b in B. In Figure 5, we illustrate an example of an unwanted cell flip.

3.9. 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.

3.9.1. BM Minimization of $Cost\left(f\right)$ over $f\in \mathit{MAP}$

Let B, $B_{+}$ be two successive sets of cells. As outlined above, we have reduced the problem to one in which we can assume that $N=card\left(B\right)=card\left(B_{+}\right)$, so that there is no cell division during the interframe. Write $B=\{b_1,\dots ,b_N\}$. For short, denote $W\left(j\right)\subset B_{+}$ instead of $W\left(b_j\right)$ the target window of cell $b_j$. We seek to minimize $cost\left(f\right)$ over all registrations $f\in \mathit{MAP}$. Let $\mathit{BM}$ be a BM with sites $S=\{1,\dots ,N\}$ and stochastic neurons $\{U_1,\dots ,U_N\}$. At time t, the random state $Z_j\left(t\right)$ of $U_j$ will be some cell $z_j$ belonging to the target window $W\left(j\right)$ and the random configuration $Z\left(t\right)=\{Z_1\left(t\right),\dots ,Z_N\left(t\right)\}$ of the whole $\mathit{BM}$ belongs to the configurations set $\mathit{CONF}=W\left(1\right)\times \dots \times W\left(N\right)$.

To any configuration $z=\{z_1,\dots ,z_N\}\in \mathit{CONF}$, we associate a unique cell registration $f\in \mathit{MAP}$ defined by $f\left(b_j\right)=z_j$ for all j, denoted by $f=map\left(z\right)$. This determines a bijection $z\mapsto f=map\left(z\right)$ from $\mathit{CONF}$ onto $\mathit{MAP}$. The inverse of $map:\mathit{CONF}\to \mathit{MAP}$ will be called $range:\mathit{MAP}\to \mathit{CONF}$, and is defined by $z=range\left(f\right)$, when $z_j=f\left(b_j\right)$ for all j.

3.9.2. BM Energy Function $E\left(z\right)$

We now define the energy function $E\left(z\right)\ge 0$ of our BM for all $z\in \mathit{CONF}$. Denote $E^{*}={minimize}_{z\in \mathit{CONF}}E\left(z\right)$. Since $f\mapsto z=range\left(f\right)$ is a bijection from $\mathit{MAP}$ to $\mathit{CONF}$, we must have

$$E^{*}=\underset{z\in \mathit{CONF}}{minimize}E\left(z\right)=\underset{f\in \mathit{MAP}}{minimize}E(range\left(f\right)).$$

Our goal is to minimize $cost\left(f\right)$, and we know that BM simulations should roughly minimize $E\left(z\right)$ over all $z\in \mathit{CONF}$. Thus, we define the BM energy function $E\left(z\right)$ by forcing

$$cost\left(f\right)=E(range(f\left)\right)$$

(13)

for any registration mapping $f\in \mathit{MAP}$, which—due to the preceding subsection—is equivalent to

$$E\left(z\right)=cost(map(z\left)\right)$$

(14)

for all configurations $z\in \mathit{CONF}$. The next subsection will explicitly express the energy $E\left(z\right)$ in terms of cliques of neurons. Due to (13) and (14), we have

$$E^{*}=\underset{f\in \mathit{MAP}}{minimize}cost\left(f\right)=\underset{z\in \mathit{CONF}}{minimize}E\left(z\right).$$

For large time t, the BM stochastic configuration $Z\left(t\right)$ tends with high probability to concentrate on configurations $z\in \mathit{CONF}$, which roughly minimize $E\left(z\right)$. The random registration $F^t=map\left(Z\left(t\right)\right)$ will belong to $\mathit{MAP}$ and verify $Z\left(t\right)=range\left(F^t\right)$, so that $E\left(Z\left(t\right)\right)=E(range\left(F^t\right))=cost\left(F^t\right)$. Consequently, for large t—with high probability—the random mapping $F^t=map\left(Z\left(t\right)\right)$ will have a value of the cost functional $cost\left(F^t\right)$ close to ${minimize}_{f\in \mathit{MAP}}cost\left(f\right)$.

3.9.3. Cliques of Interactive Neurons

The BM energy function $E\left(z\right)$ just defined turns out to involve only three sets of small cliques: (i) ${\mathit{CL}}_1$ is the set of all singletons $K=\left\{i\right\}$, with $i=1\dots N$. (ii) ${\mathit{CL}}_2$ is the set of all pairs $K=\{i,j\}$ such that cells $b_i$ and $b_j$ are neighbors in B. (iii) ${\mathit{CL}}_3$ is the set of all triplets $K=\{i,j,k\}$ such that cells $b_j$ and $b_k$ are both neighbors of $b_i$ in B. Denote $\mathit{CLQ}={\mathit{CL}}_1\cup {\mathit{CL}}_2\cup {\mathit{CL}}_3$ as the set of all cliques for our BM.

Cliques in ${\mathit{CL}}_1$. For each clique $K=\left\{i\right\}$ in ${\mathit{CL}}_1$, and each $z\in \mathit{CONF}$, define its energy $J_{\mathrm{match},K}\left(z\right)=J_{\mathrm{match},K}\left(z_i\right)$ by

$$J_{\mathrm{match},K}\left(z\right)=-\frac{1}{N}logLIK(b_i,z_i)\mathrm{for}\mathrm{all}z\in \mathit{ZW},$$

where LIK is given by (12). Set $J_{\mathrm{match},K}\equiv 0$ for K in ${\mathit{CL}}_2\cup {\mathit{CL}}_3$. For all $z\in \mathit{CONF}$, define the energy $E_{\mathrm{match}}\left(z\right)$ by

$$E_{\mathrm{match}}\left(z\right)=\sum _{K\in \mathit{CLQ}}{J}_{\mathrm{match},K}\left(z\right)=\sum _{K\in {\mathit{CL}}_1}{J}_{\mathrm{match},K}\left(z\right),$$

which implies that the registration $f=map\left(z\right)$ verifies $match\left(f\right)=E_{\mathrm{match}}\left(z\right)$.

Cliques in ${\mathit{CL}}_2$. For all $z\in \mathit{CONF}$, all cliques $K=\{i,j\}$ in ${\mathit{CL}}_2$, define the clique energies $J_{\mathrm{over},K}\left(z\right)=J_{\mathrm{over},K}(z_i,z_j)$ and $J_{\mathrm{stab},K}\left(z\right)=J_{\mathrm{stab},K}(z_i,z_j)$ by $J_{\mathrm{over},K}\left(z\right)=1_{z_i=z_j}/N$ and

$$J_{\mathrm{stab},K}\left(z\right)=\frac{1}{N\left|G_i\right|\left|G_j\right|}{1}_{b_j\sim b_i}{1}_{z_j\neg \sim z_i},$$

where $\left|G_i\right|$ and $\left|G_j\right|$ are the numbers of neighbors in B for cells $z_i$ and $z_j$, respectively. Set $J_{\mathrm{over},K}=J_{\mathrm{stab},K}\equiv 0$ for K in ${\mathit{CL}}_1\cup {\mathit{CL}}_3$. Define the two energy functions

$$\begin{array}{cc}\hfill E_{\mathrm{over}}\left(z\right)& =\sum _{K\in \mathit{CLQ}}{J}_{\mathrm{over},K}\left(z\right)=\sum _{K\in {\mathit{CL}}_2}{J}_{\mathrm{over},K}\left(z\right),\hfill \\ \hfill E_{\mathrm{stab}}\left(z\right)& =\sum _{K\in \mathit{CLQ}}{J}_{\mathrm{stab},K}\left(z\right)=\sum _{K\in {\mathit{CL}}_2}{J}_{\mathrm{stab},K}\left(z\right),\hfill \end{array}$$

which implies that $f=map\left(z\right)$ verifies $over\left(f\right)=E_{\mathrm{over}}\left(z\right)$ and $stab\left(f\right)=E_{\mathrm{stab}}\left(z\right)$.

Cliques in ${\mathit{CL}}_3$. For each clique $K=\{i,j,k\}$ in ${\mathit{CL}}_3$, define the clique energy $J_{\mathrm{flip},K}$ by

$$J_{\mathrm{flip},K}\left(z\right)=J_{\mathrm{flip}}^{i,j,k}\left(z\right)=\frac{1}{N{\left|G_i\right|}^2}FLIP(f^{i,j,k},b_i,b_j,b_k),$$

where $f^{i,j,k}$ is any registration mapping $b_i$, $b_j$, $b_k$ onto $z_i$, $z_j$, $z_k$. The indicator FLIP was defined in Section 3.8.2. Set $J_{\mathrm{flip},K}\equiv 0$ for K in ${\mathit{CL}}_1\cup {\mathit{CL}}_2$. Define the energy

$$E_{\mathrm{flip}}\left(z\right)=\sum _{K\in \mathit{CLQ}}{J}_{\mathrm{flip},K}\left(z\right)=\sum _{K\in {\mathit{CL}}_3}{J}_{\mathrm{flip},K}\left(z\right),$$

which implies that $f=F\left(z\right)$ verifies $flip\left(f\right)=E_{\mathrm{flip}}\left(z\right)$.

Finally, define the clique energy $J_K$ for all $K\in \mathit{CLQ}$ by the linear combination

$$J_K={\lambda }_{\mathrm{match}}{J}_{\mathrm{match},K}+{\lambda }_{\mathrm{over}}{J}_{\mathrm{over},K}+{\lambda }_{\mathrm{stab}}{J}_{\mathrm{stab},K}+{\lambda }_{\mathrm{flip}}{J}_{\mathrm{flip},K}.$$

Summing this relation over all $K\in \mathit{CLQ}$ yields

$$\sum _{K\in \mathit{CLQ}}{J}_K={\lambda }_{\mathrm{match}}{E}_{\mathrm{match}}+{\lambda }_{\mathrm{over}}{E}_{\mathrm{over}}+{\lambda }_{\mathrm{stab}}{E}_{\mathrm{stab}}+{\lambda }_{\mathrm{flip}}{E}_{\mathrm{flip}}.$$

(15)

Define then the final BM energy function $z\mapsto E\left(z\right)$ by

$$E\left(z\right)=\sum _{K\in \mathit{CLQ}}{J}_K\left(z\right)\mathrm{for}\mathrm{all}z\mathrm{in}\mathit{CONF}.$$

(16)

For any $z\in \mathit{CONF}$, the associated registration $f=map\left(z\right)$ verifies $match\left(f\right)=E_{\mathrm{match}}\left(z\right),over\left(f\right)=E_{\mathrm{over}}\left(z\right)$, $stab\left(f\right)=E_{\mathrm{stab}}\left(z\right)$, $flip\left(f\right)=E_{\mathrm{flip}}\left(z\right)$. By weighted linear combination of these equalities, and, due to (15), we obtain for all configurations $z\in \mathit{CONF}$, $E\left(z\right)=cost\left(f\right)$ when $f=map\left(z\right)$ or, equivalently, when $z=range\left(f\right)$.

3.9.4. Test Set of 100 Synthetic Image Pairs

As shown above, the minimization of $cost\left(f\right)$ over all registrations $f\in \mathit{MAP}$ is equivalent to seeking BM configurations $z\in \mathit{CONF}$ with minimal energy $E\left(z\right)$. We have implemented this minimization of $E\left(z\right)$ 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. We have tested this specialized registration algorithm on a set of 100 pairs of successive images of simulated cell colonies exhibiting no cell divisions. These 100 image pairs were extracted from the benchmark set BENCH6 of synthetic image sequence described in Section 2.1. The 100 pairs of cell sets B, $B_{+}$ had sizes $N=card\left(B\right)=card\left(B_{+}\right)$ ranging from 80 to 100 cells. For each test pair B, $B_{+}$, each target window $W\left(j\right)$ typically contained 30 to 40 cells. The set $\mathit{CONF}$ of configurations had huge cardinality ranging from ${10}^{130}$ to ${10}^{160}$. However, the average number of neighbors of a cell was around 4 to 5.

3.9.5. Implementation of BM Minimization for $Cost\left(f\right)$

The numbers ${\mathit{clq}}_1$, ${\mathit{clq}}_2$, ${\mathit{clq}}_3$ of cliques in ${\mathit{CL}}_1$, ${\mathit{CL}}_2$, ${\mathit{CL}}_3$ have the following rough ranges $80\le clq1\le 100$, $160\le {\mathit{clq}}_2\le 250$, and $450\le {\mathit{clq}}_3\le 600$. For $k=1,2,3$, denote $\mathit{val}\left(k\right)$ the numbers of non-zero values for $J_K\left(z\right)$ when z runs through $\mathit{CONF}$ and K runs through all cliques of cardinality k. One easily checks the rough upper bounds $\mathit{val}\left(1\right)<4000$; $\mathit{val}\left(2\right)<\mathrm{200,000}$; $\mathit{val}\left(3\right)<\mathrm{300,000}$. Hence, to automatically register B to $B_{+}$, one could pre-compute and store all the possible values of $J_K\left(z\right)$ for all cliques $K\in {\mathit{CL}}_1\cup {\mathit{CL}}_2\cup {\mathit{CL}}_3$ and all the configurations $z\in \mathit{CONF}$. This accelerates the key computing steps of the asynchronous BM dynamics, namely, for the evaluation of energy change $\mathsf{\Delta }E=E\left(z^{\prime }\right)-E\left(z\right)$, when configurations z and $z^{\prime }$ differ at only one site $j\in S$. Indeed, the single site modification $z_j\to z_j^{\prime }$ affects only the energy values $J_K\left(z\right)$ for the very small number $r\left(j\right)$ of cliques K, which contain the site j. In our benchmark sets of synthetic images, one had $r\left(j\right)<24$ for all $j\in S$. Hence, the computation of $\mathsf{\Delta }E$ was fast since it requires retrieving at most 24 pairs of pre-computed $J_K\left(z\right)$, $J_K\left(z^{\prime }\right)$, and evaluating the 24 differences $J_K\left(z^{\prime }\right)-J_K\left(z\right)$. Another practical acceleration step is to replace the ubiquitous computations of probabilities $p\left(t\right)=exp(-D/\mathit{Temp}(t\left)\right)$ by simply testing the value $-D/\mathit{Temp}\left(t\right)$ against 100 precomputed logarithmic thresholds.

In our implementation of ABM dynamics, we used virtual temperature schemes such as $\mathit{Temp}\left(t\right)=50·{\rho }^t$ with $0.995\le \rho \le 0.999$. The BM simulation was stopped when the stochastic energy $E\left(Z\right(t\left)\right)$ had remained roughly stable during the last N steps. Since all target windows $W\left(j\right)$ had cardinality smaller 40, the initial configuration $Z\left(0\right)=x$ was computed via

$$x_j=\underset{y\in W\left(j\right)}{argmax}LIK(b_j,y)\mathrm{for}j=1,\dots ,N,$$

where the likelihoods LIK were defined by (12).

3.9.6. Weight Calibration

For the pair of successive synthetic images J, $J_{+}$ displayed in Figure 4, we have $N=card\left(B\right)=card\left(B_{+}\right)=513$ cells. The ground truth registration f is known by construction; we used it to apply the weight calibration described in Section 3.2. We set the meta-parameter $\gamma$ to ${10}^{10}$ and obtained the vector of weights

$${\mathsf{\Lambda }}^{*}=[{\lambda }_{\mathrm{match}}^{*},{\lambda }_{\mathrm{over}}^{*},{\lambda }_{\mathrm{stab}}^{*},{\lambda }_{\mathrm{flip}}^{*}]=[110,300,300,290].$$

(17)

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\left(f\right)$ defined above. This correctly parametrized the BM energy function $E\left(z\right)$. We then simulated the BM stochastic dynamics to minimize the BM energy $E\left(Z\right(t\left)\right)$.

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\to B_{+}$ was known by construction and the stochastic minimization of the BM energy generated an estimated cells registration $f^{\prime }:B\to B_{+}$. For each pair B, $B_{+}$ in the considered set of 100 images, the accuracy of this automatically computed registration $f^{\prime }$ was evaluated by the percentage of cells $b\in B$ such that $f^{\prime }\left(b\right)=f\left(b\right)$. When $card\left(B\right)=N$, our BM has N stochastic neurons, and the asynchronous BM dynamics proceeds by successive epochs. Each epoch is a sequence of N single site updates of the BM configuration. For each one of our 100 simulations of BM asynchronous dynamics, the number of epochs ranged from 250 to 450.

The average computing time was about eight minutes per epoch, which entailed a computing time ranging from 30 to 50 min for each one of our 100 automatic registrations $f^{\prime }:B\to B_{+}$ reported here. (We specify the hardware used to carry out these computations in Appendix B). Each image contains about 100 to 150 cells. Consequently, the runtime for the algorithm is approximately 20 s 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 [99,100,101,102] 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

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\left(f\right)$, which was parametrized by the vector of optimized weights ${\mathsf{\Lambda }}^{*}$ in (17). The average registration accuracy was 99%.

Registration Accuracy	Number of Frames
$acc=100\%$	55 frames out of 100
$99\%\ge \mathit{acc}>$ 97%	40 frames out of 100
$96\%\ge \mathit{acc}>$ 94.5%	5 frames out of 100

. 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%.

4. Results

In this section, we report results for the registration for cell dynamics involving growth, motion, and cell divisions.

4.1. Tests of Cell Registration Algorithms on Synthetic Data

We now consider more generic long synthetic image sequences of simulated cell colonies, with a small interframe duration 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 these data are shown in Figure 1. We describe these data in more detail in Section 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 [103] (also used, e.g., in [76]) 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 ad hoc 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 [104] 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 Figure 6. We provide additional information regarding our approach for the segmentation of individual bacteria cells in the appendix (see Appendix 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 Section 3.3 to identify all the short lineages. For the two successive images in COL1, the discovered short lineages are shown in Figure 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 $\mathit{redB}$ and ${\mathit{redB}}_{+}$. We can then apply our tracking algorithm (see Section 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

Table 4. Cost function weights for parent–children pairing in the COL1 images displayed in Figure 6.

Weights	${\mathit{\lambda }}_{\mathbf{cen}}$	${\mathit{\lambda }}_{\mathbf{siz}}$	${\mathit{\lambda }}_{\mathbf{ang}}$	${\mathit{\lambda }}_{\mathbf{gap}}$	${\mathit{\lambda }}_{\mathbf{dev}}$	${\mathit{\lambda }}_{\mathbf{rat}}$	${\mathit{\lambda }}_{\mathbf{rank}}$	${\mathit{\lambda }}_{\mathbf{over}}$
Value	3	7	100	$0.8$	4	$0.01$	$0.01$	600

.

The BM temperature scheme was $\mathit{Temp}\left(t\right)=2000{\left(0.995\right)}^t$, with the number of epochs capped at 5000. We illustrate our COL1 automatic registration results in Figure 7 (right pair of images). Here, if cell $b\in \mathit{redB}$ has been automatically registered onto cell $b_{+}\in {\mathit{redB}}_{+}$, b, $b_{+}$ share the same color. The cells colored in white in ${\mathit{redB}}_{+}$ are cells which the registration algorithm did not succeed in matching to some cell in $\mathit{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 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 $\mathit{redB}$ and ${\mathit{redB}}_{+}$, obtained by removing “manually” the true parent–children triplets determined by visual inspection. For the short image sequence COL2, results are displayed in Figure 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 “manually” removing 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 Table 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%.

Table 5. Cell tracking accuracy for the short image sequence COL2 in Figure 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_{+})\in \mathit{redB}\times {\mathit{redB}}_{+}$.

Task	Accuracy
	$\{{\mathbf{t}}_{\mathbf{0}},{\mathbf{t}}_{\mathbf{1}}\}$		$\{{\mathbf{t}}_{\mathbf{1}},{\mathbf{t}}_{\mathbf{2}}\}$		$\{{\mathbf{t}}_{\mathbf{2}},{\mathbf{t}}_{\mathbf{3}}\}$
correctly detected parents	15/19	79%	20/21	95%	7/10	70%
correctly detected children	35/38	92%	32/42	76%	14/20	70%
correct parent–children triplets	15/19	78%	16/21	76%	7/10	70%
correctly registered cell pairs	36/36	100%	44/49	90%	76/80	95%

Table 5. Cell tracking accuracy for the short image sequence COL2 in Figure 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_{+})\in \mathit{redB}\times {\mathit{redB}}_{+}$.

Task	Accuracy
	$\{{\mathbf{t}}_{\mathbf{0}},{\mathbf{t}}_{\mathbf{1}}\}$		$\{{\mathbf{t}}_{\mathbf{1}},{\mathbf{t}}_{\mathbf{2}}\}$		$\{{\mathbf{t}}_{\mathbf{2}},{\mathbf{t}}_{\mathbf{3}}\}$
correctly detected parents	15/19	79%	20/21	95%	7/10	70%
correctly detected children	35/38	92%	32/42	76%	14/20	70%
correct parent–children triplets	15/19	78%	16/21	76%	7/10	70%
correctly registered cell pairs	36/36	100%	44/49	90%	76/80	95%

5. 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, cell 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 for 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 [90]. 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 [105,106,107,108] seeking to minimize spin–glasses energies. For these clique-based 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 [109] 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.