1.1.

Existing Approaches

The existing automated image analysis techniques for cell segmentation and tracking can be classified into two categories: tracking by model evolution and tracking by detection. The former category detects and tracks the cells simultaneously by representing the cells with mathematical models, whereas the latter category first detects cell candidates in all the frames independently and then links the detected candidates into tracks.

In model evolution algorithms, the mathematically represented cell contours are smoothly evolved from one time frame to the next one, using a velocity term defined by the content of the “target frame.”^10–14 In general, model evolution algorithms require a high imaging frequency to allow sufficient cell overlap between frames and accurate cell contour/boundary detection. These approaches are not applicable to our images due to our low imaging frequency (one day) and unclear cell boundaries (tight cell junctions).

In tracking by detection algorithms, cells are first detected in all the frames independently using thresholding of pixel intensities,^15,16 edge detection,¹⁷ image restoration,^18–20 or template matching.²¹ Kanade et al.¹⁸ applied image restoration to segment the cells on a 2-D flat surface from phase-contrast images, as shown in Fig. 1(a). Due to the different culture environments, the morphologies of our cells within angiogenic vessels are different [see Fig. 1(b)]. Moreover, the thickness of the MFDs is $120\mu \mathrm{m}$, whereas the thickness of a cell on a flat surface is around 3 to $20\mu \mathrm{m}$. Due to the artifacts of phase-contrast microscopy, imaging thicker specimens produces severe halo effects that reduce the visibility of our cell boundaries unlike the data by Kanade et al.¹⁸ In addition to that, some cells may be slightly out of focus. Therefore, the image restoration technique is not applicable for our images. A different approach is template matching; a template is compared to patches of a gray-level image by means of a match statistics evaluated for each patch.²² Euclidean distance, mean square difference, and support vector machine (SVM) can work as the template matching criteria for cell detection.^21,23–26 Since we can easily obtain cell templates from end-point phase-contrast images by aligning the latter with their corresponding end-point confocal images, we apply convolutional neural networks (CNN) to detect cell candidates from the time-lapse phase-contrast images.

Fig. 1

Examples of cells (a) on a 2-D surface and (b) within 3-D angiogenic vessels (red arrows point to cells).

Subsequently, association algorithms, such as “nearest-neighborhood” association,^27,28 mean-shift process,^29–31 the Hungarian method,³² Viterbi algorithm,³³ Kalman filtering,³⁴ particle filtering,³⁵ multiple hypothesis tracking (MHT),^18,36,37 and active structured learning,³⁸ can be used to determine the most likely cell correspondence between frames. These techniques are either applied to images of cells migrating on a 2-D surface or 3-D in vivo florescence images, which are far different from our 3-D in vitro application. In the works by Ong et al.,⁶ stained ECs and vessels were tracked using Bayesian filtering from time-lapse confocal images of 30-min intervals. However, our images are acquired once a day, which is much longer than most existing tracking systems. Therefore, tracking these cells is a challenge that has not been previously addressed.

1.2.

Summary of Contribution

In this paper, we focus on tracking by detection algorithm to track the migrating ECs within angiogenic vessels formed in 3-D MFDs by linking the coarse time-lapse 2-D phase-contrast images and 3-D end-point confocal images. Since the ECs in 3-D MFDs can appear and disappear from our 2-D phase-contrast images, the main challenge is to accurately identify the cell candidates and track cell division and migration in 3-D (in focus/out of focus) over image sequences.

Our proposed system includes image registration, vessel segmentation, cell detection, and multiple hypothesis Kalman filtering, to track cell proliferation and migration within angiogenic vessels in different frames. As mentioned, we obtain time-lapse phase-contrast images and end-point confocal images of the angiogenic vessels from the experiments. The stained cells are easily recognized in the end-point confocal images. Therefore, instead of starting with the first time point like other cell tracking applications, we apply backward Kalman filter by initializing our tracks with the accurate cell locations obtained from the confocal images at the final time point where all proliferation has occurred. Furthermore, we incorporate biological knowledge to add constraints and to estimate the track probability during cell association:

