Cell Segmentation and Tracking in Microfluidic Platform

In this research, we have concentrated on trajectory extraction based on image segmentation and data association in order to provide an economic and complete solution for rapid microfluidic cell migration experiments. We applied region scalable active contour model to segment the individual cells and then employed the ellipse fitting technique to process touching cells. Subsequently, we have also introduced a topology based technique to associate the cells between consecutive frames. This scheme achieves satisfactory segmentation and tracking results on the datasets acquired by our microfluidic platform.


Introduction
Cell migration plays an important role in many biomedical fields, such as drug test and disease diagnosis [1].Traditionally, cell migration is observed under Boyden chamber or Tran swell assays and other cell migration assays.However, these conventional methods lack of chemical gradients control and capability for quantitative analysis and often require large amounts of reagents and cell samples.Compared with the abovementioned methods, the microfluidic devices provide a more satisfactory platform for quantitative cell migration due to its capability of configuring precise and stable chemical concentration gradients, lower cells and reagents consumption, and the potential for high-throughput experiments [2][3][4].
High-throughput of images makes the manual observation of cell migration a labor-requiring and time-consuming process, and the accuracy of manual tracking highly depends on the experience and judgment of the individual researchers.Therefore, an effective automatic multiple objects tracking system is essential to conduct the quantitative analysis.
Most cell tracking techniques are composed of two phases, detecting and segmenting cells frame by frame, and then associating the detected same cells over two or more consecutive frames.A large number of segmentation methods have been introduced in the past decades, and many of them are still receiving intensive attention from medical image analysis community, such as watershed [5,6], edge detector [7], split and merge [7][8][9], re-gion growth [10][11][12], and some clustering methods [13][14][15][16].
In this research, we have employed the active contour model, which was first introduced by Kass et al. [17].Given an initial contour, this method would evolve towards image features such as object boundaries, and the evolution will continue until the energy functional reaches the minimal.However, the original model is parametric and would fail when topology changes happen [1].To handle the topology changes, Caselles et al. proposed the geodesic active contours with the flexibility of topology, which evolves the contours under a level-set framework [18].With the level-set framework, the cell tracking group in Carnegie Mellon University has proposed several important improvements of active contour models to distinguish touching cells, though these methods require large computational time and memory load [19][20][21][22][23].
The objects association has been a focus of researchers in some scientific fields such as radar system and video surveillance.Multiple Hypothesis Tracking (MHT) [24] and Joint Probabilistic Data Association Filter (JPDAF) [25] are two well-known examples of using Recursive Bayesian estimation, which is an effective method in objects association.While MHT needs the construction of exhaustive hypothesis set to select the optimal trajectory, even with the pruning techniques, this procedure requires substantial computation time and memory space [26].In contrast, JPDAF is a simpler and suboptimal approach that demands only fixed computational resources per iteration [27].Data association, in general, can be regarded as an optimal assignment problem and could be resolved by Hungarian algorithm [28].We refer to [29] for a background study on data association techniques.

Device Fabrication and Cell Preparation
We have employed the similar fabrication of microfluidic devices with the same standard soft-lithography protocol described in [30].The pattern was designed in a computer and printed into a transparency film to make a mask.A silicon wafer was coated with a thin layer of photo resist using a spinner.The master was finished by patterning the design on the wafer through the mask by UV processing.Liquid PDMS was poured on the master and cured in an oven.The PDMS replica was then peeled off and bonded to a glass slide by plasma treating to make the microfluidic device.The device was then coated by fibronection for one hour and blocked by BAS for another hour before starting the cell experiment to help cell bonding and migration on the substrate in the microfluidic channel.
The neutrophils were isolated from human whole blood using the gradient density centrifugation method.The cells were cultured in an incubator before loading in the microfluidic device.A 10 nM IL-8 gradient was generated in the microfluidic channel.The device was then put under the microscope and time-lapse image acquisition and further analysis was done by the custom-developed program.

Region-Based Active Contour Models
For an image function I(x,y) in the image domain, the CV model proposed by Chan and Vese [31] is: where C is any possible curve, inside(C) and outside(C) are two regions inside and outside the contours, and c 1 and c 2 are the average image intensity of inside(C) and outside(C), respectively.The first two terms in Equation (1) are called as "global fitting energy", which will have the minimum values if curve C is the real boundary of an object.
Since the CV model is piecewise constant and do not contain any local information, therefore, the optimal constants c 1 and c 2 might be significantly different from the real image data if the intensities inside or outside the curve C are inhomogeneous.Considering such a situation happens commonly when an image is captured by timelapse microscopy, the region-based active contour model is a more desirable adaptation in our research.In the region-based active contour model, the global fitting energy is replaced by a Region-Scalable Fitting energy (RSF), which contains local intensity information.Assume Ω 1 = outside(C) and Ω 2 = inside(C), the RSF for each pixel x ∈ Ω is defined as: where λ 1 and λ 2 are positive constants, and f 1 (x) and f 2 (x) are the approximate image intensities in Ω 1 and Ω 2 .The intensities I(y) are taken into account in the fitting energy come from the region centered at pixel x, the size of which is under the control of the kernel function K.A Gaussian kernel was chosen in [32].
Propagating the energy ε Fit x on the entire image, it derives the following: The second term in Equation ( 3) is a penalty term to smooth the contour C. Since this model is parametric, it is necessary to translate the parametric active contour to a geometric active contour, which is more desirable to deal with topology changes [18,33].By applying Heaviside and Dirac functions, numerical approximation of the evolution of the level set function can be written as: where f 1 and f 2 are For more details of the active contour model, we would refer to [32].