Since one phase-contrast image is recorded each day for each slot, biologists can conduct many experiments in parallel, where cell viabilities are well-maintained. Moreover, by backward tracking each cell in the 3-D network structure, cell lineage and trajectory information, linking the chemical and physical characteristics of a cell can be obtained in an automated manner.

The rest of this paper is structured as follows. In Sec. 2, we explain our proposed automated cell tracking system. In Sec. 3, we present numerical results on our experimental phase-contrast images. We offer concluding remarks in Sec. 4.

2.1.

Image Registration

We acquired the experimental images manually at different days; hence, the angiogenic vessels are misaligned across the image sequences as shown in Fig. 3. We spatially register all the image sequences by transforming them into a common coordinate system, so that we can associate and track all the cells.

Fig. 3

Experimental phase-contrast images of two successive frames. The trapezoidal post is designed to retain the gel within the microfluidic channel. Lines 1 to 4 can be detected and represented using Hough transform. The red points are the correspondence points used during image registration.

To register these images, we first extract salient features from all the images. We choose the intersections of the straight post edges (lines 1 to 4) as features, which are labeled in red in Fig. 3. To automatically detect these features, we apply the Hough transform to identify the straight lines representing the post edges. From these detected lines, the corresponding intersections between images are obtained. We then derive the scale, linear translation, and rotation parameters based on the intersections’ coordinates by iterative closest point to align all the images. The pixels below the identified post edges are marked to zero to avoid any post features to be detected as vessels or cells in the later stage.

2.2.

Vessel Segmentation

We aim to investigate the behavior of cells within the gel. Hence, we do not need to segment cells under the cell–gel interface [see Fig. 4(a)]. Based on the observations, we fit a parabola to estimate the cell–gel interface in the experimental images. The parabola equation has the following form:

Fig. 4

Examples of vessel segmentation: (a) cell–gel interface, (b) the segmented binary vessel shape, and (c) vessel shape with centerline.

Next, we convert the processed experimental images into binary images with the global image threshold computed using Otsu’s⁴⁰ method. We obtain the segmented binary vessel shape [see Fig. 4(b)] through the following series of morphological operations in MATLAB™:²²

These parameters are selected based on a set of training images, and moderate changes in these parameters do not lead to significant change in the binary vessel shapes.

We represent the segmented vessels shape with MAT by a centerline [see Fig. 4(c)] and radii.³⁹ We can reconstruct the original vessel shape for visualization with the obtained MAT data set. Since the cells typically remain within the vessel, we can eliminate cell candidates outside the vessel during cell detection.

2.3.

Cell Detection

CNN is a widely used deep learning approach in pattern-recognition and image-recognition problems.⁴¹ We construct a CNN architecture by stacking multiple layers (see Fig. 5) and train the CNN network to distinguish cell and noncell templates. The resulting CNN classifier is then used to detect the cell candidates from the time-lapse phase-contrast images.

Fig. 5

Architecture of the CNN classifier for cell detection.

We use the end-point confocal images, where we can easily recognize the cells, to select our cell templates to ensure they are informative. By aligning the phase-contrast images with their corresponding confocal images, we can distinguish cell coordinates and noncell coordinates in the phase-contrast images. These coordinates are cropped as our templates, as shown in Fig. 6. After rotating the templates by 90 deg, 180 deg, and 270 deg, respectively, flipping them horizontally and vertically, and adding Gaussian noise, we obtain 3228 cell templates and 3372 noncell templates to train the CNN classifier. Each template has a size of $100\times 100$, which is sufficiently large to fit each individual cell.

Fig. 6

Examples of (a) cell and (b) noncell templates from phase-contrast images to train the CNN classifier. Size of each template is $100\times 100$.

The architecture of a CNN varies depending on the types and numbers of layers included. Some common types of layers include convolution, ReLU, pooling, dropout, fully connected, softmax, and classification layer.

The convolution layers are the core building blocks of a CNN, which extract different features of the input.^41,42 An image becomes a stack of filtered images after a convolution layer. ReLU is the abbreviation of rectified linear unit, which implements the function $y=\max (x,0)$ to set any negative input value to zero. It increases the nonlinear properties of the overall network without affecting the receptive fields of the convolution layer.⁴³ The pooling layer serves to progressively reduce the spatial size of the representation, to reduce the number of parameters and amount of computation in the network and, hence, also to control overfitting. Max-pooling is the most common pooling operation, which partitions the input image into a set of nonoverlapping rectangles, and for each such subregion, outputs the maximum. After several convolution, ReLU, and max-pooling layers, the fully connected layer combines all of the features learned by the previous layers to classify the images.⁴⁴ A dropout layer, which randomly sets input elements to zero with a given probability, follows the fully connected layer to reduce overfitting.⁴⁵ The fully connected layer usually uses the softmax activation function (softmax layer) for a classification problem.⁴⁶ The classification layer is the final layer, which uses the probabilities returned by the softmax activation function for each input to assign it to one of the classes.

The types and number of layers included in a CNN classifier depend on the training data. To get the most suitable architecture of our CNN classifier, we vary not only the number of layers but also the parameters of each layer. For the data at hand, we apply fivefold cross validation to get the optimal CNN architecture. The training templates are randomly partitioned into five equal-sized subsamples, where one subsample is retained as the validation data for testing the classifier and the remaining four subsamples are used as training data. We repeat the cross-validation process for five times to obtain our CNN classifier, which yields a classification accuracy of 92.39%, as shown in Fig. 5. The first convolutional layer has 20 filters with a filter size of $5\times 5$, and the second convolutional layer has 20 filters with a filter size of $10\times 10$. The pool size of both max-pooling layers is $2\times 2$. The output sizes of the first and second fully connected layers are 10 and 2, respectively.

We then apply the resulting CNN classifier to classify the pixels in our experimental phase-contrast images through a sliding window manner. The centers of the clusters, which are classified as positive class, are selected as the position of cell nuclei.

2.4.

Multiple Hypothesis Kalman Filtering

We combine backward Kalman filtering with MHT to associate and track the cell candidates within the vessels over time. Instead of starting with the image at the first time point, we initialize our tracks with the end-point confocal image, which provide the accurate cell locations and where all proliferation has occurred. Furthermore, we combine biological knowledge, such as vessel information and cell proliferation rate, to add constraints when estimating the probabilities in cell association.

We denote each individual cell state by ${\mathbf{x}}_k={[x_k,y_k,{\dot{x}}_k,{\dot{y}}_k]}^T$, where $x_k$ and $y_k$ are the estimated cell centroid positions in the $x$- and $y$-axes, respectively, and ${\dot{x}}_k$ and ${\dot{y}}_k$ are their corresponding velocities. The system covariance at time $k$ is denoted by ${\mathbf{P}}_k$. The observation state is the cell position in $x$- and $y$-axes obtained in cell detection process, denoted by ${\mathbf{z}}_k={[x_k,y_k]}^T$.

The backward Kalman filter consists of two steps: backward prediction and update. The prediction equations for the prior state estimate ${\widehat{\mathbf{x}}}_{k-1|k}$ and corresponding system covariance estimate ${\mathbf{P}}_{k-1|k}$ are

We assume a constant velocity model; hence, we have

Backward Kalman filter with validation gating is a popular approach to update the cell state during cell association. In each image, we extract multiple cells, which will be updated to multiple tracks. To associate an observation to a particular track, we employ a validation gate as an association threshold of whether to accept or reject an observation-to-target association. Observations that fall within the validation gate will be updated. As shown in Fig. 7, the blue and red ellipse are the validation gate for cell 1 and cell 2, respectively. Specifically, we compute the Mahalanobis distance between each observation to each state estimate as

Fig. 7

Backward Kalman filter with validation gating in data association: blue and red ellipses are the validation gate for tracks 1 and 2, respectively. Track 1 is associated with observation 1 and track 2 is associated with observation 2 in this plot.

The alternative way for multiple observations problem is to apply MHT, where the association of observations can be delayed till sufficient information is available. The flowchart in Fig. 8 shows a diagram of the proposed cell tracking system.

Fig. 8