Splitting Touching Cells by Ellipse Fitting
In general, the segmentation of touching cells is a challenging task [34,35].When the traced cells enter into a blob, the boundaries of the contacting cells are blurred, and most segmentation algorithms would fail in finding the edges of cells in this situation.
In our research, we first select a region larger than the blob where touching happens, then apply the ellipse fitting [36,37] to estimate the features of contacting cells.obvious that the cells could be presented by two ellipses in this case.Therefore, splitting the cells is equivalent to estimate the parameters of the ellipses which best fit the real data.
To apply the ellipse fitting, first, a set of seeds for the closed contours produced by the Region-Scalable Fitting based Active Control model are generated by ultimate erosion, which is a binary operator in mathematical morphology.The ultimate erosion repeat eroding an object until the object disappears, while the residual points are considered as seeds.Figures 1(b) and (c) illustrate the beginning and result of an ultimate erosion process, respectively.For convenience, we denote the pixels on the contours as C(x c ,y c ), and the seeds as S i , where i = 1, 2, ..., N and N is the number of seeds.For each of seeds S i , we then sort C(x c ,y c ) into increasing order according to the distance from C(x c ,y c ) to S i .At the end of sorting process, there would be N increasing order lists of C(x c ,y c ), denoted as C i sorted (x c ,y c ).First M elements are selected from C i sorted (x c ,y c ) into a set C i 1 (x c ,y c ), and then incrementally append more elements into the set where k means at the kth stage.In each stage, an ellipse is fitted to the pixel coordinates in set C i k (x c ,y c ) by direct least squares fitting of ellipse [38].After processing all stages, a number of candidate ellipses, as shown in Figure 1(d), are produced and the best fitting ellipse will be selected from them.The obtained 4 best fitting ellipses of the seeds in Figure 1(c Two essential features are taken into account in our criterion of fitness, the first term of Equation ( 7) rewards the ellipse belonging to the region of object, while and the second term penalizes the ellipse out of the region of object.The ellipse with the highest value of this criterion will be chosen as the best fitting from a list of candidate ellipses.By performing the selection process to all of N increasing order lists, we would obtain N best fitting ellipses from N seeds.In the ultimate erosion process, since the number of seeds would likely be more than the touching cells, thus it is necessary to make a decision that which of the best fitting ellipses represents the true cell.As a solution, we have arranged the N best fitting ellipses in decrease order by the fitting criterion values, and eliminated all ellipses that have an overlap over 60% with the previously selected ellipse.The rest two ellipses representing the touching cells are shown in Figure 1(f).

Data Association Using Graph Theory
An association process after successful cell segmentations is to link the corresponding cells between two consecutive frames.In our cell migration experiment, we have focused on investigating the slower moving cells which are bonded to the glass substrate.Since these unadhered cells are not desirable and can be discarded, we employ a data association method based on graph theory [39], which is less costly in computational complexity and more suitable to address the adhered cells in our experiments.
The migration speeds of cells are different in our experiments.The cells moving faster are regarded as active cells, while the ones with smaller displacement between two consecutive frames are classified as lazy type.Zhang et al. has presented a real-life example to explain the idea of this approach [39].For example, if the neighbors 1, 2, 3, 4, and 5 of A and B in Figure 2(a) are already identified in Figure 2(b), then A, B and X, Y can be matched correctly according to the topology information of their identified neighbors.fied neighbors is an essential procedure to discriminate the un-matched objects.Naturally, the objects moving slower are more likely to be recognized, therefore a good association strategy is to identify these lazy type cells first.In the situation that cells almost remain their positions and shapes, a nearest neighborhood search is an effective method in despite of its simple basis.Assume that N cells (T i ; i = 1, ..., N) have been tracked up to the last frame, and M cells (C i ; i = 1, ..., M) have been detected in the current frame.A cost function is introduced here to present the similarity between T i and C i :

As the example shown in
where G i is the maximum Euclidean distance that a cell i can move between two consecutive frames, Ld i and Sd i are the maximum differences of perimeters and areas for the possible matched cells between two consecutive frames, respectively.If the distance d ij is greater than G i , then Cost ij is set to zero and the correspondence will be ignored.For each of cells in the current frame, we select the track with the highest Cost value among the N tracks in the last frame and then assign the cell the same label with the track.Since the assignment between a cell and a track is a one to one relationship, a process of optimizing is essential when more than one cell tend to be associated with the same track.Hungarian algorithm [28] is an effective solution to obtain the optimal assignment in our case, and is applied in this study.In our experiments, most of cells nearly maintain their positions between two frames.Therefore, lower threshold G i can significantly reduce the computational time since majority of cells are unnecessary to be considered.The marked lazy cells are presented in Figure 3.

Active Cells Association by Graph Theory
After having successfully tracked lazy cells, which can provide the essential topology information about the neighbors of the unmatched active cells, we now will focus on the unmatched cells pairs.The term of an unmatched cells pair represents an unmatched cell in the current frame and the one in previous frame.
Two phases are correlated with the linking of two unmatched cells pairs in consecutive frames.Firstly, a search region is assigned to each of the unmatched cells pairs.If the neighbors in the search region of an unmatched cells pair are the same in the consecutive frames, the two unmatched cells pairs are claimed as associated.In other words, the unmatched cells pair should have the same number of neighbors that are similar in directional positions.The directional position where F k and F k+1 denote the index of frames, R(i) and R ' (j) compose an unmatched cells pair, k is the index of the correlated neighbors in the neighborhood group, and Angle indicates the degree value of neighbor in a space centered on the unmatched cells pair. Figure 4 illustrates an example of Angle F k ,R(i) (k) which is represented by D(k).
Then, for the rest of unmatched cells, the neighbors of the unmatched cells pair with the same label are called as Share Source Neighborhoods (SSN).The likelihood of a cell pair is evaluated by [39]: ( ( ), '( )) where Q dist is the measurement of internal attributes of the pair, such as the similarity of size and location.Q SSN is used to measure the likelihood of the topology between their neighbors.Dist and Size are the differences of distance and size between R(i) and R'(j).α k and β k are the correlated neighbors from the set of SSN.Besides, constants δ D , δ S , and δ θ are used to set the sensitivity of Q to each factors.Cells pairs with largest Q values would be matched.This matching process will be repeated until no pair can be matched.

Experimental Results
In this research, we have applied our new Cell Segmentation and Tracking system to two different sets of data, which were recorded by our microfluidic platform.Table 1 illustrates that although temporal efficiency of the Region-Scalable Fitting (RSF) model is lower than those of classic algorithms, its accuracy after ellipse fitting is generally higher than those of other methods in our experiments.Considering our application is not realtime required, we select RSF model to provide better segmentation results.Compared with other methods, RSF model achieves not only higher segmentation accuracy but also better segmented contours.The results of our experiment, as presented in Figure 5, show that the watershed transform results in over-segmentation inside the cells; the graph cuts fails to detect some cells and some obtained contours are incorrect; the edge detector could not obtain closed contours if the edges are not salient; while the contours produced by the RSF model are smooth and closed.
The trajectories of DataSet 1 and DataSet 2 produced by data association based on graph theory are presented in Figures 6 and 7.The overall tracking accuracy on DataSet 1 and DataSet 2 are 86.7% and 91.04%, respectively.

Conclusions
In this research, we have conducted study on trajectory extraction based on image segmentation and data associ-    ation in order to provide a solution for rapid microfluidic cell immigration experiments.We applied region scalable active contour model to segment the individual cells and the ellipse fitting technique to process touching cells.
We have also introduced a topology based technique to associate the cells between consecutive frames.By applying our new Cell Segmentation and Tracking system to two different sets of data recorded by microfluidic device, we have came up with some encouraging outcome.The overall tracking accuracy on two sets of data is 86.7% and 91.04%, respectively.

Figure 1 (Figure 1 .
Figure 1.(a) Contour of two touching cells; (b) Binary image of touching cells; (c) Result of ultimate erosion; (d) The set of candidate ellipses; (e) Best fitting ellipses of all the seeds; (f) Best fitting ellipses after overlapping deletion.

where
) are presented in Figure 1(e).To rank the fitness of the ellipses, we have adopted the fol-Ellipse(C i k (x c ,y c )) denotes the ellipse fitted by the point C i k (x c ,y c ), ϕ object represents the touching objects, and a and b are the weights.

Figure 3 .
Figure 3. (a) All marked cells in the previous frame; (b) Marked lazy cells in the current frame.ismeasured by calculating the difference of degree values of each correlated neighborhood pair[39]:

Figure 4 .
Figure 4.The unmatched cells pair is centered in the axes.D(k) represents the angle of each correlated neighbor.

Figure 5 .
Figure 5.The comparison of results by different segmentation techniques on Dataset 1.(a) Segmented by watershed transform on the gradient magnitude; (b) Segmented by graph cuts; (c) Segmented by canny edge detector; (d) Segmented by RSF model.

Figure 6 .
Figure 6.The results of applying the new Cell Segmentation and Tracking system to DateSet 1: (a) At frame 3; (b) At frame 14; (c) At frame 38.

Figure 7 .
Figure 7.The results of applying the new Cell Segmentation and Tracking system to DateSet 2: (a) At frame 2; (b) At frame 22; (c) At frame 45.