Flowchart of multiple hypothesis Kalman filtering. Cell state, system covariance, and probabilities of all the tracks are recursively updated over time to form multiple hypotheses. The cell associations are obtained by selecting the most probable hypotheses.

At the end point, which is the initial time point of the backward Kalman filter, we create one track for each cell. Each cell state and covariance estimate at time $k$ are backpropagated to time $k-1$ using Eqs. (2) and (3). If there are $N_{k-1}$ observations which fall within the validation gate, MHT updates all the observations to the same track by generating $N_{k-1}$ new branches to the track. Figure 9 provides an example of the MHT process for one cell in two successive time steps. Within the blue ellipse, the purple diamond represents the predicted position of one track at time $k-1$. The green triangles indicate the observations that may be associated with this track with probability $p_1$, $p_2$, and $p_3$, respectively. The red pentagon depicts the hypothesis (with probability $p_0$) that the track is coming from the cell–gel interface. All the possible hypotheses in step $k-1$ are maintained and considered independently to generate new hypotheses in the next time step, as indicated by the red ellipses in Fig. 9. Since the number of branches for each track increases exponentially, to keep the computation under control, we maintain $M$ branches for each time by including biological knowledge and pruning the tracks with low probability.

Fig. 9

MHT for one cell in two successive time steps.

Let us consider cell $i$ at time $k$. This cell might be due to the migration or proliferation of an existing cell labeled as $j$ at time $k-1$. It is also possible that this cell is in the cell–gel interface or out of focus (not in the focal plane) at time $k-1$. The corresponding probabilities are calculated in different ways.

We employ the Mahalanobis distance to calculate the probability that track $i$ stems from the migration or proliferation of an existing cell $j$. If the distance is smaller than the threshold gate value, the probability that the track $i$ is associated with observation $j$ is calculated as

We assume that the cells nearer to the cell–gel interface (larger $y$-position) are more likely to be new tracks. The probability that track $i$ is due to migration from the cell–gel interface is computed as

We include biological knowledge to reduce the number of the new tracks in each step. Since the doubling time of our ECs is 15 to 48 h, we assume that one cell can at most split once within one day. In other words, we allow at most two cells at time $k$ to be associated with one same observation at time $k-1$ during merging. In our angiogenic vessels, cells cannot overtake each other due to the existence of tight cell junctions, so we eliminate the tracks where overtaking happens during association. Meanwhile, we discard the tracks with low probability to keep the number of tracks generation computations under control.

As mentioned, we keep $M$ branches and obtain their corresponding cell state estimates and system covariance estimates with backward Kalman update for each time step. Typically, $M=4$ is selected since the remaining branches have negligible probabilities. Each of these new tracks is considered independently and used to generate new predictions for next time step. We repeat this tracking process until $k=1$ and obtain $N_0{M}^{t-1}$ hypotheses in total, where $N_0$ is the number of the cells obtained from the end-point confocal images and $t$ is the total time steps. Last, the tracks are selected as the most probable hypotheses.

3.1.

Image Registration Results

With our proposed registration algorithm, we aligned all the phase-contrast image sequences and their corresponding end-point confocal images. We overlay the end-point phase-contrast images and their corresponding end-point confocal images to illustrate the registration results (see Fig. 10).

Fig. 10

Illustration of image registration results: (a) and (d) are the end-point phase-contrast images, (b) and (e) are their corresponding confocal images, and (c) and (f) show the coregistered confocal and phase-contrast images.

Figures 10(a) and 10(d) are the original end-point experimental images obtained from phase-contrast microscopy. Figures 10(b) and 10(e) are their corresponding confocal images. Here, the blue blobs are the cell nuclei stained by Hoechst, and the green structures are the actin stained by Phalloidin. We can see that the size and position of the vessels are different in the phase-contrast images and confocal images. The aligned phase-contrast and confocal images are shown in Figs. 10(c) and 10(f), where the resized phase-contrast images are shown in the background. We marked the cell nuclei stained in the confocal images as green blobs for clearer identification. The angiogenic structures are well-aligned in Fig. 10, suggesting our image registration algorithm is accurate.

3.2.

Cell Detection Results

We consider the same examples as in Sec. 3.1 to illustrate our cell detection results. Figures 11(a) and 11(e) show the results for the proposed CNN-based cell detection algorithm. Figures 11(b) and 11(f) show the results from partial least square regression (PLSR) approach with 30 PLS components as explained in Ref. 47. Figures 11(c) and 11(g) show the cell detection results for principal component analysis (PCA) in combination with SVM classification. We also trained an SVM classifier based on the original high-dimensional template data to distinguish the cells from the background in our experimental images, for comparison. The results are shown in Figs. 11(d) and 11(h). The red stars $(*)$ label the location of the detected cell nuclei for all these images.

Fig. 11

Cell detection results of examples in Fig. 10 from different approaches: (a) and (e) from CNN classifier, (b) and (f) from PLSR,⁴⁷ (c) and (g) from PCA in conjunction with SVM, and (d) and (h) from an SVM classifier. Red stars $(*)$ indicate the locations of the detected cell nuclei.

To evaluate the performance of these approaches quantitatively, precision and recall measures are calculated. Precision is the fraction of detected cells that are actual cells, and recall is the fraction of the actual cells that are detected

Table 1

Confusion matrix for a binary classifier.48

To assess the different cell detection approaches, we also compute $F$-score $F_1$, defined as the harmonic mean of precision and recall

To validate the cell detection results, we use the end-point phase-contrast images, matching the estimated cell locations with their corresponding confocal images. We analyzed 18 end-point phase-contrast images in this paper. Since the phase-contrast image is a 2-D slice and the confocal image is in 3-D with multiple slices, we choose the ground truth in two ways. For the first case, the actual positives and negatives are extracted from all the slices of confocal images. However, since our phase-contrast image is in 2-D, we use only one slice of the confocal image to provide the ground truth for the second case. We choose the slice that has a similar focal plane as the 2-D phase-contrast image. The TP, FP, and FN values for the four approaches are provided in Table 2. Their corresponding precision, recall, and $F_1$ values are listed in Table 3.

Table 2

TP, FP, and FN values for different cell detection approaches.

Table 3

Classification performance for different cell detection approaches.

For case 1, all the precision values are above 80%, suggesting that most of the detected cells are actual cells. The recall for CNN reaches 85.7%, which is higher than for the other three approaches, indicating that CNN detects most of the actual cells.

The recall of CNN is limited to about 85.7%, which can be explained as follows. The height of our device is around $120\mu \mathrm{m}$ and there are cells migrating throughout this range. We only acquire the 2-D image of the angiogenic vessels in focus. Some cells, which are external to the vessels, are out of focus and almost invisible in the phase-contrast images. These cells are not detected by our algorithms. However, these out-of-focus cells are observed in the 3-D confocal images. This scenario can be seen in Fig. 10.

For case 2, the recall for all the four approaches increases as shown in Table 3. Particularly, the recall for the proposed cell detection algorithm (CNN) increases to 91.2%. From the values of $F$-score in Table 3, which balances recall and precision, we can also conclude that CNN is a preferable algorithm to detect the ECs within the 3-D angiogenic vessels from our experimental phase-contrast images. We hypothesize that the features extracted from CNN are more suitable to represent the cell templates.

3.3.

Tracking Results

Since our cell association approach is detection based, we use the cell detection results from the proposed approach based on CNN as the cell locations during cell association.

We provide two examples of the association results in Fig. 12. The results for the proposed tracking system are shown in Figs. 12(a), 12(b), 12(e), and 12(f). We also apply backward Kalman filter with validation gating for cell association [see Figs. 12(c), 12(d), 12(g), and 12(h)]. Each image presents the tracking results for two successive frames and the background is the experimental phase-contrast image. We can see that the proposed multicell tracking approach produces better tracking results than backward Kalman filter with validation gating.

Fig. 12

Cell tracking examples: (a)–(d) for example 1 and (e)–(h) for example 2. (a), (b) and (e), (f): backward Kalman filter with MHT, (c), (d) and (g), (h): backward Kalman filter with validation gating. Each image shows the outcomes of the tracking algorithms for two successive frames, where the phase-contrast image is shown in the background.

We applied the proposed system to track cells in five sprouting slots, where each slot has four to five frames (one frame per day). We estimated the association accuracies (the percentage of the correct associations in the total associations) for the proposed approach, the nearest-neighborhood association, the backward Kalman filter with validation gating, and the Hungarian method³² (see Table 4). We excluded in this comparison methods in the literature that require more information about the cells than merely the cell locations, since such information is unavailable for our data. Since these methods cannot be applied to our data, we did not assess them. We obtained the ground truth through manual association, which is the traditional method of analyzing the cell migratory behaviors.

Table 4

Accuracy (percentage of correct associations) for different cell association approaches. I, nearest-neighborhood association; II, Kalman filter with validation gating; III, Hungarian method; and IV, Kalman filter with MHT.

Since the sampling interval is one day, cells migrate over a relatively long distance between frames. Using nearest-neighborhood association, the new tracks at time $k$ are typically associated wrongly to cells that are near to the cell–gel interface at time $k-1$. In the backward Kalman filter with validation gating approach, we predict the track backward each time step via a constant velocity motion model. In this model, we assume that the track velocity is a constant during one time step, and we update it at every time step. This approach yields more reliable association results, improving the association accuracy from 26.7% to 60.1%. In the Hungarian method, we first predict the new positions of the cells at time $k$ using a constant velocity model and then estimate the Euclidean distances of the predicted cells at time $k$ and the cells at time $k-1$, to find the associations with minimal cost. However, cell migration and proliferation are not considered since it only allows for one-to-one matching. Thus, only 53.0% association accuracy is obtained. In our approach, we combine MHT with backward Kalman filter, which delays the decision making, and determine the associations once sufficient information is available. Moreover, biological knowledge is incorporated when considering the new cells from the cell–gel interface and reducing the number of new tracks. As a result, we can track not only the migration and proliferation of the existing cells but also the new cells from cell–gel interface and out-of-focus plane (see Fig. 12). This explains why the association accuracy increases to 86.4% for our approach. Compared with the association accuracy (85.9%) based on PLSR in Ref. 47, we can also conclude that the increment in cell detection accuracy due to the CNN improves the association accuracy.

Cell trajectory information can be obtained from the tracking results. Moreover, the system can automatically generate cell lineage plots (see Fig. 13), showing the history of the cell proliferation and cell migration into the gel, with timestamps of when it is in focus and out of focus. Such cell lineage plot provides crucial information (e.g., the number of cells, division time, growth fraction, etc.) to biologists who are interested in studying cell migration behaviors under different conditions.

Fig. 13

Cell lineage plot: (a) for example 1 and (b) for example 2 in Fig. 12. The blue lines represent the cells that migrate within the 3-D gel; the red stars represent the incoming cells from the cell–gel interface, whereas the green stars indicate the cells migrating out of focus.

The cell lineage plots in Figs. 13(a) and 13(b) correspond to the two examples in Fig. 12. The blue lines indicate that cells migrate within the 3-D gel. The branch points of these blue lines represent cell proliferation. The red stars illustrate the incoming cells from the cell–gel interface. The green stars indicate that cells migrate out of focus. The cell lineage plot illustrates the history of cell migration and proliferation in a compact and effective manner.

By comparing the cell lineage plots obtained under different experimental conditions, biologists can easily figure out the influences of the chemical factors on cell division and cell migration. For instance, Fig. 13(a) is from angiogenic experiment under positive sphingosine-1-phosphate (S1P) gradient, and Fig. 13(b) is from angiogenic experiment without S1P.⁹ By comparing these two figures, one can conclude that positive S1P gradient induces more recruitment of cells from the monolayer and stimulates cell proliferation. The cell lineage plots provide quantitative information about the influence of S1P on angiogenesis. Such quantitative data are an important step toward better understanding and new theories of angiogenesis; hence, the proposed tools probably will enable biologists to make new discoveries concerning angiogenesis. These data can also be used to develop mathematical models to predict the cell migration under different conditions, leading ultimately to a better understanding cell migration in angiogenesis.