Object Tracking and Data Analysis
with Time-Lapse Video of in Vitro
MTLn3 Cell Lines
Research Supervisors
Dr. Ir. Fons Verbeek
Imagery & Media
Section Imaging & Bioinformatics
Leiden Institute of Advanced Computer Science (LIACS)
Leiden University
Prof. Dr. Bob van der Water
Section Toxicology
Leiden/Amsterdam Center of Drug Research (LACDR)
Leiden University
Master Thesis Report
Kuan Yan
Section Imaging & Bioinformatics
Leiden Institute of Advanced Computer Science (LIACS)
Leiden University
Abstract
This paper considers the problem the improvement and application of the KDE Mean Shift
tracking algorithm and data analysis in a migration study of MTLn3 cells. The aim is to convert cell
migration videos into numeric description of changes in cell behavior. The choice of KDE Mean
Shift Tracking is based on its robust and prominent performance in time-lapse studies compared
to other algorithms. Morphology and motility measurements are defined by considering both
statistics and true biological translation. The detection of significant behavior change relies on
both feature selection and hypothesis testing in the feature space.
The feedback from the “wet-lab” to our tracking analysis and results indicate that the accuracy of
cell migration analysis is increased significantly and labor time is reduced enormously (over
300%). In addition, it is also believed that the current feature measurements reveal several
important behavior measurements, which cannot be determined by manual observation.
Supported by the results in cellular analysis we intend to extend our experiments into molecular
analysis such as focal adhesions in the near future.
This project is conducted under joint operation of Leiden Institute of Advanced Computer Science
(LIACS) Leiden University with Leiden/Amsterdam Center of Drug Research (LACDR) Leiden
University.
-2-
Table of Content
Abstract............................................................................................................................................ 2
Table of Content .............................................................................................................................. 3
1 Introduction ................................................................................................................................... 5
2 Video Data Analysis ..................................................................................................................... 7
3 Image Segmentation and Cell Labeling ....................................................................................... 9
3.1 Region Based Segmentation ............................................................................................... 10
3.1.1 Isodata Segmentation ................................................................................................... 10
3.1.2 Local Adaptive Segmentation ....................................................................................... 11
3.1.3 Region Based Segmentation in Cell Images................................................................. 12
3.2 Edge Based Segmentation .................................................................................................. 14
3.2.1 Derivative of Gaussian Edge Detector .......................................................................... 14
3.2.2 Laplacian of Gaussian Edge Detector........................................................................... 15
3.2.3 Canny Edge Detector .................................................................................................... 16
3.2.4 Edge Based Segmentation in Cell Images.................................................................... 17
3.3 Region Growth Segmentation.............................................................................................. 19
3.3.1 Region Growth Mechanism ........................................................................................... 19
3.3.2 Different Similarity Criterions......................................................................................... 20
3.3.3 Region Growth Segmentation in Cell Images ............................................................... 20
3.4 Object Labeling and Filtering ............................................................................................... 22
3.4.1 Region Based Labeling ................................................................................................. 22
3.4.2 Object Filtering .............................................................................................................. 22
4 Object Tracking Algorithms ........................................................................................................ 25
4.1 KDE Mean-Shift Tracking Algorithm .................................................................................... 26
4.1.1 KDE Mean-Shift using Steepest AscentOptimization ................................................... 26
4.1.2 KDE Mean-Shift Tracking in Cell Imaging..................................................................... 28
4.2 KDE Pairwise Alignment Tracking ....................................................................................... 31
4.2.1 Pairwise Alignment Algorithm ....................................................................................... 31
4.2.2 KDE Pairwise Alignment Tracking ................................................................................ 33
4.2.3 KDE Pairwise Alignment Tracking in Cell Image .......................................................... 34
4.2.4 Future Development...................................................................................................... 36
4.3 Overlap Ratio Based Tracking ............................................................................................. 37
4.3.1 Overlap Ratio Based Tracking ...................................................................................... 37
4.3.2 Overlap Ratio Based Tracking in Cell Image ................................................................ 37
4.4 Parametric Motion Model Tracking ...................................................................................... 40
4.5 Artificial Object Test ............................................................................................................. 41
4.5.1 Virtual Cell Factory ........................................................................................................ 41
4.5.2 Error Estimation in Artificial Object Test........................................................................ 45
4.5.3 Perturbation Study in Artificial Object Test.................................................................... 46
4.6 Measurements in Cell Migration Video ................................................................................ 49
4.6.1 Static Features .............................................................................................................. 49
4.6.2 Continuous Features ..................................................................................................... 52
4.6.3 Summarization Approaches .......................................................................................... 55
4.6.4 Noise Suppression in Measurement ............................................................................. 57
4.7 Measurement Bias in Time Resolution ................................................................................ 60
4.7.1 Time Resolution Experiment ......................................................................................... 61
4.7.2 Weighted Correlation..................................................................................................... 62
4.7.3 Sampling Rate............................................................................................................... 65
5 MTLn3 Cell Study........................................................................................................................ 67
5.1 High-Throughput Data Analysis ........................................................................................... 67
5.1.1 Feature Selection .......................................................................................................... 67
5.1.2 Statistical Hypothesis Test ............................................................................................ 68
5.1.3 High-Throughput Data Analysis in Cell Migration Study ............................................... 69
5.2 MTLn3 Cell Controlled Experiment....................................................................................... 72
-3-
5.2.1 Experiment Protocol ...................................................................................................... 73
5.2.2 Experiment Conclusion ................................................................................................. 73
5.3 Large Volume siRNA Screening using MTLn3 plGFP.......................................................... 75
6 Conclusion .................................................................................................................................. 77
7 Discussion .................................................................................................................................. 78
Appendix I: Workflow of Video Data Analysis ............................................................................... 80
Appendix II: Index List of Feature Measurements......................................................................... 81
Appendix III: Artificial Object Test Data Sheet .............................................................................. 82
III-1 Artificial Object Test and Perturbation Experiment............................................................. 82
Appendix IV: MTLn3 Cell Migration Data Sheet............................................................................ 84
Appendix V: Large Volume siRNA Screening Results .................................................................. 85
Appendix VI: Experiment Protocol of siRNA Screening ................................................................ 88
Reference ...................................................................................................................................... 89
-4-
1 Introduction
This report is contributed to the first stage of a cancer research project. This project is conducted
under collaborated operation of Section Imaging & Bioinformatics LIACS and Section Toxicology
LACDR, which is supervised by Dr. Ir. Fons Verbeek and Prof. Dr. Bob van der Water. The main
goal of this project is to establish the correlation between phenotypic behavior and focal adhesion
behavior.
Phenotypic research, contrary to genotypic research, is a study of how transfection affects
cytomorphology and cell migration. Compared with genotypic research, phenotypic research has
to major differences. First, genotypic behavior must be combined with special environmental
condition in order to result in certain phenotypic behavior. Second, certain phenotypic behavior
may be a result of multiple-transfection. However, sole study of phenotype may not reveal the
true mechanism of RNA translation. Thus, it must be combined with integrin protein.
Focal adhesion is believed to be an important integrin signaling kinase, which is found at the cell
[4]
membrane where the cytoskeleton interacts with proteins of the extracellular matrix . It is also
believed that a correlation existing between phenotypic behavior and focal adhesion behavior.
If such correlation can be established, it may reveal the molecular mechanism of cellular
communication. Current study of both phenotype and focal adhesion relies mostly on manual
observation and judgment. Therefore, object tracking is employed to provide objective
measurements and possibility of high-throughput.
Object tracking is a prosperous research field in computer vision. The time of its exact appearing
is ambiguous and oldest article can be dated back to early 1992. Object tracking has been
successfully applied in video surveillance [8], traffic control [9], military intelligence, etc. The
applications of object tracking in biological and medical domain appear in 1999[10]. Recent study
of blood cell experiment [12] shows a positive feedback to object tracking.
In this report, we will demonstrate the performance estimation of existing tracking algorithms in
complicated time-lapse video of living cells. The performance of each tracking algorithms is
tested by an artificial object test. The major function of this artificial test is to create objects that
simulating the cell migration. Tracking algorithms are applied to the video created from these
artificial objects. Tracking biases of each algorithm are studied under different parameter settings
for artificial objects. The tracking bias of the time-lapse recoding movie is done by specialist
visual confirmation.
Kernel Density Estimation (KDE) Mean Shift tracking shows nearly ideal performance in both
artificial tests and cell movies. In our artificial test, it has an accuracy of 96% in optimal parameter
settings, where cells are uniformly distributed, and 85% in random parameter settings, where
cells are randomly distributed. The tracking result of cell movies shows that KDE provides an
accuracy of 75%~85% by specialist visual confirmation. Measurements given by tracking
algorithm show a positive outcome in preserving global trends of cell migration.
Yet, the fundamental principle of current object tracking algorithms assumes a significant
discrimination between aligned subjects and non-aligned subjects while such trait is not
guaranteed in biological imaging. Therefore, our study also includes finding the factors that may
affect the tracking accuracy in biological domain while tracking bias or performance estimation is
the main criterion being involved.
With cell trajectories, feature measurements of cell migration can be extracted from raw video
data. These feature measurements are proven valuable in MTLn3 cell study. Several predefined
experiments such as EGFR upregulation and S178A downregulation are confirmed by current
measurements. By implementing feature selection and statistical test with current feature
-5-
measurement, full automatic and high-throughput data analysis of cell video is showing its great
potential in cell migration study.
The major sections of report are given as follows:
Chapter 2, Video Data Analysis, focuses on background information of time-lapse video analysis.
It describes the major image analysis principle, which will be involved in later section.
Chapter 3, Image Labeling, focuses on introducing popular segmentation algorithm. Image
labeling procedure is defined as the step to label all objects in image based on segmentation
result. Region based segmentation, edge based segmentation, and region growth segmentation
is introduced in this chapter.
Chapter 4, Object Tracking Algorithm, demonstrate our research of several existing tracking
algorithms and error estimations replying on artificial object tests. Tracking analysis of real video
(MTLn3 cell) is also given in this chapter.
Chapter 4, Section 4.1~4.4 focus on four tracking algorithms under different principles. It covers
1
affine transformation tracking , KDE mean shift tracking, overlap tracking, and KDE pairwise
tracking.
Chapter 4, Section 4.5, Measurements in Cell Migration Video, will introduce the motility and
morphology measurement employed in our tracking analysis.
Chapter 4, Section 4.6, Artificial Object Test, will describe the principle of artificial object test.
Under different parameter setting, three tracking algorithm will be compared. The exclusion of
affine transformation tracking is due to its limitation by definition, which is discussed in chapter 4.
Perturbation test on parameter setting demonstrates which factors may cause tracking bias.
Chapter 4, Section 4.7, Measurement Bias due to Time Resolution, is an important study of bias
due to sampling rate. Sampling rate is a critical factor in time-lapse video analysis study. This
chapter demonstrates how sampling rate affects measurements.
Chapter 5, section 5.1 MTLn3 Cell Test, shows several sample of applying tracking analysis in
MTLn3 cell migration study. Well-known experiments such as EGF stimulation, EGFR
upregulation, and S178A downregulation will be introduced and tracking analysis results are
studied and compared. The measurement results show nearly identical traits described in existing
literatures related to these experiments.
Chapter 5, section 5.2, High-throughput Data Analysis, is a study of full automatic cell video
analysis. With the help of tracking analysis and feature measurement, we are able to collect
numeric estimation of cell migration. Combing with statistical tools such as feature selection and
statistic test, significant cell behavior changes can be determined automatically.
Chapter 5, section 5.3, Large Volume siRNA Screening using MTLn3 plGFP, is a study of
applicability of tracking analysis. A demonstration of high-throughput data analysis with siRNA
screening experiment will be given in this section.
1
Also known as rigid object tracking
-6-
2 Video Data Analysis
Object or video tracking shows value applications in video surveillance, traffic control, military
intelligence, etc. Combing with pattern recognition and other statistic approaches, important
information can be extracted from video data. We now implement video data analysis for cell
migration study. In order to standardize procedure (see figure 1), we define the following steps as
the guideline for our research:
• Image Enhancing
• Image Segmenting
• Object Labeling
• Object Tracking
• Feature Measuring
• Data Analysis
Figure 1
Image Enhancing is the step to prepare raw image for further processing. Live cell image
contains a certain level of noise, which may not be easily improved by microscope settings.
Inconsistent staining is a common noise in fluorescence imaging. Figure 2 and Figure 3 shows
two cancer cells cropped from the same image. Threshold segmentation method may extract
foreground as figure 2, but not as figure 3. Only partial of cell in figure 3 has a significant intensity
value. Digital processing must be implemented in for further improvement.
-7-
Figure 2
Figure 3
Segmentation is an important step in image analysis. In order to label object from image, pixels
belonging to foreground must be first segmented from other pixels. Unlike segmentation in real
world image, microscopy image has a significant discrimination between pixels of foreground and
background. Our study shows that a sensitive segmentation method is the key to obtain
maximum tracking efficiency.
The major purpose of Object Labeling is to label isolated foreground region defined by
segmentation. Each object will be given a unique label. However, labeling method must be
consistent for each frame. The order of labeling may be critical for tracking algorithm such as
KDE pairwise tracking.
Object Tracking is the major aspect of research. The basic idea of tracking algorithm is to build up
alignment between objects in subsequential frames. Most tracking algorithms rely on similarity
measurements. If a strong connection can be built between two objects in subsequential frames,
an alignment will be made. However, high similarity of statistic measurement may not have
biology significant. The chance of find false candidates in cell video, which is due to high similarity
of objects, is considerable significant than any other tracking problem.
Feature Measuring is the purpose of tracking. If sensitive numeric measurement can be extracted
from cell tracking result, a more objective conclusion can be drawn. Compared with specialist
observation, which is relative subjective and doubetful, feature measurement based on video
tracking is indeed a major improvement in cellular and molecular video analysis.
Data Analysis is the step of finding implicit structure from feature measurements. Pattern
recognition and other statistic analysis are employed as major tools. We expect these feature
measurements may eventually establish the correlation between phenotypic behavior and focal
adhesion behavior.
-8-
3 Image Segmentation and Cell Labeling
The image segmentation is the operations to define the foreground from raw image. Identifying
foreground in image is non-trivial. There is no universal solution to all problems. Different
segmentation algorithm may only be suitable for several particular problems. When foreground is
defined, unique label is given to each isolated foreground region. Thus, potential regions of cells
can be derived from raw image. We will briefly introduce three types of image segmentation
algorithms.
• Region Based Segmentation
• Edge Based Segmentation
• Region Growth Segmentation
[16]
Region based segmentation algorithm is regularly used in microscopy image analysis . It
assumes that there exists a significant intensity boundary between foreground and background.
With a global threshold or local threshold value, foreground can be defined. Mechanism such as
iterative mean calculation of isodata threshold may provide automatic finding of intensity
threshold.
Edge based segmentation is a special segmentation method based on edge detection[16]. It
determines the significant gradient changes in intensity value. Laplacian filter of Gaussian edge
detection and Canny edge detection is commonly used. Based on a closed contour path, object
labeling can be given to each isolated contour path. However, convolution algorithm may result in
discontinuous edge. Thus, a path finding method is always needed for edge based segmentation.
Another approach based on pixel similarity is the region growth segmentation[16]. Region growth
algorithm intends to find similar pixels around seeds based on predefined similarity criterions. The
dynamic programming grows from one pixel to its neighbor pixels until criterions are no longer
fulfilled. The differentiation of gradient of intensity and intensity difference are frequently used as
similarity criterions.
When the image is segmented into foreground and background, each isolated region will be given
a unique label. The procedure of labeling is an important step before tracking.
-9-
3.1 Region Based Segmentation
Region Based Segmentation[16] is a segmentation technique given a predefined threshold value.
Pixel with intensity higher than threshold will be defined as foreground, otherwise background.
The binary labeling will be used for object labeling (see figure 4).
Pseudo Code: When foreground has a higher intensity than background
For each pixel in image
If current pixel has intensity larger than threshold value
Pixel is labeled as foreground
Else
Pixel is labeled as background
For each pixel
Label pixel as
background
false
Intensity>threshold
true
Label pixel as
foreground
Figure 4
3.1.1 Isodata Segmentation
Isodata segmentation [5][16] is defined as automatic threshold finding procedure (see figure 5). The
basic principle is identical as region based segmentation. The threshold is calculated based on
iteratively segmented image. A new threshold is defined as the mean value given average
intensity in foreground and background. The dynamic programming will stop when threshold
calculation is convergent.
Pseudo Code:
Given an initial threshold (also known as pervious threshold)
Do {
For each pixel in image
If current pixel has intensity larger than threshold value
Pixel is labeled as foreground
Else
Pixel is labeled as background
For each pixel in image based on labeling result
Calculate average intensity for foreground based on binary mask
Calculate average intensity for background based on binary mask
New threshold = (foreground average intensity + background average intensity)/2
} While (new threshold is not equal to pervious threshold)
- 10 -
Single Threshold Segmentation
Procedure on given threshold
All pixels are
labeled
Define new Threshold=(average
foreground + average background)/2
Calculate average intensity for
foreground and background
New threshold=Old threshold
false
true
Labeling
convergent
Figure 5
3.1.2 Local Adaptive Segmentation
Instead of looking for global significant[16], local adaptive segmentation uses neighbor pixels
around targeted pixel to find local significant (see figure 6). The segmentation will only be
performed for the center pixel.
Pseudo Code:
For each pixel in image
Define neighborhood region by 8-connectivity based on current pixel
Perform segmentation on neighborhood region
Labeling only current pixel by segmentation result
- 11 -
For each region
around pixel
Use threshold
from pervious
pixel
Monochromatic
Is Polychromatic
Polychromatic
Perform segmentation based on
neighborhood region
Labeling Center
Pixel
Figure 6
3.1.3 Region Based Segmentation in Cell Images
The table 1 gives a sample of converting a raw cell image into binary mask using region based
segmentation and its ramifications. Although, it appears that region based segmentation shows a
nearly prefect binary mask, it requires operator guiding for each images. Manually choosing
threshold is not considered as high-throughput method. Thus, automatic approach like isodata
segmentation is more preferred in cell image analysis.
Table 1
Raw Image
Region based Segmentation
Isodata Segmentation
Local Adaptive Segmentation with Isodata
- 12 -
- 13 -
3.2 Edge Based Segmentation
Edge based segmentation[16] follows the principle of finding isolated contour path (closed edge). A
binary mask of object can be defined by assign all pixels within its edges. Edge of object can be
determined by edge detection method.
An edge is a set of connected pixels that are on the boundary between two regions. Edge
detection assumes that significant gradient changes will occur at edge (regularly in intensity
space). However, edge detection is susceptible to noise. Noise suppression is required before
edge detection can be performed. Gaussian blurring is a frequently used noise suppression
procedure in edge detection. We will briefly introduce the following edge detection methods:
• Derivative of Gaussian (DoG) Edge Detector
• Laplacian of Gaussian (LoG) Edge Detector
• Canny Edge Detector
Still, edge detection procedures do not guarantee a continuous contour around object. An edge
linking operation (path finding) is needed when edge detection procedure is accomplished.
[19]
Different edge linking schema can be used and we will use Boie-Cox algorithm and hard c[20]
Mean algorithm.
3.2.1 Derivative of Gaussian Edge Detector
DoG edge detector[12][16] is a gradient-based edge detection procedure. With convolution
operation, DoG emphasizes pixels with a significant gradient change. The convolution operator is
the spatial approximation of frequency multiplication operator in frequency domain. Wrapper
frequency (also known as masking frequency) must be transformed into spatial matrix based on
either square distance or circle distance. Convolution operator[16] is required for responding
calculation.
Pseudo Code:
Convolution operation is performed for each pixel
Each significant responding of DoG is a potential edge
Magnitude of x-direction and y-direction is calculated
For each magnitude of pixel
If magnitude is higher than threshold
It is an edge
Else
It is not an edge
Each discontinuous contour is connect by edge linking schema
The filter can be described as follow:
G (r ) =
1
2πσ
−
e
DoG → G ' (r ) =
r2
2σ 2
,r =
−r 2
−r
2πσ
x2 + y2
3
e 2σ
2
- 14 -
Figure 7
Edge pixels can be determined by threshold segmentation on filtered image. A common approach
is to label all pixels as edge when its edge strength is above a threshold (see figure 7).
3.2.2 Laplacian of Gaussian Edge Detector
LoG edge detector[3][12][16] is designed for edge detection based on LoG responding (see Figure 8).
The simplest edge detection calculation is to threshold the LoG output at zero crossing (see
Figure 9). These boundaries can then be easily detected and marked in single pass.
Pseudo Code:
Convolution operation is performed for each pixel
Each significant responding of LoG is a potential edge
For each responding of pixel
If responding is a zero-crossing
Pixel is an edge
Else
Pixel is not an edge
Each discontinuous contour is connect by edge linking schema
The Laplacian filter of Gaussian kernel defined as:
G (r ) =
1
2πσ
−
e
r2
2σ 2
,r
= x2 + y2
∂ I ∂ I
L(x, y ) = 2 + 2 , LoG → G" (r ) =
∂x
∂y
2
2
−1
2πσ 3
Figure 8
Response from LoG (see Figure 9)
- 15 -
−r2
e
2σ 2
⎡
r2 ⎤
⎢1 − 2 ⎥
⎣ 2σ ⎦
Figure 9
3.2.3 Canny Edge Detector
Canny edge detector[3][12] [16] [18] is a robust method, which widely used in autonomous system. Its
aim was to discover the optimal edge detection algorithm. In this situation, an "optimal" edge
detector means:
•
•
•
Good detection - the algorithm should mark as many real edges in the image as possible.
Good localization - edges marked should be as close as possible to the edge in the real
image.
Minimal response - a given edge in the image should only be marked once, and where
possible, image noise should not create false edges.
To satisfy these requirements Canny used the calculus of variations - a technique which finds the
function which optimizes a given functional. The optimal function in canny edge detector is
described by the sum of four exponential terms, but can be approximated by the first derivative of
a Gaussian.
Step 1: Noise Suppression
Because the Canny edge detector uses a filter based on the first derivative of a Gaussian, it is
susceptible to noise present on raw unprocessed image data, so to begin with the raw image is
convolved with a Gaussian filter. The result is as a slightly blurred version of the original, which is
not affected by a single noisy pixel to any significant degree.
Step 2: Calculating intensity gradient of the image
An edge in an image may point in a variety of directions, so the Canny algorithm uses four filters
to detect horizontal, vertical and diagonal edges in the blurred image. The edge detection
operator (Roberts, Prewitt, Sobel for example) returns a value for the first derivative in the
horizontal direction (Gy) and the vertical direction (Gx). From this the edge gradient and direction
can be determined:
2
1⎤
⎡1
⎡1 0 − 1 ⎤
G = G x2 + G y2 , G x = ⎢⎢2 0 − 2⎥⎥ * A , G y = ⎢⎢ 0
0
0 ⎥⎥ * A
⎢⎣− 1 − 2 − 1⎥⎦
⎢⎣1 0 − 1⎥⎦
⎛G ⎞
θ = arctan⎜⎜ x ⎟⎟
⎝ Gy ⎠
- 16 -
The edge direction angle is rounded to one of four angles representing vertical, horizontal and the
two diagonals (0, 45, 90 and 135 degrees for example). An edge point is defined to be a point
whose strength is local maximum in the direction of gradient.
Step 3: Non-maximum suppression
The edge points determined in Step 2 give rise to ridges in the gradient magnitude image. The
algorithm then tracks along the top of these ridges and sets to zero all pixels that are not actually
on the ridge top to give a thin line in the output, a process known as nonmaximal suppression.
The ridge pixels are then thresholded using two thresholds, T1 and T2, with T1<T2. Ridge pixels
with values greater than T2 are said to be “strong” edge pixels. Ridge pixels with values between
T1 and T2 are said to be “weak” edge pixels.
Step 4: Enclosing Contour
Finally, the algorithm performs edge linking by incorporating the weak pixels that are 8-connected
to the strong pixels.
3.2.4 Edge Based Segmentation in Cell Images
The segmentation results of three methods are shown in Table 2. DoG and LoG shows less
desired results. Edge detection methods are too sensitive to noise. Small structure is a common
noise that requiring image to be blurred before edge detection. In addition, the structure of cell
causes the detection of extra edge, which is in fact not required. Canny edge detector is more
resistant to noise with its non-maximum suppression. Less significant gradient change will be
removed.
We believe that canny filter may be a robust approach for segmentation problem. Still, different
edge linking schema may result in different contours or even introduce further error, which is
beyond the scope of this report.
Table 2
Original Image
Edge by DoG Filter
Edge by LoG Filter
Edge by Canny Filter
- 17 -
- 18 -
3.3 Region Growth Segmentation
3.3.1 Region Growth Mechanism
Region growth segmentation assumes a certain level of similarity among pixels belonging to
same object. The algorithm will start from one prior pixel and grows until predefined similarity
criterion is no longer fulfilled[12][16]. The similarity criterion (also known as homogeneous criterion)
is regularly defined on pixel value considerations or on the anticipated size or shape of the
object[12]. Figure 10 shows a successful segmentation sample on right kidney from a CT scan
result based on intensity similarity.
Figure 10
Pseudo Code:
Interface Homogeneity_Detector
If condition fulfills
Return True
Else
Return False
Function Region_Growth
With given seed position
For each pixels in neighbor
If Homogeneity_Detector (current pixel) is fulfilled
Pixel is given the same label from seed
Perform Region_Growth (current pixel)
Else
Stop
Moving to next pixel
- 19 -
Stop growing if no pixel fulfills similarity criterion
Move to a new isolated seed
To locate the start position (seed) automatically is non-trivial in cancer study. Moreover, starting
from single pixel may not provide too many information in cell image. In our application, region
growth segmentation will use a less optimal threshold segmentation result as an initial scenario
and starts growing around the known region.
3.3.2 Different Similarity Criterions
Several techniques are derived from basic intensity similarity criterions include:
• Intensity difference
• Intensity gradient
• Probability threshold
Intensity difference based similarity criterion assumes that all pixels in interesting region share a
range of intensity value. Any pixel connected to initial region will be given the same label if its
intensity is in the predefined range (see Figure 6)[12].
Pseudo Code:
Function Intensity_Range (current pixel) interface Homogeneity_Detector
If intensity value for current pixel is in given range
Return true
Else
Return false
Intensity gradient is a hybrid technique relies on edge detection and intensity information.
[14]
Significant gradient change suggests a nature boundary of interesting region .
Pseudo Code:
Function Intensity_Range (current pixel) interface Homogeneity_Detector
If intensity gradient value for current pixel is in given range
Return true
Else
Return false
Probability threshold cuts off pixels with probability estimation lower than threshold given density
estimation in feature spaces[15]. A predefined initial region is required. In addition, a density kernel
must be defined.
Pseudo Code:
Function Intensity_Range (current pixel) interface Homogeneity_Detector
If probability value for current pixel is larger than given threshold
Return true
Else
Return false
3.3.3 Region Growth Segmentation in Cell Images
Region growth segmentation requires pixels belonging to same object share a strong feature,
such as very similar intensity value. However, cell image has many noises, which is difficult to
find such strong feature. Our current implementation of region growth segmentation using
gradient changes shows no significant improvement for segmentation result (see table 3).
- 20 -
Table 3
Original Image
Local Adaptive with Isodata
Region Growth Segmentation
- 21 -
3.4 Object Labeling and Filtering
3.4.1 Region Based Labeling
Segmentation procedure results in a binary image that representing foreground and background.
Yet it cannot be used for tracking. Each object in foreground must be identified. A common
approach is to give unique label to each separated foreground region (or interested region). Our
version of object labeling is given as following (see figure 11):
Pseudo Code:
Region_Based_Growing_Labeling(a pixel, a link list)
If current pixel is a foreground
Register the current pixel to current link list
Set current pixel to background
Image_Labeling( a binary image)
For each pixel in segmented image
If current pixel is a foreground
Region_Based_Growing_Labeling(current pixel, a link list for object)
For each pixel around current pixel
The current pixel is set to neighbor pixel and labeling is perform recursively
If no neighbor pixel is foreground
Current link list is terminated and a new object is created
Figure 11 A sample of labeling program growing from one pixel to the rest
However, due to intensity and many other factors, segmentation regularly results in a less desired
result of foreground. Therefore, objects extracted from segmented image will be filtered and all
abnormal objects will be removed before tracking starts.
3.4.2 Object Filtering
According to the summarization of practical error from segmentation result, we define the
following rule of object filtering. For clarification, some filtering rules are implemented during the
performance of tracking while others are performed right after segmentation (see figure 12):
- 22 -
Figure 12 a flow chart of filtering
Filtering after Segmentation
Small objects are often dirt or died cells. It has no research value and causes bias in feature
measurement. We have encountered both immobile and fast moving noises in experiment
(section 5.2). They both result in an inconsistent measurement against initial experiment design.
Therefore, they must be filtered.
Extreme large objects are commonly clustered cell due to less-optimal segmentation result. Using
separation method such as watershed, extreme large object may be avoided. Yet, the result of
separation may consist of bias due to different distance transformation.
Cells touching the edge of image are deformed. In addition, these cells are moving either away or
into the screen. It leads to a dilemma for tracking algorithm due to both deformation and
unpredictable direction of motion. To avoid tracking bias, cell-on-edge is often removed.
Filtering during Tracking
Due to high similarity between objects in microscopy image, it is essential to minimize the number
of possible candidate. This is accomplished by implementing candidate filtering. Based on our
research on artificial object (section 4.6) and MTLn3 (section 5.2), we notice that the most efficient
way to limited the candidate is to set a distance radius around the predecessor cell. A predefined
distance is given for each aligning. All cells that is outside of radius will be discarded from
candidate.
- 23 -
Filtering mechanism shows clearly improvement on tracking accuracy, yet it may also result in
problems such as early termination of trajectory or measurement bias (see chapter 4). Since bias
of both tracking and measurement are dynamically linked, an interesting research may be
conducted on the optimal seeking of both tracking and measurement. In next chapter, we will
introduce tracking algorithm.
- 24 -
4 Object Tracking Algorithms
In order to perform phenotypic analysis with time-lapse video (see introduction), object tracking is
employed. The fundamental principle of object tracking is to build up alignment between two
objects in subsequential frames. Therefore, most tracking algorithm relies on similarity
measurement. This chapter will introduce four tracking algorithms based on different principle of
similarity measurements. These algorithms include:
• KDE Mean Shift Tracking Algorithm[1]
• KDE Pairwise Alignment Algorithm
• Overlap Ratio Tracking Algorithm
• Affine Transformation Based Algorithm[2]
However, live cell video does not share several common characteristics of real world video. Thus,
the performance of each tracking algorithm in time-lapse video will be compared to find a robust
and sensitive tracking algorithm in live cell video. In addition, result of the chosen tracking
algorithm will be demonstrated in this chapter.
- 25 -
4.1 KDE Mean-Shift Tracking Algorithm
Kernel Density Estimator with Mean Shift Tracking[1] is a popular kernel based tracking algorithm
in computer vision. Kernel Based Mean-Shift tracking is considered as real world application of
localization approach. It is assumed next the support of desired modules, which should provide
(also see figure 13) (a) detection and localization in the initial frame of the objects to track
(targets), and (b) periodic analysis of each object to account for possible updates of the target
modules. The KDE Mean Shift Tracking consists of two steps:
1. The gradient of the nonparametric density function is estimated by a predefined kernel
density estimator.
2. Starting from each sample point, the mean shift procedure iteratively finds a path along
the gradient direction away from the valleys and towards the nearest peak.
Figure 13 Workflow Diagram of KDE Mean Shift Tracking Algorithm
The nearest peak, which is the maximum possibility estimation based on nonparametric density
estimator (model), is the theoretical center of true aligned object. Object, whose center of mass is
closest to the peak, will be the true aligned object. The method is introduced given as following.
4.1.1 KDE Mean-Shift using Steepest Ascent Optimization
Given n data points xn in the d-dimensional space Rd, the kernel density estimator with kernel
function K(x) and a window bandwidth h, is given by
1
fˆn ( x ) = d
nh
n
⎛ x − xi ⎞
⎟
h ⎠
∑ K ⎜⎝
Where the d-variate kernel K(x) is nonnegative and integrates to one. A widely used class of
kernels is the radially symmetric kernels
- 26 -
( )
K ( x ) = c k ,d k x
2
Where the function k(x) is called the profile of the kernel Normalization constant
c k ,d makes K(x)
integrate to one. The density estimator can be rewritten as
c k ,d
fˆh ,k ( x ) = d
nh
Where
⎛ x − xi
⎜
k
∑
⎜ h
i =1 ⎝
n
2
⎞
⎟
⎟
⎠
c k ,d is the normalization constant
Commonly used kernels are the Epanechnikov kernel
(
⎧ 1 −1
⎪ c d (d + 2) 1 − x
K E (x ) = ⎨ 2
⎪0
⎩
2
)
0 ≤ x ≤1
x >1
And the multivariate Gaussian kernel
K N ( x ) = (2π )
−
d
2
∑
−
1
2
e
1
− ∑ −1 x
2
2
∑ is the covariance coefficient of module.
, where
The standard mean shift algorithm is a steepest ascent procedure, which requires estimation of
the density gradient:
⎛ x − xi
∇fˆh ,k ( x ) = d + 2 ∑ ( xi − x )g ⎜
⎜ h
nh
i =1
⎝
⎤
⎡ n
⎛ x − xi 2 ⎞
⎟
⎥
⎢ ∑ xi g ⎜
⎜ h
⎟
i =1
⎢
⎝
⎠ − x⎥
= c k ,d fˆh ,G ( x )⎢
⎥
2
n
⎛
⎞
−
x
x
⎥
⎢
i
⎜
⎟
⎥
⎢ ∑ g⎜ h
⎟
⎠
⎦
⎣ i =1 ⎝
2c k , d
n
2
⎞
⎟
⎟
⎠
, where g ( x ) = − k ' ( x ) can in turn be used as
profile to define a kernel G(x). The kernel K(x) is called the shadow of G(x).
fˆh ,G ( x ) is the density
estimation with the kernel G. c k , d is the normalization coefficient. The last term is the mean shift
⎛ x − xi 2 ⎞
∑i =1 xi g ⎜⎜ h ⎟⎟
⎠−x
⎝
m( x ) =
2
⎛ x−x ⎞
n
∑i =1 g ⎜⎜ h i ⎟⎟
⎠
⎝
n
which is proportional to the normalized density gradient and always points toward the steepest
ascent direction of the density function. The standard mean shift algorithm iteratively performs
•
( )
Computation of the mean shift vector m x
k +1
k
( )
• Updating the current position x
= x +m x
Until reaching the stationary point, which is the center of candidate cluster (See Figure 14).
k
- 27 -
k
Figure 14
4.1.2 KDE Mean-Shift Tracking in Cell Imaging
With object annotation, potential region of cells are segmented into binary mask and labeled with
unique number. Each cell is converted into nonparametric kernel estimator (module). Steepest
descent is performed on the gradient space of density estimation. The illustration of tracking one
cell in two subsequential images is given in figure 15 and 16. The initial cell is segmented from
raw image and converted into density estimator.
Each extracted cell will be converted into a three dimensional module. Each dimension is defined
as (1) All x-coordinate of binary mask (2) All y-coordinate of binary mask (3) All intensity value of
binary mask. The probability estimation is calculated based on following equation. For detail,
please prefer section 4.1.1 KDE Mean-Shift using Steepest Ascent Optimization.
c k ,d
fˆh ,k ( x ) = d
nh
⎛ x − xi
k⎜
∑
⎜ h
i =1 ⎝
n
2
⎞
⎟
⎟
⎠
Object to be tracked
The next frame image (using raw image)
Figure 15
Figure 16
The figure 17 is the surface plotting of the possibility estimation given a normalized model, which
is converted from one labeled cell. Each peak represents the density estimation for each
candidate by given module. The closest and highest peak will be the true alignment. However,
normalized models are not implemented in our application. The reason is that normalization of
- 28 -
initial models will introduce remote candidate. In fact, the idea of mean shift is to measure the
closest shifting of distribution. Thus, local optimization is preferred.
Figure 17
The table 4 shows two samples of cell tracking in (a) video with low cell density and high
magnification level (b) video with high cell density and low magnification level. Trajectories of two
samples are approaching to 90% of tracking accuracy according to specialist visual observation.
Table 4
First Frame of Video (a)
First Frame of Video (b)
Last Frame of Video (a)
Last Frame of Video (b)
- 29 -
Trajectories of (a)
Trajectories of (b)
- 30 -
4.2 KDE Pairwise Alignment Tracking
KDE Mean Shift tracking algorithm may result in combining of two cell trajectories into one. Such
error will occur when two cells are nearly identical. In addition, cell migration, such as cell
differentiation, will separate one object into two true alignments that both of them must be tracked.
KDE Mean Shift tracking can only have one true candidate by definition, thus cannot solve such
problem. In order to solve both of these potential problems, we propose an experimental tracking
algorithm using pairwise alignment principle. The algorithm can be illustrated as following.
#pseudo code
Labeling all objects m in each frames
For each m in frames
Calculate sum density pairwise distance in frame m and m+1
Perform backtracking to find optimal through pairwise distance matrix with largest total
score
Adding new path
The current KDE pairwise alignment relies on cumulative density estimation as its distance table.
However, the chance of two identical objects is much higher in cell tracking, thus statistical
significance may not represent biological meaning. In KDE pairwise alignment tracking, this
drawback is becoming a decisive issue.
4.2.1 Pairwise Alignment Algorithm
Pairwise alignment algorithm [5] (see Figure 19) uses a distance transformation matrix of two sets
of objects. The distance may not necessary to be the Euclidean distance. It can also be KDE
(kernel density estimation, MSE (mean square error) estimation, or other distance functions. The
pairwise algorithm aims to find out the alignment path with highest score (prefect matching). We
will briefly explain the pairwise alignment algorithm in this section.
For gene sequence matching, a global pairwise alignment[5] can be calculated (also see figure 18):
Figure 18
INPUT: Two sequences S = s1...sn and T = t1...tm (n and m are approximately the same)
- 31 -
QUESTION: Find an optimal alignment (see figure 16).
Notation: Let σ(a, b) be the score (weight) of the alignment of character a with character b
(including spaces).
Lemma: Let V(i, j) be the optimal alignment score of S1...i and T1...j (0 ≤ i ≤ n, 0 ≤ j ≤ m).
V(A,B) has the following properties:
Recurrence Programming
for (1 ≤ i ≤ n,1 ≤ j ≤ m )
⎧V (i − 1, j − 1) + σ (S i , T j )
⎪
V (i, j ) = max ⎨V (i − 1, j ) + σ (S i ,− )
⎪V (i, j − 1) + σ (−, T )
j
⎩
Given σ (Si , T j ) is a mismatching penalty, it will be zero for a true matching
σ (Si ,− ) is a deletion penalty, it stands for a deletion of nucleotide
Given σ (−, T j ) is an insertion penalty, it stands for an insertion of nucleotide
Given
In normal application, a long sequence of deletions or insertions will be given a sequential gap
penalty, which using an incremental penalty calculation. The sequential gap penalty is calculated
as following, where open is the gap opening penalty, ext is the gap extension penalty, and n is
the number of gap extended.
penalty = open + ext ⋅ n
If no incremental gap penalty is used, the pairwise alignment will provide global alignment,
otherwise local alignment[5]. Local alignment is a special type of pairwise alignment, which long
sequence of gaps is discouraged. Global alignment will always return a sequence with the same
length of targeted sequence. Local alignment will return a sequence with the highest matching
score with a minimum number of gaps.
Figure 19
- 32 -
In order to find a prefect match, the backtracking algorithm is used to find the path with maximum
match score in distance table. The backtracking[13] is described as following (see figure 20):
Figure 20
1. Starting at Root, your options are A and B. You choose A.
2. At A, your options are C and D. You choose C.
3. C is bad. Go back to A.
4. At A, you have already tried C, and it failed. Try D.
5. D is bad. Go back to A.
6. At A, you have no options left to try. Go back to Root.
7. At Root, you have already tried A. Try B.
8. At B, your options are E and F. Try E.
9. E is good
4.2.2 KDE Pairwise Alignment Tracking
As pairwise algorithm used in gene matching, distance table and gap penalty is introduce for KDE
Pairwise Tracking. However, local alignment is preferred in current situation while global
alignment introduces many high matching with no biological meaning. In high identical object set
such as cell, the chance of finding two identical objects is relatively higher than other tracking
scenario. Table 5 shows a sample of cumulative density estimation. The density estimation is
calculated based on kernel density estimator.
c k ,d
fˆh ,k ( x ) = d
nh
Table 5
obj1
obj1 7.8
obj2 5.5
obj3 4.3
obj4 5
obj5 3.1
⎛ x − xi
⎜
k
∑
⎜ h
i =1 ⎝
n
obj2
7.5
6.5
4.4
5.1
3.1
obj3
6.5
6.8
4.6
5.3
3.1
2
obj4
5.8
6.5
4.4
5.3
3.1
⎞
⎟
⎟
⎠
obj5
5.4
5.3
4.3
5.1
5.1
1
(x ) = (2π )− 2 ∑ − 2 e
d
where K N
obj6
3.2
3.5
4.3
5.1
3.1
obj7
3.2
3.3
4.3
5.1
3.1
1
− ∑ −1 x
2
2
obj8
3.2
3.3
4.3
5.1
3.1
For pairwise alignment tracking, the gap penalty is calculated as following. Unless for very high
diversity, there will be very little double matching in final sequence. Therefore, we discourage the
gap by adding incremental penalty for each extension of gaps.
penalty = ext n
- 33 -
Table 6
All Trajectory
(Trajectory of object 3)
obj3 trajectory
160
140
120
100
80
obj3 trajectory
60
40
20
0
0
(Trajectory of object 4)
100
200
300
400
500
600
700
800
(Trajectory of object 5)
obj5 trajectory
obj4 trajectory
250
350
300
200
250
150
200
obj4 trajectory
obj5 trajectory
150
100
100
50
50
0
0
0
200
400
600
800
1000
1200
1400
1600
0
100
200
300
400
500
600
700
800
900
1000
The experiment (table 6) suggests that pairwise may still introduce remotely located but identical
object. Therefore, a strong distance penalty is introduced for remote objects. Since cumulative
density estimation is always ranged from 0 and 1, the distance penalty will set -1 to all objects,
which is too far from current object.
⎧kde( x, xi ) where x i is the model for object i
⎩− 1 if distance of object is higher than threshold
ωi (x ) = ⎨
Sample of distance limited table (see Table 7)
Table 7
obj1 obj2 obj3 obj4 obj5 obj6
obj1 7.8
7.5
6.5
5.8
-1
-1
obj2 5.5
6.5
6.8
6.5
5.3
-1
obj3 -1
4.4
4.6
4.4
4.3
-1
obj4 -1
5.1
5.3
5.3
5.1
5.1
obj5 -1
-1
-1
-1
5.1
-1
obj7
-1
-1
-1
-1
-1
obj8
-1
-1
-1
-1
-1
4.2.3 KDE Pairwise Alignment Tracking in Cell Image
KDE Pairwise Alignment algorithm shows potential in preserving cell migration information (see
table 8). Compared with KDE Mean Shift algorithm, KDE Pairwise Alignment algorithm shows
slight improvement when cells are dividing or escaping from screen. However, our
implementation suggests that the runtime of KDE Pairwise Alignment algorithm is significantly
higher than KDE Mean Shift tracking, which by average is nearly 40 times. It may be a major
drawback of this approach.
- 34 -
Table 8
First Frame of Video (a)
First Frame of Video (b)
Trajectories of (a) KDE Mean Shift
Trajectories of (b) KDE Mean Shift
Trajectories of (a) KDE Pairwise
Trajectories of (b) KDE Pairwise
- 35 -
4.2.4 Future Development
Pairwise alignment algorithm relies on distance matrixes (mutation table). To find such score
system is nontrivial in video data. Popular mutation table for gene matching such as BLOSUM62
is calculated based is on possibility of each nucleotide base appearing in two manually picked
sequences with 62% of similarity. However, mutation in video is highly dynamic. Such statistic
score requires an iterative calculation based on each alignment.
Yet, recursively calculation of mutation table increases computation costs. Our study suggests
that even without considering multiple alignment, the average tracking time of pairwise alignment
is nearly 40 times of KDE mean shift. It may be the major disadvantage of pairwise alignment.
Further study of score system will be the major direction of pairwise alignment tracking.
Regardless of its immaturity, pairwise alignment provides some advantages compared with other
tracking algorithm. Pairwise alignment tracking makes it possible for one object to have multiple
alignments. It provides the potential of capturing important cell activity such as cell differentiation.
Tracking of multiple predecessors are very important in cell migration study. Current tracking
algorithms, which relying on statistical similarity between alignment objects, show less robustness
to such problem. Algorithm such as KDE mean shift has no implicit possibility of performing
alignment of multiple predecessors pairwise alignment tracking may be a breakthrough in this
aspect.
- 36 -
4.3 Overlap Ratio Based Tracking
Cell accomplishes its motion by elongating and shrinking its lamellipodium and cytoskeleton
structure, which is also known as amoebic motion model. Such elongation guarantees that partial
of cell body is always fixed. Meaning an overlap between objects in subsequential images can
always be found. Overlap ratio tracking algorithm calculates the overlap ratio between two cells in
subsequential frames. Thus, alignments can be made based on highest overlap ratio.
The time resolution is a major fraud in overlap-based tracking. If time resolution is too low,
missing trajectory or early termination will occur since no overlap can be determined. Therefore, a
threshold time resolution must be defined and video must at least provide logical connection.
Meaning the time interval must be small enough to provide correct overlap. In order to provide
error-tolerance in less sensitive situation, a color-histogram based similarity measurement is
introduced.
4.3.1 Overlap Ratio Based Tracking
Overlap ratio can be measured by following formula:
Given two objects in adjacent frames, obj1 and obj2
OverlapRatio(obj1, obj 2) =
count ( pixel ∈ (mask obj1 I mask obj 2 ))
count (mask obj1 )
However, if a small object is a subset of a large object in subsequential frame, wrong trajectory is
frequently made between a small object and a large object. In case of objects with different size,
uncorrected overlap ratio may have the following bias:
Obj1
Large
Small
Obj2
Small
Large
Overlap
Low
High
Therefore, we introduce a reverse-correction of overlap measurement. The idea of reversecorrection is to avoid alignment between large object and small object by requiring overlap
information from both objects. The corrected overlap ratio is defined as the average of overlap
ratio and reversed overlap ratio.
Re verseOverlapRatio(obj 2, obj1) =
OverlapRatioCorrected =
count ( pixel ∈ (mask obj1 I mask obj 2 ))
count (mask obj 2 )
OverlapRatio(obj1, obj 2) + Re verseOverlapRatio(obj 2, obj1)
2
A 100% overlap ratio may occur when the mask of obj1 is a subset of obj2, which may have no
true connection in video. The reverse overlap ratio may only be 10%. The corrected overlap ratio
is 55%. True alignment will only occur when two overlap ratios are high. If overlap ratio is 80% for
obj1 and obj2 while reverse overlap ratio is 70%, the corrected overlap ratio will be 75%. The
corrected overlap ratio significantly reduces the chance of subset-alignment.
4.3.2 Overlap Ratio Based Tracking in Cell Image
The table 9 shows two samples of cell tracking in (a) video with low cell density and high
magnification level (b) video with high cell density and low magnification level. Trajectories of two
samples are approaching to 85% of tracking accuracy according to specialist visual confirmation.
- 37 -
Table 9
First Frame of Video (a)
First Frame of Video (b)
Last Frame of Video (a)
Last Frame of Video (b)
Trajectories of (a)
Trajectories of (b)
- 38 -
- 39 -
4.4 Parametric Motion Model Tracking
The Parametric Motion Model Tracking assumes that motion of object will result in rigid transform,
which can be described by affine transformation. It is nearly true for all rigid objects with no detail
motion behaviors or detail structure changing[2].
The Parametric Motion Model Tracking presumes linear motion models such as affine
transformation. The motion model in this tracking algorithm does not describe the trajectory but
the transformation of objects in 2D projection of 3D world. This approach requires high quality
image with clear boundary of object. Normally, the tracking region is defined by rectangle or circle.
Irregular shape of interesting region is very rare since it causes computation difficulty in affine
transformation. The major usage of this approach is to track solid object with no detail structure
changes, such as car, football, missile, or airplane. It also prefers implicit linear motion model.
If image changes are approximately affine, correspondence for alignment can be achieved by
treating an image as a function of an affine transformation given by parameters a= (a0, a1, a2, a3,
a4, a5) T, where
⎡a ⎤ ⎡ a
f ( x, a ) = ⎢ 0 ⎥ + ⎢ 1
⎣ a3 ⎦ ⎣a 4
a 2 ⎤ ⎡ x − xc ⎤
a 5 ⎥⎦ ⎢⎣ y − y c ⎥⎦
Parametric motion model tracking is proven insensitive to shape changes in cell image. Given the
amoebic motion, cell motion cannot be described by rigid motion model. Cell structure change is
more pronounced than rigid object, which cannot be measured by rigid transform. Moreover, most
parametric motion model presumes either a linear or quadratic motion model, which is also not
guaranteed in cell image. However, it remains useful for cells with less pronounced shape
changes and no self-controlled motility behavior. The blood cell is a typical cell type that can be
tracked by parametric motion model.
Since all four tracking principles have been introduced, we will start to compare the performance
of each algorithm in next section, meaning tracking bias. However, only detail comparison of KDE
Mean Shift tracking and Overlap Ratio tracking will be given. The idea of ruling out parametric
motion model and KDE Pairwise Alignment tracking is due to their limitation and immaturity by
principle.
Estimation of tracking bias in real video contains certain level of bias, which is mainly due to
random cell migration and error in manual observation. In addition, real cells are difficult to control
and do not always correctly react towards treatments. In order to estimate tracking bias
objectively, an artificial object test is employed in our research. In next section, we will give a
detail introduction of this artificial object test.
- 40 -
4.5 Artificial Object Test
The error estimation in tracking algorithm is non-trivial and it is done mostly by specialist manual
observation. It only provides abstract and individual impression of the tracking performance. No
statistical estimation of error can be obtained. In addition, cell experiments may result in certain
level of random bias. In order to estimate the performance of tracking algorithms, a controlled test
is designed.
Artificial object test creates videos of objects, which simulating basic cell migration, given userdefined parameter settings. These settings allow user to create video of artificial cells with nearly
optimal or controllable behaviors. In addition, automatic error estimation can be done since true
trajectories of artificial cells are known in advance. The error estimation of artificial object test is
the decisive factors of our later tracking algorithm selection.
The artificial object test consists of two parts. (1) Artificial cells are generated by the core program
named virtual cell factory. Artificial video will be generated based on self-managed artificial cells.
(2) Trajectories found by tracking algorithms will be compared with true trajectories. During this
chapter, we will give detail explanation of our artificial object test.
4.5.1 Virtual Cell Factory
Our version of artificial object test, named virtual cell factory, is designed to simulate two major
behaviors of cell migration, position shifts and shape changes. Parameter settings are defined to
create a user-controlled experiment condition. Virtual cells will be created given initial parameter
settings. After being created, virtual cells will be self-managed based on their internal behavior
simulators, thus, movie can be made from these artificial cells. (see figure 21)
The following parameter settings are defined in our version of implementation.
1. Travel distance (maximum and minimum)
2. X-axis extension (maximum and minimum)
3. Y-axis extension (maximum and minimum)
4. Initial size (x-axis length and y-axis length)
5. Number of Cells
6. Uniformly distributing cells
7. Perfect Connectivity
- 41 -
Figure 21
Travel Distance
Travel distance controls how quick artificial cell shifts position. It is used to simulate the motility
ability of cell. Since the relation between lamellipodium extension and cell velocity is unknown,
the settings of velocity and shape are given separately.
X-axis Extension and Y-axis Extension
X-axis extension and Y-axis extension control how artificial cell extends in both directions. It
simulates lamellipodium extension of cells.
Initial Size
The initial size of artificial cell is used to simulate microscopy magnification level. High
magnification level often results in larger cell size.
Number of Cells
Studying video with few cells has no statistical significance. Therefore, the artificial cell test must
be able to create video with high cell density. User is allowed to create 200 cells at maximum but
the initial size must be properly chosen. Extremely high cell density will disable motility behavior
due to programming limitation.
Uniformly Distributing Cells
This setting allows user to prevent initial masking error. All cells are evenly located in the image
with maximum distance and nearly zero overlapping. The chance of connected cells is nearly
- 42 -
impossible expect extremely high density of cells. If not set to true, cells are randomly located in
image. Yet, in either even or random distribution, cell will not touch the border of screen.
Prefect Connectivity
Uniformly distributing cells may not prevent cells from later overlapping. After long-term moving, it
is highly possible that cell trajectory will cross each other. Cells are completely isolated from
others during entire video. However, similar to uniformly distributing cells option, this option
cannot work with extremely high density of cells. Extremely high cell density will prohibit motility
ability.
Due to programming complicity, rotation and extension by angle is not allowed in current setting.
Amoebic motion model does suggest that cell motion is completed by extending lamellipodium.
However, it does not explicitly suggest that the motion direction is on the same axis of
lamellipodium extension. Still, artificial cell test is designed to test tracking algorithm in an optimal
situation. Introducing further complicity may not be necessary.
The artificial cells created by the artificial cell factory have the following traits:
1. Artificial cell is given a numeric range for its elongation behavior
2. Elongation of artificial cell will be only in one direction
3. Velocity of artificial cell has a predefined linear correlation
4. No other cell activity will appear, such as division
5. Artificial cells will not suddenly appear or disappear, they remains in the screen
6. Cells will not crash into each other and will keep distance from each other
7. Cells can be evenly distributed in the image
Artificial cell test allows us to create controlled cell behaviors to determine which factors will result
in tracking bias. The images show two visualized trajectories of performing tracking algorithm on
artificial videos (See Figure 22 and Figure 23).
49 artificial cells
9 artificial cells
49 artificial cells trajectory
9 artificial cells trajectory
- 43 -
Figure 22
Figure 23
In later experiment, we will use virtual cell factory to create iterative experiment with designed
sets parameter settings. The idea is to find out, numerically, the critical factors that directly
resulting in tracking bias. We will introduce the error estimation method in next section.
- 44 -
4.5.2 Error Estimation in Artificial Object Test
To compare the performance of each tracking algorithm, we must first define the error estimation
for tracking biases. A mean-error based score system is designed in our experiments. The score
is calculated based on the average percentage of correct tracked cells in each path times the
ratio between the number of true paths S t → counting (true) and the number of founded
paths S f → counting ( founded ) . The path ratio is given:
PathRatio =
Sf
St
The error of each tracked path (Tracked Percentage) is given:
Tracked % =
N ( tracked cell in true path )
N ( cell in true path )
The final score of tracking performance is given:
score = PathRatio ⋅ Tracked %
To obtain a high score, the tracking algorithm must return trajectories that fulfill the following
conditions:
• Finding all possible trajectories (path) of objects
• Finding all correct alignments for each paths
Calculation Example:
We have a video has a true path of 20 objects and length of 30 images. The tracking algorithm is
performed on this video. For each tracked path, the average error rate is 15%. Therefore, the
correct tracked percentage is 85%. Only 18 true paths are found by tracking algorithm. Thus, the
final score is calculated:
score = 85% ⋅
18
= 76.5%
20
The following table is a sample of score for each tracking algorithm. In section 4.5.3, a serial of
iterative perturbation experiments will be conducted. We expect these experiments to providing
not only the robustness estimation for each tracking algorithm, but also finding the decisive
factors of successful tracking.
Tracking Algorithm
KDE Mean Shift Tracking
Overlap Tracking
Pairwise Tracking
Standard Mean Error of 160 tests
85%
76%
52%
When artificial object videos are created, tracking is performed independently on videos. The
trajectories of tracking algorithms will be compared with the true trajectories recorded in advance.
With the help of artificial object test, we can not only compare the tracking performance of each
algorithm, but also perform a deep study of the factors that may cause tracking biases. In next
section, we will demonstrate a perturbation study of tracking algorithm using artificial object test.
- 45 -
4.5.3 Perturbation Study in Artificial Object Test
16 groups with 10 iterative tests for each are used to estimate the performance of three
algorithms. The purpose of repeated test is to provide objective error estimation. The settings for
groups are manually set and paired in order to perform the later perturbation test. 20 tests will be
given to one set of perturbation setting (see Appendix III-1).
The mean error of 160 tests is given as follow. For detail data, please referring to datasheet.
Tracking Algorithm
Standard Mean Error
KDE Mean Shift Tracking
85%
Overlap Tracking
76%
Pairwise Tracking
52%
The overall perform suggests that KDE mean shift tracking has the highest correct score while
pairwise tracking is the lowest. KDE mean shift tracking and overlap tracking have similar
performance but detail test suggests that overlap tracking relies on optimal situation while the
KDE mean shift tracking does not. Moreover, due to the immaturity of pairwise tracking algorithm,
it will not be involved in perturbation analysis section.
As discussed in the very beginning of this report, several factors will causes the tracking
algorithm to fail. These factors are involved in the following perturbation analysis. Several factors
are believed causing more error in tracking algorithm. These factors include:
1. motion speed
2. cell density
3. cell remote
2
4. cell separation
Experiment protocols are given as following:
• 20 tests are equally divided into two groups. One group is given cell speed from 5 to 10
pixels. Another is given cell speed from 10 to 20 pixels.
• 20 tests are equally divided into two groups. One group is given 5 cells in total. Another is
given 20 cells in total.
• 20 tests are equally divided into two groups. Cells in one group are randomly located.
Cells in another group are located remotely.
• 20 tests are equally divided into two groups. Cells in one group are randomly separated.
Cells in another group are located separately.
Each artificial cell is given following common setting.
• Cells are given the same initial size (width 45 and height 30)
• The total length of video is identical (30 frames)
2
The cell separation is not the cell differentiation by biological definition. We define the cell
separation as the masking condition of cell connectivity. It refers to the overlap ratio between two
non-alignment objects.
- 46 -
Two parameter
settings subjected
to target factor
10 iterative
tracking based on
parameter setting
1
10 iterative
tracking based on
parameter setting
2
Error estimation
based on true
trajectories
Error estimation
based on true
trajectories
Does tracking
error significantly
increase?
false
This factor is not
decisive
true
This factor is
decisive
Figure 24
For detail of tracking error, please refer Appendix III-1: Artificial Object Test and Perturbation Test.
Based on tracking results, we determine the following decisive factors that affecting tracking
accuracy.
Motion Speed
High motion speed causes incorrect tracking. The differentiation ratio for both tracking algorithms
is around 6% (see Table 10). Overlap tracking is less affected by motion speed, yet global
performance of overlap tracking is still lower than KDE mean shift. Time resolution is a
manifestation of motion speed setting. An extremely high motion speed may also be the result of
a low time resolution. The problem of time resolution will be introduced in section 4.7 Time
Resolution Study.
Table 10
t(5/10)
t(10/20) Diff %
KDE
88.45
82.685
6.5178
Overlap 78.232
73.784
5.6857
Cell Density
The cell density has a stable effect on tracking algorithm. Both tracking algorithms show that error
is increased by 13% (see Table 11). More cells means higher chance of identity and worse
connectivity condition, which causes both density based algorithm and non-density based
algorithm to fail. However, the differentiation ratio suggests that overlap tracking is less affected
by cell density. In addition, compared with motion speed, cell density has more affect on tracking.
- 47 -
Table 11
KDE
Overlap
c(5)
91.06
81.394
c(20)
78.632
70.465
Diff %
13.6481
13.4273
Cell Remote
The initial separation of cell has a critical role in KDE mean shift algorithm. Accuracy of KDE
mean shift will significantly drop due to bad initial cell location (See Table 12). The overlap
tracking is less affected by the cell remote. Therefore, a correct initial mask is critical for a
successful tracking.
Table 12
r(true)
r(false) Diff %
KDE
91.936
76.718
16.55282
Overlap 80.556
71.461
11.29028
Cell Separation
The later separation of cells during entire motion is a critical factor for overlap ratio tracking.
Compared with initial separation (cell remote), both tracking algorithms are less affected by later
separation (see Table 13). Yet, our further research on the combination of cell remote and cell
separation reveals that the collaborated affect of both conditions will be a decisive factor (see
Table 14).
Table 13
S(true)
s(false) Diff %
KDE
87.874
82.031
6.649293
Overlap 79.831
72.186
9.57648
Combination of Cell Remote and Cell Separation
Combination study of both cell remote and cell separation suggests that in optimal situation while
cells are initially separated and later not connected. KDE mean shift tracking will provides
accuracy around 96% and overlap ratio tracking will provides accuracy around 83% (see Table
14).
Table 14
r(t) 3 &s(t) 4 r(t)&s(f) r(f)&s(f) r(f)&s(t)
KDE
96.449
89.082
74.217
79.637
Overlap 83.598
77.515
66.856
76.066
Based on the general performance and perturbation test, we believe that KDE Mean Shift
tracking is indeed a robust tracking algorithm in cell tracking. Our next step is to define several
morphology and motility estimations for cell migration study. The definition of feature
measurements will be discussed in next section.
3
4
Abbr. cell remote is set to true
Abbr. cell separation is set to true
- 48 -
4.6 Measurements in Cell Migration Video
With the help of tracking algorithm, we can now perform statistic measurement on cell migration.
Definitions of several statistic measurements (features) will be introduced in this section. These
features are employed to preserve motility and morphology changes of cell migration. The current
features used in our research are divided into two categories:
1. Static features
2. Continuous features
Static features are the statistic measurements being extracted from each cell in each static frame.
These features include area, perimeter, coordination of mass center, etc. The reason of such
naming is due to their limitations of representing correlation between subsequential frames.
Information related to motility or continuous morphology cannot be directly observed from static
features, instead continuous features are introduced.
Continuous features are the statistic changes extracted from static features. These features
include velocity, motion linearity, etc. They are specially designed for preserving motility and
morphology changes during cell migration.
Since both static feature and continuous feature provide valuable information on cell migration,
data analysis will be performed on both raw and generalized features.
4.6.1 Static Features
The Static Features are the measurements result extract directly from static frames. The
measurements are independent and no correlations are set for subsequential frames. The Static
Features defined include:
Raw Feature
Description
Area
The size of cell
Perimeter
The perimeter of cell
Center of Mass
The mass center of the binary mask of cell
Extension
Extension measures how much the shape differs from the circle
Dispersion
Dispersion is the minimum extension that can be attained by uniform
compression of the shape
Elongation
Elongation measures how much the shape must be compressed along its
long axis in order to minimize the extension
Circle
It is the ratio of the area of the shape to the area of a circle (the most
Compactness
compact shape) having the same perimeter
The mean of Static Features preserve global information such as the average size of cells, or the
average elongation of cells. However, the static features do not preserve information about
diversification of cell migration.
Area
The area of cell is measured based on the total number of pixels in the binary mask of cell. The
area measurement of cell is identical to the size of cell. Frequently, area is also referred as zero
order moment, which is given by
sy −1 sx −1
μ pq = ∑∑ (x − x ) p ( y − y )q F (x, y ) ,where p=0 and q=0
0
0
By simplification, the area of binary mask is given by μ 00
sy −1 sx −1
= ∑∑ F (x, y )
0
- 49 -
0
The F(x,y) is binary function that F(x,y)=1 if pixel at (x,y) belongs to object or F(x,y)=0. x is the x
coordinate of mass center and y is the y coordinate of mass center. The F(x,y) can also be the
intensity function, which will result in a intensity-weighted calculation. It may be useful when
intensity value of cell contain a biological translation.
Perimeter
Perimeter of cell is calculated from the length of edge of the binary mask. With different choice of
connectivity, the result of perimeter may different from each other. Normally, perimeter obtained
from 4-connected is longer than from 8-connected. In our version, we will use 8-connected.
Center of Mass
By the definition of computer vision, the center of mass is the position given by:
sy −1 sx −1
μ pq = ∑∑ (x − x ) p ( y − y )q F (x, y )
0
0
μ10
μ 00
μ
MC _ Y = 01
μ 00
The x-coordinate of mass center is MC _ X =
F(x, y) = b (x, y) binary mask
The y-coordinate of mass center is
F(x, y) = b (x, y) binary mask
We here assume that the mass center in binary mask is overlapping with the physical mass
center of cells. Thus, the center of mass will be used as the reference point for distance and
motion measurements.
Extension
Extension is a moment invariant measurement feature. Extension measures how much the shape
differs from the circle. The calculation of extension is based on normalized moment of mask.
sy −1 sx −1
μ pq = ∑∑ (x − x ) p ( y − y )q F (x, y )
0
0
Normalized Moment
η pq =
μ pq
p+q
,r =
+1
r
2
μ 00
φ1 = η 02 + η 20 , φ 2 = (η 02 − η 20 )2 + 4η112
(
)
(
λ1 = 2π φ1 + φ 2 , λ 2 = 2π φ1 − φ 2
)
The extension is defined by E = ln λ1
The major advantage of moment invariant is its native scale-free behavior. Regardless of how
microscopy magnification will change, the moment invariant remains consistent unless affected
by the mask detail revealed by higher resolution.
Dispersion
Dispersion is a moment invariant measurement feature. Dispersion is the minimum extension that
can be attained by uniform compression of the shape. It is invariant to stretching, compressing or
shearing the shape in any direction. The calculation of dispersion is based on normalized moment
of mask.
- 50 -
sy −1 sx −1
μ pq = ∑∑ (x − x ) p ( y − y )q F (x, y )
0
0
Normalized Moment
η pq =
μ pq
p+q
,r =
+1
r
2
μ 00
φ1 = η 02 + η 20 , φ 2 = (η 02 − η 20 )2 + 4η112
(
)
(
λ1 = 2π φ1 + φ 2 , λ 2 = 2π φ1 − φ 2
The extension is defined by D =
)
1
ln λ1λ2
2
The major advantage of moment invariant is its native scale-free behavior. Regardless of how
microscopy magnification will change, the moment invariant remains consistent unless affected
by the mask detail revealed by higher resolution.
Elongation
Elongation is a moment invariant measurement feature. Elongation measures how much the
shape must be compressed along its long axis in order to minimize the extension. The calculation
of elongation is based on normalized moment of mask.
sy −1 sx −1
μ pq = ∑∑ (x − x ) p ( y − y )q F (x, y )
0
0
Normalized Moment
η pq =
μ pq
p+q
,r =
+1
r
2
μ 00
φ1 = η 02 + η 20 , φ 2 = (η 02 − η 20 )2 + 4η112
(
)
(
λ1 = 2π φ1 + φ 2 , λ 2 = 2π φ1 − φ 2
The extension is defined by L =
)
1 λ1
ln
2 λ2
The major advantage of moment invariant is its native scale-free behavior. Regardless of how
microscopy magnification will change, the moment invariant remains consistent. The figure 25
shows samples of extension, dispersion, and elongation measurements for different typical
shapes of cells.
- 51 -
Figure 25
Circle Compactness
The circle compactness or shape factor[16] is a common shape measure. It is the ratio of the area
of the shape to the area of a circle (the most compact shape) having the same perimeter. Circle
compactness is calculated given formula compactness =
Perimeter 2
.For a circle, the ratio is
4πArea
one; for a square, it is π / 4; for an infinitely long and narrow shape, it is zero (see figure 26).
Figure 26
4.6.2 Continuous Features
The continuous features are features preserving trends of changing or correlation. The current
continuous features include several important change rates and motion speed (see table 15).
- 52 -
Table 15
Raw Feature
Velocity
Absolute Position Shift
Trajectory Length
Motion Linearity
Description
The motion speed of cells
The position shift between the first frame and last frame
The total travel distance based on the mass center of cell in all frames
It measures the linearity of motion trajectory
Velocity
By physical definition, velocity is defined as the rate of change of position. It is given by v =
ΔS
.
Δt
The position shifting ΔS is normally calculated by Euclidian distance between two mass centers
(reference point) of the same object in subsequential frames. The velocity of cell in figure 24 is
calculated as v AB =
2131.32
2026.89
, v BC =
.
30
30
Absolute Position Shift
The absolution position shift of cells is used to describe the position difference between first
frame and last frame. Such position difference is calculated based on the mass center of each
cell in the first frame and last frame. The figure 23 shows a sample of U-shape motion trajectory
of one cell. The absolution position shift is 1751.55 in figure 27.
Trajectory Length
The trajectory length is defined as the total travel distance of cells. The distance of travel is
calculated based on the mass center of cells in subsequential frames. The trajectory length of cell
in figure 19 is S = 2131.32 + 2026.89 .
Motion Linearity
Motion linearity is widely used but not officially defined. Our definition of motion linearity is given
as the ratio between absolute position shift and trajectory length:
Linearity =
S AC
, with given motion trajectory (see figure 28)
S AB + S BC
For long-term running, global linearity may lose important detail or periodic trends. Special case
such as U-type or O-type motion may result in low linearity, which is in fact considered linear.
Therefore, for long-term measurement, window based calculation of motion linearity is employed
in our version of implementation.
We use the U-type motion shown above to explain how window based calculation is performed
(see figure 27 and table 16).
Table 16 Artificial Distance Measurement
SEQ# Distance
1
0
2
1233.74
3
1010.58
4
1048.50
5
1032.42
Window based calculation of global motion linearity with a window-size of three:
Linearity13 =
2113.70
= 0.9418
1233.74 + 1010.58
- 53 -
1964.83
= 0.9442
1048.5 + 1032.42
Linearity13 + Linearity 35
Linearity global =
= 0.943
2
Linearity 35 =
Standard calculation of global motion linearity:
Linearity global =
1751.55
= 0.4049
1233.74 + 1010.58 + 1048.50 + 1032.42
A standard calculation of global motion linearity will be significantly smaller than window-based
calculation of global motion linearity. Based on the discussion with biological specialists, a
window-based calculation is preferred. It is believed that window-based calculation gives more
meaningful measurement than non-window-based calculation.
1751.55
5
2
1032.42
19
64
.83
1233.74
3
211
.70
1
4
8
.5
10
10
3
50
8.
4
10
Figure 27
C
B
9
6.8
202
A
2
2131.3
Figure 28
- 54 -
These features have been frequently employed in both computer vision and microscopy
techniques[22][23]. Yet, study of these features in a time-lapse and continuous style remains atelic.
Feature measurements in time-lapse video result in a time-based signal (see figure 29) while
summarization approaches such as mean and standard deviation are used without considering
the statistic representation of these approaches in frequency domain. In next section, we will give
the biological translation of summarization approaches including mean, standard deviation, and
frequency using tracking analysis results from MTLn3 cell.
elongation
1.2
1
0.8
0.6
0.4
0.2
0
1
8
15 22 29 36
43 50 57 64 71
78 85 92 99 106 113 120
Figure 29 the elongation measurement of one cell in a 120-frames video
4.6.3 Summarization Approaches
The mean and standard deviation are popular summarization approaches in statistics. However,
both of them may only preserve partial trends of dataset. To compensate the shortage of mean
and standard deviation in time-based problem, the frequency will be employed. We will first
introduce the difference between these three approaches.
By mathematic definitions:
• The mean describes the central location of the data
• The standard deviation describes the spread of data
• The frequency measures the number of occurrences of a repeating event per unit time
By biological translation:
• A mean (or average) of area for one cell in a time-lapse video gives the true size of this
cell
• A standard deviation of area for one cell in a time-lapse video gives how significant it
changes its size
- 55 -
•
A frequency of area for one cell in a time-lapse video gives how quickly it changes its
size 5
To explain the difference between mean, standard deviation, and frequency, we introduce a
family of sine function as an optimal result from tracking analysis. The figure 30 and table 17 give
a brief idea of how these three approaches responding to this family of sine function. This
experiment shows to which factors that each of these approaches responding.
Mean, Std, and Frequency
5
y=sin(2*pi.*t)
y=2.*sin(2*pi.*t)
y=sin(4*pi.*t)
y=sin(2*pi.*t)+4
4
3
2
1
0
-1
-2
1
1.5
2
2.5
3
t (time)
3.5
4
4.5
5
Figure 30
Table 17
Function
y = sin (2πT )
y = 2 sin (2πT )
y = sin (4πT )
y = sin (2πT ) + 4
Description
A measurement result from a optimal control
group
A measurement result from a treated group with
larger lamellipodium extension changes
A measurement result from a treated group with
high frequency of lamellipodium extension
changes
A measurement result from a treated group with
longer lamellipodium in general
mean
0
std
0.707
frequency
2
0
1.414
2
0
0.707
4
4
0.707
2
The function y = sin (2πT ) is the control signal while functions y = 2 sin (2πT ) , y = sin (4πT ) ,
and y = sin (2πT ) + 4 are the treated signals.
•
•
The mean only responds to the value-shift (see y = sin (2πT ) + 4 )
The standard deviation only responds to the amplification level (see y = 2 sin (2πT ) )
5
By our definition, a continuous or gradually shape changing such as shrinking or extending in
amoebic motion will not be considered as significant shape change. However, a complete shrinkextend motion circle will be considered as significant shape change.
- 56 -
•
The frequency only responding to frequency change (see y = sin (4πT ) )
It is notable that all three summarization approaches preserving important but limited information
on cell migration. Thus, all three of them will be consulted in our research. Figure 31 shows the
procedure of converting raw feature measurement into dataset, which can be processed by data
analysis.
Feature
measurement for
cell 1 in frame 1
Feature
measurement for
cell 1 in frame n
Feature
measurement for
cell n in frame 1
Feature
measurement for
cell n in frame n
Mean, Standard
Deviation, and
Frequency
summarization for
cell 1 for each
feature
Summarized
dataset for large
volume data
analysis
Mean, Standard
Deviation, and
Frequency
summarization for
cell n for each
feature
Figure 31
With the combination of tracking algorithm and feature measurements, we can now convert raw
video into objective and meaningful numeric estimation of cell migration. In large volume data
analysis, measurement of single object in each frame will not be considered. Instead, these
features will be summarized by mean, standard deviation, and frequency for each object. Yet,
due to masking bias and time resolution problem, measurements will always contain a certain
level of bias. Even with extremely high time resolution, small oscillations in numeric estimation
are difficult to avoid. In addition, error in masking will result in exceptions in feature measurement.
These measurement errors will especially affect the estimation of frequency. However, it can be
fixed by implementing noise suppression using moving average filtering.
4.6.4 Noise Suppression in Measurement
Bias or noise in feature measurement is difficult to avoid but fixable. By considering time-lapse
measurement as a signal, these oscillation (or noise) can be suppressed by mean filtering. The
idea of mean filtering is to correct the current value based on its neighbor. The algorithm of mean
filtering in discrete domain is defined as[16]:
xi =
k
1
∑ xi + j where xi is a measurement point at time i, 2k+1 is the kernel size
2k + 1 j = − k
- 57 -
Plot of y1 and y2
1.5
y1=sin(2*pi*t)
y2=sin(2*pi*t)+0.5*rand-0.25
1
0.5
0
-0.5
-1
-1.5
0
50
100
150
200
250
Figure 32 a noise suppression using moving average filter
Plot of y1 and y3
1
y1=sin(2*pi*t)
y3=conv2(y2,meanfilter)
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
50
100
150
200
250
Figure 33
Figure 32 shows a sample of noise suppression using moving average. An artificial noisy sinesignal y2 (blue line in figure 32) is created from adding a normally distributed noise ranging from 0.25 to 0.25 to a normal sine-signal y1 (red line figure 32). The sine-signal y3 (blue line in figure
- 58 -
33) is the result of noise suppressed y2 using moving average. Compared with noisy signal y2,
suppressed signal y3 is more likely to preserve the original trends of y1. This experiment is
repeated by 1,000 times and the average of the summarization estimation of three signals is
given in table 18.
Table 18 signal sample rate 1000 Hz
mean Std
freq
Y1 0.000 0.707 0.020
Y2 0.000 0.722 0.639
Y3 0.001 0.695 0.114
Table 18 shows that a normally distributed noise will not have significant effect on mean and
standard deviation estimation. However, it significantly increased the frequency estimation. Such
error can be significantly reduced by mean filtering.
However, we also notice that noise suppression is significantly effected by sample rate of signal.
It will only work properly when signal is oversampled. An undersampled signal will cause a
notable error in mean and standard deviation estimation after noise suppression (see table 19
and table 20).
Table 19 signal sample rate 200 Hz
Mean std
freq
Y1 0.000 0.707 0.098
Y2 -0.004 0.712 0.341
Y3 -0.005 0.680 0.098
Table 20 signal sample rate 20 Hz
Mean std
freq
Y1 0.000 0.000 0.600
Y2 -0.116 0.144 0.200
Y3 -0.082 0.036 0.400
It is obvious that when sample rate is reducing (table 18 at 1000Hz, table 19 at 200Hz and table
20 at 20Hz), the three estimations of original signal y1 start to gain more bias. The time-resolution
(or sample rate) problem will be further discussed in section 4.7.
- 59 -
4.7 Measurement Bias in Time Resolution
The definition of time resolution may vary from different aspects. We define time resolution as the
number of frame taken per second for each well in plate. In this section, we focus on
measurement bias caused by time resolution.
Figure 34
Time resolution is a critical imaging condition to time-lapse tracking problem. It is a special type of
sampling problem with similar trait of sound sampling (see figure 34). In order to study the affect
of sampling rate, a time-lapse video of MTLn3 cell with extreme-high time resolution is created for
subsampling. In this section, we will introduce an experiment of estimating measurement bias
using the same video with different time resolution.
- 60 -
4.7.1 Time Resolution Experiment
To study bias of feature measurements with different time resolution, a video of extreme-high
time resolution is created. The ultimate video (mtln3-gfp-paxillin+hgf+egf_30) is a sequence of
120 images with approximate 15 cells taken at 30 seconds per frame per well. Subsampling is
employed to create data with artificially reduced sample rate. The idea of subsampling is to
extract a subset of frames from original sequence. Our rule of extracting subset is described in
figure 35:
false
j=j+1
1
2
3
Divided original
frames into n subgroups with m
frames
4
1 2 3 4
1
2
3
4
Extract one frame
from sub-group i
at position j
1 2 3 4
1
2
3
False
i=i+1
4
All sub-groups
are visited
1 2 3 4
1
2
3
Bias estimation
4
true
All positions in
sub-group are
visited
true
Mean bias
estimations is
calculated
1 2 3 4
Tracking is
performed on
subsampled video
Figure 35
Pseudo Code: For time-resolution h times lower than original
Define a window size h for all frames m
For each position j in subgroup
For each subgroup i defined by window size h
Extract frame at position j
A subset of original frames is created
Bias of measurements is estimated between full frames and subset of frames
Average of bias for each position j in subgroup is calculated
In our implementation, weighted correlation, which is a score system combining correlation
coefficient and mean square error, is used as bias estimation.
- 61 -
4.7.2 Weighted Correlation
In our implementation, the weighted correlation is defined as the ratio between correlation
coefficient and mean square error:
wc =
Correlation
MSE
The best-subsampled dataset should have high correlation and low mean square error compared
with original data, meaning preserve maximum global trends and minimum numeric difference.
Since the numeric scalar of correlation coefficient and MSE are different, both of them must be
normalized to range from zero to one given:
Value N =
Value − Min
Max − Min
One exception in weighted correlation (WC) calculation is that when MSE is equal to zero, the
WC score will be set to one since zero MSE means measurements from subsampled data has no
numeric difference from original data.
The advantage of using weighted correlation to avoid the limitation of either correlation coefficient
or mean square error. The table 21 explains the essential of combining both criteria (see also
Figure 36~39).
Summary of Correlation and MSE combination
Table 21
Correlation MSE Description
High
High Subsampled data tracking result may have same trends in changing but
values are highly different from non-subsampled data. Considering only
correlation may mislead conclusion.
High
Low Subsampled data tracking result is nearly identical to non-subsampled. This
is the best situation since subsampled data consists of low bias. It is the
major reason why both correlation and MSE are employed.
Low
High This is the worst case. Neither the subsampled data has same trends in
changing, nor is value even close. It may occur when the subsampled data
completely lost its original features. Such situation should be avoided.
Low
Low This is a special case. Subsampled tracking result may not have same
trends in changing, but its value is quite close to non-subsampled data. It
happens if non-subsampled tracking result is close to global mean of
subsampled data. However, such situation may occur rarely since unlikely
low MSE will lead into low correlation.
- 62 -
Correlation/MSE Criteria
9.00
8.00
7.00
6.00
5.00
4.00
3.00
2.00
1.00
0.00
-1.00
-2.00
1
2
3
4
5
Non-Sub
1.00
2.00
3.00
4.00
5.00
Sub1
1.4
1.6
3.4
3.6
5.4
Sub2
1.8
1.2
3.8
3.2
5.8
Sub3
2.2
0.8
4.2
2.8
6.2
Sub4
2.6
0.4
4.6
2.4
6.6
Sub5
3
0
5
2
7
Sub6
3.4
-0.4
5.4
1.6
7.4
Sub7
3.8
-0.8
5.8
1.2
7.8
Figure 36
Correlation/MSE Compare
1.2
1
0.8
0.6
0.4
0.2
0
wc normal
Sub1
1
Sub2
Sub3
Sub4
Sub5
0.219325894 0.079730253 0.034069449 0.014730053
Sub6
Sub7
0.00520286
0
Figure 37
Figure 37 shows the weighted correlation corresponding to global data trends.
- 63 -
MSE (Shifted Test)
7.00
6.00
5.00
4.00
3.00
2.00
1.00
0.00
-1.00
Non-Sub
1
2
3
4
5
1.00
2.00
3.00
4.00
5.00
Sub1
1.4
2.4
3.4
4.4
5.4
Sub2
1.80
2.80
3.80
4.80
5.80
Sub3
2.2
3.2
4.2
5.2
6.2
Sub4
2.60
3.60
4.60
5.60
6.60
Sub5
0.6
1.6
2.6
3.6
4.6
Sub6
0.20
1.20
2.20
3.20
4.20
Sub7
-0.2
0.8
1.8
2.8
3.8
Sub8
-0.60
0.40
1.40
2.40
3.40
Figure 38
wc normal
1.2
1
0.8
0.6
0.4
0.2
0
wc normal
1
2
1
0.2
3
4
0.05185185 2.3685E-17
5
6
7
8
1
0.2
0.05185185
0
Figure 39
Figure 39 shows two peaks (sub1 and sub5) with the shortest shift distance from non-sub data.
- 64 -
4.7.3 Sampling Rate
The weighted correlation of dataset (mtln3-gfp-paxillin+hgf+egf_30) is given in figure 40.
Correlation/MSE Estimation (0-120)
1.2
1
Correlation/MSE
0.8
0.6
0.4
0.2
0
rate8
rate10
rate12
rate14
rate16
test1
1
0.86842984
0.743869814
0.003690154
0
test2
1
0.680856654
0.577925956
0.063385182
0
test3
0.939061381
1
0.478411354
0.508095131
0
test4
0.96650698
1
0.508286344
0.076730249
0
Subsample Rate
Figure 40
The diagram (figure 41) is the average weighted correlation from four videos. The rate (number)
means the subsample rate. A rate8 curve is the weighted correlation measurement with
subsample rate 1-to-8, meaning every one frame will be extracted from each eight frames.
Observable from the curve, the weighted correlation will reduce while decreasing subsample rate,
meaning bias will increase while time resolution reduces. It also implicitly suggests a linear
correlation between subsample rate and weighted correlation. Therefore, the optimal threshold
subsample rate will always be the highest subsample rate. The current data suggests 1/8~1/12
subsampling will be an optimal range of subsample rate.
The diagram (figure 41) also consists of a 0-59 curve. This is the curve of bias estimation for a
video with half-length of original video (60 frames in total). The same subsampling experiment is
performed on the half-length video. As you can see, bias is significantly increased with high time
resolution while insignificant increased with low time resolution. Our conclusion is that when time
resolution drops below a threshold, the measurement of video data will be significantly
undersampled. Such undersampling causes high error in reconstruction of original data, which
also resulting in high bias.
- 65 -
Mean Correlation/MSE
1.2
1
Mean Correl/MSE
0.8
0.6
0.4
0.2
0
rate8
rate10
rate12
rate14
0-120
0.97639209
0.887321624
0.577123367
0.162975179
rate16
0
0-59
0.739105235
0.708234578
0.597895185
0.144162856
0.023665797
Subsample Rate
Figure 41
In this section, we describe our experiment on estimate measurement bias due to time resolution.
We demonstrate the reason why 6 min is chosen as the maximum time resolution for MTLn3 cells.
Although it shows that 6 min is a clear threshold for unbias measurement, the threshold is celloriented. Meaning if a desired threshold needs being found in cell lines, the same experiment
must always be performed.
In chapter 4, we introduced several tracking algorithms and our implementation of feature
measurements. We also demonstrate how time resolution affects features measurements (bias
estimation). In next chapter, we will demonstrate a real application of our tracking analysis and
high-throughput data analysis with siRNA screening of MTLn3.
- 66 -
5 MTLn3 Cell Study
In this chapter, we will first introduce and test the credibility of our high-throughput data analysis
with real video of MTLn3 with three well-understood experiments. Then, the same analysis will be
performed on video data of experimental siRNA screening without preknowledge of
corresponding reaction, meaning double blind test.
5.1 High-Throughput Data Analysis
The idea of high-throughput data analysis is to employ computational technologies to accelerate
or fully automate the processing, quantification and analysis of large amounts of high-informationcontent biomedical imagery. Modern image analysis systems augment an observer's ability to
make measurements from a large or complex set of images, by improving accuracy, objectivity,
or speed. A well-developed analysis system may completely replace the specialist observation.
The study of cell migration with tracking algorithm is classified under Phenotypic Research.
A straightforward method to study behavior change is the differential ratio. It is calculated based
on increasing or decreasing percentage in each feature between control samples and treated
samples. The formula can be generalized as:
diff =
f (treat ) − f (ctrl )
Where f (treat) is a certain feature measurement of treated group and
f (ctrl )
f (ctrl) is a certain feature measurement of control group
However, differential ratio relying on global average of each feature only shows the numeric
increasing or decreasing between control and treated samples. Since global average may not
preserve the true trend, two other computational technologies being employed in our research.
These two technologies include feature selection and statistical test. After a brief explanation, we
will demonstrate the implementation of these two techniques and the reason why both of them
must be consulted.
5.1.1 Feature Selection
Feature selection is the technique of selecting a subset of relevant features for building robust
learning models. When applied in biology domain, the technique is also called discriminative gene
selection, which detects influential genes based on DNA microarray experiments. We will
implement the same principle in cell behavior discrimination. We believe the behavior
measurement providing largest variety is the significantly changed behavior between groups (see
figure 42).
plGFP
mean velocity
EGFEGF+
4
2
0
-2
0
-2
montion Linearity
-4
-1
0
1
2
3
freq extension
Figure 42
- 67 -
The most popular form of feature selection is stepwise regression. It is a greedy algorithm that
adds the best feature (or deletes the worst feature) at each round. The main control issue is
deciding when to stop the algorithm. In machine learning, this is typically done by cross validation.
In statistics, some criteria are optimized. This leads to the inherent problem of nesting. Robust
methods have been explored, such as “Branch & Bound” and “Piecewise Linear Networks”.
Forward Selection
The repeating of feature selection is to find the stable features. Forward searching with crossvalidation may result in random bias since training set is randomly picked by cross-validation
procedure. Features may swap position in ranking even these features are the true influent
features. The order of features found by feature selection does not possess important information
in our study. Therefore, we need to find the features, which most frequently appearing in feature
selection results. For detail of forward selection algorithm, please refer to [19]
Backward Selection
Backward selection is more stable compared with forward selection. The idea of backward
selection is to remove the feature that having the least impact on criterion, meaning the feature
has no significant improvement in variety in our case. Unlike forward selection, which intending to
find combination of most influent features, backward selection intends to remove redundant
features. This is also the reason why backward selection is much more stable. For detail of
backward selection algorithm, please refer to [19]
5.1.2 Statistical Hypothesis Test
Statistical hypothesis test or significant test is a method of making statistical decisions from and
about experimental data. Null-hypothesis testing just answers the question of "how well the
findings fit the possibility that chance factors alone might be responsible"[21]. This is done by
asking and answering a hypothetical question. One use is deciding whether experimental results
contain enough information to cast doubt on conventional wisdom. Two popular approaches are
the Z-Test and T-Test.
Z-Test
The Z-test is a statistical test used in inference which determines if the difference between a
sample mean and the population mean is large enough to be statistically significant, that is, if it is
unlikely to have occurred by chance. The Z-test is used primarily with standardized testing to
determine if the test scores of a particular sample of test takers are within or outside of the
standard performance of test takers.
The Z-test requires the following to be known:
•
•
•
•
σ
μ
(the standard deviation of the population)
(the mean of the population)
x (the mean of the sample)
n (the size of the sample)
In actuality, knowing the true σ of a population is unrealistic except for cases such as
standardized testing in which the entire population is known. In cases where it is impossible to
measure every member of a population, it is more realistic to use a t-test, which uses the
standard error obtained from the sample along with the t-distribution.
T-Test
A t-test is any statistical hypothesis test in which the test statistic has a Student's t distribution if
the null hypothesis is true. It is applied when sample sizes are small enough that using an
assumption of normality and the associated z-test leads to incorrect inference. Among the most
frequently used t tests are:
- 68 -
•
A test of the null hypothesis is that the means of two normally distributed populations are
equal. Given two data sets, each characterized by its mean, standard deviation and
number of data points, we can use some kind of t test to determine whether the means
are distinct, provided that the underlying distributions can be assumed to be normal. All
such tests are usually called Student's t tests, though strictly speaking that name should
only be used if the variances of the two populations are also assumed to be equal; the
form of the test used when this assumption is dropped is sometimes called Welch's t test.
There are different versions of the t test depending on whether the two samples are
o independent of each other (e.g., individuals randomly assigned into two groups),
or
o paired, so that each member of one sample has a unique relationship with a
particular member of the other sample (e.g., the same people measured before
and after an intervention, or IQ test scores of a husband and wife).
o If the calculated p-value is below the threshold chosen for statistical significance
(usually the 0.05 level or 0.005 for high-restriction), then the null hypothesis that
usually states that the two groups do not differ is rejected in favor of an
alternative hypothesis, which typically states that the groups do differ.
• A test of whether the mean of a normally distributed population has a value specified in a
null hypothesis.
• A test of whether the slope of a regression line differs significantly from zero.
Once a t-value is determined, a P value can be found using a table of values from Student's tdistribution. We choose a P value of 0.05 due to the existing measurement bias. The t-value is
calculated given:
t=
x−μ
σ
where
x is the sample mean, μ is the population mean, and σ is the population
n
standard deviation
In our experiment, we choose T-Test instead of Z-Test since we cannot normalize the data and
we have no knowledge of the true standard deviation of population, which are both critical to ZTest. It leaves T-Test as a reasonable alternative in this problem.
By studying both right-tail T-Test and left-tail T-Test, we are able to separate the feature into
either significant increasing or decreasing. Appendix V: Table 26 is a sample of visualized T-Test
result. The red cell indicates a significant increasing feature and a yellow cell indicates a
significantly decreasing feature. We notice that instead of showing significant changing (both-tail
T-Test), combining right-tail and left-tail results gives a more meaningful translation.
5.1.3 High-Throughput Data Analysis in Cell Migration Study
Both feature selection and statistical T-Test are consulted since we looking for high variety and
significant difference between two datasets. To use both techniques as criterions, we introduce a
weighted average. For duplicated experiment with control datasets c of size m and treated
datasets t of size n, the decisive score of feature x is
⎞
⎞⎛ m n
⎛ m n
⎜ ∑∑ ttest (ci ( x ), t j ( x )) ⎟⎜ ∑∑ fsel (ci ( x ), t j ( x )) ⎟
⎟
⎟⎜ i =1 j =1
⎜ i =1 j =1
score( x ) = ⎜
⎟
⎜
⎟
m⋅n
m⋅n
⎟⎟
⎟⎟⎜⎜
⎜⎜
⎠
⎠⎝
⎝
The ‘ttest’ and ‘fsel’ are corresponding to T-Test and feature selection, which returning either one
or zero. Thus, it converts our score system into a weighted percentage on the occurrence of
features.
- 69 -
Function
T-Test
One
Hypothesis is rejected, there is a
significant difference between two
datasets
Feature is an informative feature,
which providing high variety between
two datasets
Feature
Selection
Zero
Hypothesis is not rejected, there is no
significant difference between two
datasets
Feature is an redundant feature, which
providing insufficient variety between two
datasets
The consideration of both criterions appears to be redundant since it is practically believed that
high variety is identical to significant difference. It is true in experiment such as gene profiling.
However, the biological translations of these two criterions in cell migration study are both
practically and theoretically different from gene profiling. Employing both criterions provides an
advantage in our study.
The major reason is to capture two types of behavior changes. A type-one change indicates a
significant value shift (see figure 43) and a type-two change indicates a significant distribution
change (see figure 44).
1.5
1
y1=rand(200,1)
y2=rand(200,1)+0.5
y1=rand(200,1)
y2=0.5*rand(200,1)+0.25
0.9
0.8
0.7
1
0.6
0.5
0.4
0.5
0.3
0.2
0.1
0
0
20
40
60
80
100
120
140
160
180
0
200
Figure 43
0
20
40
60
80
100
120
140
160
180
200
Figure 44
Figure 45 indicates a real type-one change while elongation of cell has significant reduced due to
siRNA4 downregulation (further experiment in section 5.3). Figure 46 indicates a real type-two
change while the size of cells is more widely spreading in treated group, meaning various in the
size is much larger in treated group.
siRNA4 Mean Elongation
0.9
ctrl
treated
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
20
40
60
80
100
120
140
160
Figure 45
- 70 -
180
siRNA4 Mean Area
900
ctrl
treated
800
700
600
500
400
300
200
100
0
20
40
60
80
100
120
140
160
180
Figure 46
The combination of both feature selection and hypothesis test provides an accurate and
reasonable biological translation of behavior measurement. In next two sections, 5.2 and 5.3, two
quantitive experiments of cell migration study using our tracking analysis will be presented.
- 71 -
5.2 MTLn3 Cell Controlled Experiment
KDE tracking algorithm has been confirmed by artificial object test. In this section, its tracking
ability will be tested in real experiment with well-documented cell reaction. Using well-designed
MTLn3 siRNA screening, we intend to confirm the repetition and reliability of phenotypic research
based on tracking analysis
The MTLn3 cells, also known as the mammary adenocarcinoma primary tumor, are used as study
subject of siRNA screening by toxicology section of LACDR. Two cell lines of original MTLn3 cells
are used in our tracking analysis. The first cell line, known as plGFP cell line (see figure 47), is a
cell line that both cytoplasm and nucleolus respond to GFP. The second cell line, known as GFPPaxillin cell line (see figure 48), is a cell line that only cytoplasm responds to GFP.
plGFP cell lines
GFP-Paxillin cell lines
Figure 47
Figure 48
Tracking analysis will be applied on video data collected from both cell lines with several wellstudied mutation. The tracking study of MTLn3 cells is expected to confirm the following issues:
1. EGF 6 stimulation
2. S178A downregulation
3. EGFR 7 or Erb1 Upregulation
Five cell groups are defined:
Group
plGFP
GFP-EGFR
GFP-Paxillin
GFP-Paxillin S178A
GFP-Paxillin S178A EGFR
Description
It has no transfection. It is also used as the control group of GFPEGFR.
GFP-EGFR, which is also known as GFP-ErbB1, is a siRNA
transfection upregulating ErB1.
It has no mutation. It is also used as the control group of GFPPaxillin S178A and GFP-Paxillin S178A EGFR.
The DNA sequence-S178A is downregulated.
The DNA sequence-S178A is downregulated, but with EGFR
upregulated.
Experiment-1: EGF stimulation
6
Epidermal Growth Factor (EGF) is a growth factor that plays an important role in the regulation
of cell growth, proliferation, and differentiation.
7
The epidermal growth factor receptor EGFR is the cell-surface receptor for members of the
epidermal growth factor family (EGF-family) of extracellular protein ligands.
- 72 -
In order to confirm the effect of EGF stimulation, all five groups of cells are treated with EGF
stimulation. A higher cell motility and rapid action are expected with EGF stimulation.
Experiment-2: S178A downregulation
To confirm the correlation between EGFR and S178A, the GFP-Paxillin S178A and the GFPPaxillin S178A EGFR groups have a S178A downregulated. The GFP-Paxillin S178A EGFR has
both S178A downregulation and EGFR upregulation.
Experiment-3: EGFR upregulation
To confirm the effect of EGFR upregulation, GFP-EGFR group and GFP-Paxillin S178A EGFR
group is treated with EGFR upregulation. It is believed that EGFR upregulation may result in
uncontrollable cell division, which is a preposition of cancer.
5.2.1 Experiment Protocol
To reduce cost, three experiments are combined into one plate for scanning. The well content
and locations are given in figure 49 and table 22. When videos are taken from the control group,
EGF stimulation will be added into the same plate and videos will be taken again. Dual positions
are taken from each well, which resulting in a 3 minutes/frame movie in a total length of 90
minutes for control group and 120 minutes for EGF treatment (1.2ul per 2ml). Since the
stimulation effect of EGF is stabilized after 60~120 minutes, there is no significant improvement
for taking video longer than 120 minutes for EGF stimulation.
Figure 49
Table 22
Videos are created in size of 512X512 pixels and tracking analysis is performed on all videos.
5.2.2 Experiment Conclusion
Based on the tracking analysis result of video data (see table 23), we have the following
conclusion:
- 73 -
•
•
•
Compared with control samples, EGF treated samples shows strong expression on the
motility ability and extension of lamellipodium. The stimulation effect is weak effect on
GFP EGFR and GP S178A. The significant decrease in velocity of GP S178A EGFR is
due to several rapid-moving noise, thus measurement cannot be trusted. (see table 23)
Compared with GFP-Paxillin WT and GFP-Paxillin S178A, S178A downregulation (GFPPaxillin S178A) shows strong suppression of motility ability. EGFR upregulation (GFPPaxillin S178A EGFR) has a weak compensation effect on S178A downregulation, which
is an evident of existing correlation between S178A and EGFR. (see table 24)
Compared with plGFP and GFP-EGFR, EGFR upregulation (GFP-EGFR) shows strong
expression on both motility ability and lamellipodium extension frequency. However, the
effect on expression of lamellipodium is insignificant. (see table 25)
Table 23 Control vs. EGF Stimulated Samples
'std
'mean
'area'
extension' velocity'
plGFP
-0.177
-0.128
0.459
GFP EGFR
-0.095
-0.064
0.107
'motionLinearity'
0.046
0.012
GP WT
GP S178A
EGFR
-0.031
-0.094
0.201
0.011
0.235
0.283
-0.106
-0.023
GP S178A
-0.003
0.086
0.188
0.045
freq
area
0.040
0.108
0.026
0.268
0.100
freq
extension
0.632
0.337
0.121
0.010
-0.034
Table 24 Differential ratio Table GFP-PAX WT vs. S178A-EGFR and S178A (Partial)
GP WT ctrl
Std
mean
freq
freq
mean
Diff Ratio
extension'
velocity'
area
extension extension
'area'
GFP Paxillin
-0.420
-0.079
0.020
0.012
-0.324
S178A EGFR -0.354
GFP Paxillin
-0.295
-0.215
-0.296
0.295
0.135
-0.091
S178A
Std
mean
freq
freq
mean
GP WT ctrl
'area'
extension'
velocity'
area
extension extension
T-Test
GFP Paxillin
-1
0
0
0
-1
S178A EGFR -1
GFP Paxillin
-1
-1
-1
1
1
0
S178A
Table 25 plGFP vs. GFP EGFR (Differential ratio and T-Test)
std
std
mean
plGFP ctrl Diff
Ratio
'area' extension' elongation' velocity'
GFP EGFR
-0.169
0.041
0.061
0.841
std
std
mean
plGFP ctrl T-Test
'area' extension' elongation' velocity'
GFP EGFR T-Test
-1
0
0
1
at 0.05
- 74 -
5.3 Large Volume siRNA Screening using MTLn3 plGFP
Large volume siRNA screening or siRNA transfection experiment is a designed experiment to
reveal the function of targeted siRNA. Genetic microarray and cell microarray are two common
approaches used in high-throughput siRNA transfection experiment. The major difference
between genetic microarray and cell microarray is that genetic microarray intends to find
genotypic correlations while cell microarray intends to find phenotypic differences.
In current biology experiment, it is difficult to conduct large volume phenotypic analysis due to the
limitation in analysis software. The phenotypic analysis is performed by either manual observation
or quantitive image analysis [24]. Although considering as automatic approach, quantitive image
analysis shows its limitation in motility estimation and morphology differentiation. In fact,
quantitive image analysis is principally identical to genetic microarray, which only providing a
relative abstract idea of measurements. In general, quantitive image analysis uses averaged
measurement based on all cell in each frame. However, we can avoid such overall estimation by
targeting each cell and looking directly into its motility and morphology estimation during the
entire video. Such action can now be performed by our tracking analysis and large volume data
analysis.
Although both artificial object test and real experiment test show positive feedback on KDE based
tracking analysis, measurement bias is nearly unavoidable in cross-sample experiments. One
major reason is due to the rarity of repetition in biological experiment. Slightly change in
experimental protocol may result in different cell phenotypic behavior.
To assure the usability of tracking analysis in later siRNA screening, we design a small-scale
siRNA screening experiment. The small-scale siRNA screening experiment is designed to
demonstrate the applicability of our tracking analysis in large volume video data. Table 26 shows
the list of siRNA transfection (downregulation) chosen in this experiment. The decision is made
on current literatures while most siRNA have known influence in cell proliferation, differentiation
and transformation. The experiment protocol is given in Appendix VI: Experiment Protocol of
siRNA Screening.
Table 26
label
GFPsiRNA
ctrlsiRNA#2
siRNA1
siRNA2
siRNA3
siRNA4
siRNA5
siRNA6
siRNA7
siRNA8
paxsiRNA
target
GFP
siControl #2
EGFR
LOC364403
PTK2
ILK
FOSL1
VIL2
RDX
LOC305679
PXN-PREDICTED
notes
SRC
Protein tyrosine kinase 2, FAK
VCL
Paxillin
Based on Appendix V, we can directly observe the following changes:
• Downregulation of GFPsiRNA (GFP related) suppresses the speed that cell changes its
shape. Lamellipodium extension is expressed and cell is larger than control group.
Motion linearity is higher than control group.
• Downregulation of siRNA1 (EGFR related) suppresses the extension of lamellipodium,
but increases the size of cells while compared with control groups.
- 75 -
•
•
•
•
•
•
•
•
Downregulation of siRNA2 (SRC related) suppresses the extension of lamellipodium. It
also decreases the speed that cell changes its shape. Cell motion is linear but motion
velocity is lower than control group. Size of cell is larger than control group.
Downregulation of siRNA3 (PTK2 related) suppresses the motion velocity and motion is
non-linear. Shapes of cells are less smooth, which consists of small lamellipodium
extension. Cells are larger than control group.
Downregulation of siRNA4 (ILK related) suppresses the speed that cell changes its
shape.
Downregulation of siRNA5 (FOSL1 related) expresses the extension of lamellipodium, a
higher speed that cell changes its shape can be observed. Cell shapes are more irregular
and larger than control group.
Downregulation of siRNA6 (Vilin2 related) suppresses the extension of lamellipodium.
There is a significant decreasing of lamellipodium activity. Motion linearity of cell is higher
than control group.
Downregulation of siRNA7 (Radixin related) shows no significant effect on phenotypic
behavior of cells.
Downregulation of siRNA8 (Vinculin related) suppresses the extension of lamellipodium.
Cell motion velocity is lower than control.
Downregulation of paxsiRNA (Paxillin related) suppresses the extension of lamellipodium.
There is a significant decreasing of lamellipodium activity. Motion linearity of cell is
increased and cell is larger than control.
Based on known literatures, we can confirm most of our observation from tracking analysis:
• Mutations involving EGFR (siRNA1) could lead to its constant activation, which could
result in uncontrolled cell division. Consequently, mutations of EGFR have been identified
in several types of cancer, and it is the target of an expanding class of anticancer
therapies. [27]
• Mutations in SRC-protein (siRNA2) could be involved in the malignant progression of
colon cancer. [28]
• Activation of PTK2 (siRNA3) may be an important early step in cell growth and
intracellular signal transduction pathways triggered in response to certain neural peptides
or to cell interactions with the extracellular matrix. [26]
• ILK (siRNA4) has been associated with multiple cellular functions including cell migration,
cell proliferation, cell-adhesions, signal transduction. [29]
• FOSL1 (siRNA5) has been implicated as regulators of cell proliferation, differentiation,
and transformation. [30]
• VIL2 (siRNA6) as a member of the ERM protein family, this protein serves as an
intermediate between the plasma membrane and the actin cytoskeleton. It plays a key
role in cell surface structure adhesion, migration, and organization. [33]
• Radixin (siRNA7) is a cytoskeletal protein that may be important in linking actin to the
plasma membrane. It is highly similar in sequence to both ezrin and moesin. [32]
• Vinculin (siRNA8) is a membrane-cytoskeletal protein in focal adhesion plaques that is
involved in linkage of integrin adhesion molecules to the actin cytoskeleton. [31]
By both visual and literature confirmation from our small-scale siRNA screening, tracking analysis
reveals and identify several important hypothesis. We also notice that compared with plGFP, the
Erb1 cell line shows a consistently over-expressing and is not a good study subject for further
screening. The 6-mintues rule is indeed cell line oriented and it is fulfilled by current siRNA
screening. Tracking error is less than 20% on average. The total processing time is 12 hours for
48 movies with length of 120 frames and 100 cells per movie on averages. For 520 movies (the
size of final screening (140 selected siRNAs), it takes less than 6 days while normally it would
take at least one month to manually investigate the movies. The most important improvement is
that we can now obtain a numeric estimation of cellular phenotypic behavior.
- 76 -
6 Conclusion
In this report, we involved four different tracking algorithms. Tracking performances of these
algorithms are compared. KDE tracking is proven less sensitive to shape changes, time
resolution and connectivity error. The error estimation of tracking suggests that KDE tracking is
indeed a robust and sophisticated algorithm. The feedback from Toxicology Section, LACDR
shows that duration of post-experiment analysis is reduced from one month to one week. The
accuracy of data analysis is increased from 30% to 80%. Several important measurements, which
cannot be observed manually, are captured by tracking analysis.
During the study of tracking algorithm, we notice that time-resolution is playing a decisive role in
tracking success and measurement bias. However, no universal rule has been found to
determine the proper time-resolution. The controlled experiment in section 4.7 shows that timeresolution is highly experiment-oriented. Different cell lines may require a re-estimation of a
proper time-resolution. In addition, time-resolution of experiment such as stimulation must be
estimated based on the most excited culture instead of control culture. Our current rule of thumb
for MTLn3 cell lines is to use movie with time interval less than 6 minutes per frame.
Using trajectories obtained from KDE Mean Shift tracking, morphology and motility
measurements can be made. Current measurements are proven accurate in estimating cell
migration. Small controlled experiment such as section 5.2 shows that current measurements
provide an accurate estimation of treatment such as EGF stimulation or S178A downregulation.
In addition, cell behavior can be directly translated from numeric estimation. Significant behavior
change can be easily determined from feature measurements. Large volume data analysis is
proven useful for siRNA screening.
Feature selection and hypothesis test are potential solutions to full automatic cell migration study
in large volume video data. The section 5.3 shows that the combination of feature selection and
hypothesis test provides rapid estimation of significant behavior changes. Current version of our
automatic cell migration analysis, which is under evaluation by toxicology section, LADRC, has
shown great potential in double blind large volume siRNA screening. By the feedbacks from
toxicology section, our automatic approach has tremendously reduced the labor hours and the
conclusion errors. The processing time of cell migration study in siRNA screening is reduced from
one month to one week. Accuracy of conclusion is increased from 30% to 80%. In addition, most
procedures are performed automatically without human interference.
The tracking analysis with MTLn3 cell strongly gives us sufficient confidence to continue further
research on molecular level tracking. The correlation between integrin kinase and cell migration
will be our study subject in future.
- 77 -
7 Discussion
KDE Mean Shift tracking shows a reduced performance on later focal adhesion study. We believe
that KDE Mean Shift is limited to only cellular level tracking. For molecular level tracking, a new
algorithm must be implemented. Study of failure samples in each tracking algorithm shows that
the following factors may be the potential reason of most error in tracking and measurement.
High Similarity of Cells
Most tracking algorithm relies on certain level of dissimilarity between non-alignment object and
similarity between alignment object. However, such principle is not guaranteed in biological
imagining, especially in microscopy imagining. The size, shape and intensity of cells are not
unique and the chance of two identical cells can be extremely high when noise suppression
techniques are used. The reason why KDE mean shift tracking survived is unclear, but it does
provide a certain level of accurate tracking compared with other existing algorithms.
Bias due to Time Resolution
Sampling problem is quite common in time-based, especially time-lapse data. Since
undersampling is practically unavoidable due to microscopy mechanism, we can only reduce the
bias caused by improper time resolution. Our research shows that choice of time resolution is
highly cell-oriented and experiment-oriented. There is no universal and automatic solution for time
resolution choosing. Specialist interferences are required in finding a proper time resolution.
Stochastic Structure Change
Stochastic structure change problem is a combination of between time resolution and cell
activities. Cell migration results in a highly deformation of shape and intensity distribution,
together with low time resolution, it significant increases the dissimilarity between alignment cells
and potentially forces the tracking algorithm to looking for less desire alignment. A high time
resolution is a proper solution for this problem.
Segmentation Limitation
Illumination is the major problem in microscopy imagining. Uneven illumination and inconsistent
intensity of objects are two fields that we have to correct before tracking. Uneven illumination can
be corrected by background subtraction. Yet, it is not easy to correct inconsistent intensity of
objects. Local adaptive threshold shows significant improvement for such problem. Still, certain
level of operator guidance is required, such as deciding the neighborhood size. Automatic
approach can only be applied to videos with following conditions:
1. Microscope flat field illumination
2. Cells have a consistent intensity spreading
3. Cells have a maximum separation between each other
Limitation of Sampling Percentage in Each Well
A common well contains over 50,000 cells while only 100~150 cells can be observed by
microscopy. The sampling percentage of each well is relative low. For siRNA transfection
experiment, only 60%~70% can be transfected while 60~70% of them are correctly transfected.
Meaning for 50,000 cells, there will be at least 15,000 cells that are completely not transfected,
24,500 cells are correctly transfected, and 10,500 cells are incorrectly transfected. Based on
these number, the chance that the observed cells (150 cells) falling in any one of the group is still
quite high. An arbitrary solution is to increase the number of observed cells by reducing the
magnification level. However, magnification level is limited by the basic principle of phenotypic
research since cell structure must be observed. Meaning minimum magnification level is around
20X. Such limitation can be a potential problem for post-experiment analysis.
Regarding all existing problems, tracking analysis does provide important morphology and motility
behavior measurements. Current implementation of KDE mean shift tracking provides valuable
information on cell migration of MTLn3.
- 78 -
There are several problems being partially addressed in our research. These problems will be our
further research directions.
Tracking Analysis for Focal Adhesion Kinase
Showing outstanding performance in cell tracking, KDE Mean Shift has an initial failure in focal
adhesion kinase tracking. The major reason is due to their relative small size and consistent
appearing and disappearing. However, tracking focal adhesion kinase has several advantages
that are not held by cell. Most motion in focal adhesion is nearly linear, which contrary to the
random motion in cell migration. Although change its shape, the motion direction appears to be
on the same direction of motion. Based on these two principles, we will implement particular filter
tracking as a replacement for KDE Mean Shift in study of focal adhesion kinase.
Correlation between Phenotype and Focal Adhesion
Cytomics research will be the future of cell biology. Different from genomics and proteomics,
cytomics is the study of cell systems at a single cell level. It combines all the bioinformatic
knowledge to attempt to understand the molecular architecture and functionality of the cell system.
Much of this is achieved by using molecular and microscopic techniques that allow the various
components of a cell to be visualized as they interact in vivo. By combining the current phenotype
with our potential result of focal adhesion kinase, we will be able to find a systematic way to
reveal the foundation of cytomics.
Multi-Descendent Tracking
One important cell activity is cell differentiation. It strongly violates the basic principle of nearly all
tracking algorithm since most tracking algorithm assuming there is only one true alignment. A
multi-descendent tracking algorithm will provide important information concerning cell
differentiation rate.
Cell Activity Modeling
Although widely agreed, amoebic motion model provides no correlation between cytomorphology
and cellular motility. Yet, it does suggest an implicit binding between lamellipodium extension and
cellular motility. We believe that such correlation does exist but more complicated than amoebic
motion model. Moreover, it may be mathematically modeled by periodic function.
Texture Based Cell Segmentation
Texture based segmentation is widely used in computer vision. It has several advantages in
determine homogeneous regions. Compared with region based, edge base, and region growth
segmentation, texture based segmentation shows several potential in cell segmentation problem.
However, the principle of application texture based segmentation in cell segmentation problem
remains unclear. Obtaining initial module is a decisive factor in texture based cell segmentation.
- 79 -
Appendix I: Workflow of Video Data Analysis
Biological
Experiments and
Time-Lapse Cell
Migration Video
Cell Tracking
Background
subtraction and
flat field correction
KDE Mean Shift
Tracking
Noise
Suppression
Mean Filter
Artificial Object
Test
Gaussian Filter
Quantitive Cell
Video Test
Image
segmentation
Region Based
Segmentation
Raw Feature
Measurements
Edge Based
Segmentation
Generalized
Feature
Measurements
Region Growth
Segmentation
Mathematic binary
morphology
Erosion
Overlap Ratio
Tracking
Pattern
Recognition
Dilation
Feature Selection
Statistical Test
Clustering
Discriminate and
Significance
Analysis
Object Labeling
- 80 -
T-Test
Appendix II: Index List of Feature Measurements
Index
Number
1
2
Feature Name
Feature Description
‘area’
'perimeter'
3
‘std extension'
4
‘std dispersion'
5
6
8
9
10
‘std elongation'
‘mean velocity'
‘mean
motionLinearity'
‘std compact
factor'
‘freq area’
‘freq perimeter’
11
12
‘freq extension’
‘freq dispersion’
13
14
15
16
‘freq elongation’
‘mean extension’
‘mean dispersion’
‘mean elongation’
The average size of cell in pixel
The average perimeter of cell in pixel
It measures how significant the cell changes its extension of
lamellipodium
It measures how significant the cell changes its small structure
(cytoskeleton or lamellipodium)
It measures how significant the cell changes its elongation of
lamellipodium
It measures how fast the cell moves in pixel
It measures how consistent the motion is. A straightforward, Uturn or circular motion will also be considered as linear.
It measures how significant the cell changes its shape while
being compared with a circle
It measures how fast the cell changes its size
It measures how fast the cell changes its perimeter
It measures how fast the cell changes its extension of
lamellipodium
It measures how fast the cell changes its small structures
It measures how fast the cell changes its elongation of
lamellipodium
It measures how extending the cell is on average
It measures how smooth the cell is on average
It measures how elongating the cell is on average
7
For example:
The ‘std elongation’ means the standard deviation of elongation in each frame for one cell. For
whole experiment, the mean of ‘std elongation’ for all cells will be a quite representative
estimation for the entire culture.
Feature 1 and 2 describe the size of cells on average. Upregulation of siRNA such as EGFR may
result in a increasing of these two features while compared with control groups.
Feature 3, 4 and 5 describe how significant the cell changes its shape. Downregulation of S178A
shows a sample of significant decreasing in shape compared with control group.
Feature 9, 10, 11, 12, and 13 describe the activity of cytomorphology changes. Downregulation of
S178A shows that cells has a less extended but quickly changing lamellipodium, which resulting
in a low value in feature 3, 4, and 5 but high value in feature 9, 10, 11, 12, 13.
- 81 -
Appendix III: Artificial Object Test Data Sheet
III-1 Artificial Object Test and Perturbation Experiment
KDE Mean Shift Tracking
minH:20 maxH:40 minW:30 maxW:60 H:30 W:45 Seq:30
Mint
maxT
cell
remote Separated
5
10
5
FALSE
FALSE
5
10
5
TRUE
FALSE
10
20
5
FALSE
FALSE
10
20
5
TRUE
FALSE
5
10
20
FALSE
FALSE
5
10
20
TRUE
FALSE
10
20
20
FALSE
FALSE
10
20
20
TRUE
FALSE
5
10
5
FALSE
TRUE
5
10
5
TRUE
TRUE
10
20
5
FALSE
TRUE
10
20
5
TRUE
TRUE
5
10
20
FALSE
TRUE
5
10
20
TRUE
TRUE
10
20
20
FALSE
TRUE
10
20
20
TRUE
TRUE
- 82 -
Correct%
78.601
99.407
79.5527
96.001
73.78
86.518
64.933
74.403
89.4
100
86.255
99.266
74.084
97.057
68.81
89.474
Overlapping ratio Tracking
minH:20 maxH:40 minW:30 maxW:60 H:30 W:45 Seq:30
minT
maxT
cell
remote Separated
5
10
5
FALSE
FALSE
5
10
5
TRUE
FALSE
10
20
5
FALSE
FALSE
10
20
5
TRUE
FALSE
5
10
20
FALSE
FALSE
5
10
20
TRUE
FALSE
10
20
20
FALSE
FALSE
10
20
20
TRUE
FALSE
5
10
5
FALSE
TRUE
5
10
5
TRUE
TRUE
10
20
5
FALSE
TRUE
10
20
5
TRUE
TRUE
5
10
20
FALSE
TRUE
5
10
20
TRUE
TRUE
10
20
20
FALSE
TRUE
10
20
20
TRUE
TRUE
- 83 -
Correct%
76.687
87.353
77.892
83.642
61.08
81.069
51.766
57.995
79.097
77.914
89.179
80.653
72.226
90.433
63.76
85.39
Appendix IV: MTLn3 Cell Migration Data Sheet
Table 27 Full Data Sheet (Video Data 27-Feb-2008)
plGFP
GFP
EGFR
GP WT
0.21
freq
perimete
r
0.21
freq
extensio
n
0.16
freq
dispersi
on
0.27
freq
elongati
on
0.16
0.22
0.26
0.26
0.32
0.26
'area'
'perimet
er'
'extensio
n'
'dispersi
on'
'elongati
on'
'velocity'
'motionLinearit
y'
'compac
t factor'
freq
area
ctrl
589.17
112.54
0.14
0.01
0.14
0.01
0.89
0.14
egf
484.68
98.76
0.13
0.01
0.12
0.01
0.93
0.10
ctrl
494.01
104.41
0.17
0.01
0.17
0.02
0.91
0.18
0.21
0.24
0.22
0.27
0.22
egf
447.32
96.19
0.16
0.01
0.15
0.02
0.92
0.16
0.23
0.25
0.29
0.29
0.30
ctrl
812.92
132.30
0.13
0.03
0.12
0.01
0.85
0.17
0.23
0.25
0.25
0.26
0.23
egf
787.59
128.86
0.12
0.02
0.10
0.02
0.86
0.14
0.22
0.26
0.29
0.31
0.29
GP
S178A
EGFR
ctrl
547.51
101.64
0.08
0.01
0.08
0.01
0.79
0.10
0.29
0.30
0.28
0.31
0.29
egf
676.02
119.72
0.11
0.02
0.10
0.01
0.77
0.14
0.21
0.22
0.28
0.26
0.27
GP
S178A
ctrl
573.91
107.98
0.10
0.01
0.10
0.01
0.75
0.09
0.27
0.29
0.28
0.31
0.26
egf
572.38
109.38
0.11
0.02
0.10
0.01
0.78
0.14
0.24
0.26
0.27
0.31
0.30
- 84 -
Appendix V: Large Volume siRNA Screening Results
Table 28 plGFP siRNA Screening
movie
plGFP
CTRL
GFPsi
RNA
ctrl
siRNA
2
siRNA
1
siRNA
2
siRNA
3
paxsi
RNA
siRNA
8
siRNA
7
siRNA
6
siRNA
5
siRNA
4
'area'
'perim
eter'
std
extens
ion'
std
disper
sion'
std
elong
ation'
mean
velocit
y'
mean
motio
nLine
arity'
std
comp
act
factor'
freq
area
freq
perim
eter
freq
extens
ion
freq
disper
sion
freq
elong
ation
mean
extens
ion
mean
disper
sion
mean
elong
ation
407.92
86.43
0.231
0.025
0.221
0.036
0.830
0.217
0.463
0.473
0.497
0.497
0.497
0.383
0.018
0.365
463.00
93.42
0.232
0.025
0.223
0.038
0.857
0.227
0.423
0.429
0.457
0.464
0.458
0.422
0.020
0.402
428.24
88.77
0.211
0.019
0.204
0.033
0.840
0.194
0.421
0.435
0.456
0.470
0.456
0.377
0.016
0.361
471.25
92.98
0.206
0.022
0.197
0.027
0.822
0.191
0.482
0.485
0.485
0.499
0.486
0.365
0.017
0.348
450.38
91.01
0.204
0.019
0.196
0.030
0.849
0.185
0.434
0.453
0.451
0.484
0.453
0.375
0.015
0.360
449.49
92.11
0.222
0.026
0.213
0.026
0.816
0.223
0.465
0.476
0.495
0.494
0.498
0.390
0.020
0.370
460.34
92.08
0.206
0.019
0.198
0.030
0.847
0.185
0.422
0.432
0.460
0.477
0.465
0.394
0.015
0.379
425.69
88.58
0.218
0.022
0.208
0.031
0.814
0.195
0.478
0.485
0.501
0.515
0.502
0.373
0.017
0.357
426.63
88.96
0.219
0.026
0.210
0.031
0.805
0.219
0.473
0.489
0.503
0.519
0.508
0.395
0.019
0.376
446.33
91.17
0.206
0.019
0.198
0.034
0.849
0.193
0.409
0.408
0.417
0.425
0.420
0.385
0.016
0.369
489.85
96.53
0.228
0.027
0.218
0.032
0.814
0.230
0.458
0.457
0.473
0.485
0.477
0.414
0.023
0.391
425.56
90.24
0.243
0.029
0.233
0.036
0.833
0.250
0.421
0.426
0.439
0.465
0.439
0.431
0.022
0.409
Area is measured in pixels instead of nm
Perimeter is measured in pixels instead of nm
- 85 -
Table 29 plGFP Differential Ratio
GFPsi
RNA
ctrl
siRNA
2
siRNA
1
siRNA
2
siRNA
3
paxsi
RNA
siRNA
8
siRNA
7
siRNA
6
siRNA
5
siRNA
4
'area'
'perim
eter'
std
extens
ion'
std
disper
sion'
std
elong
ation'
mean
velocit
y'
mean
motio
nLine
arity'
std
comp
act
factor'
freq
area
freq
perim
eter
freq
extens
ion
freq
disper
sion
freq
elong
ation
mean
extens
ion
mean
disper
sion
mean
elong
ation
0.135
0.081
0.005
0.000
0.008
0.055
0.033
0.045
-0.086
-0.093
-0.080
-0.066
-0.080
0.103
0.116
0.103
0.050
0.027
-0.084
-0.231
-0.075
-0.078
0.012
-0.107
-0.090
-0.080
-0.083
-0.053
-0.083
-0.014
-0.115
-0.009
0.155
0.076
-0.107
-0.122
-0.107
-0.260
-0.009
-0.122
0.042
0.026
-0.024
0.004
-0.022
-0.046
-0.060
-0.046
0.104
0.053
-0.117
-0.226
-0.114
-0.161
0.023
-0.151
-0.063
-0.041
-0.092
-0.025
-0.090
-0.021
-0.162
-0.013
0.102
0.066
-0.036
0.029
-0.034
-0.280
-0.016
0.026
0.004
0.006
-0.004
-0.006
0.002
0.020
0.122
0.015
0.129
0.065
-0.107
-0.234
-0.101
-0.158
0.021
-0.150
-0.088
-0.085
-0.074
-0.041
-0.065
0.029
-0.151
0.038
0.044
0.025
-0.057
-0.105
-0.058
-0.137
-0.019
-0.101
0.033
0.026
0.008
0.036
0.008
-0.024
-0.083
-0.021
0.046
0.029
-0.050
0.050
-0.049
-0.149
-0.030
0.010
0.021
0.036
0.013
0.044
0.022
0.032
0.044
0.032
0.094
0.055
-0.105
-0.234
-0.103
-0.043
0.023
-0.113
-0.116
-0.137
-0.160
-0.145
-0.156
0.005
-0.135
0.012
0.201
0.117
-0.013
0.100
-0.013
-0.100
-0.018
0.058
-0.011
-0.033
-0.049
-0.025
-0.041
0.082
0.296
0.071
0.043
0.044
0.052
0.150
0.053
-0.010
0.004
0.150
-0.090
-0.099
-0.116
-0.064
-0.118
0.127
0.213
0.123
Red labeling indicates an increasing ratio
Yellow labeling indicates a decreasing ratio
- 86 -
Table 30 plGFP T-Test at 0.005
GFPsi
RNA
ctrl
siRNA
2
siRNA
1
siRNA
2
siRNA
3
paxsiR
NA
siRNA
8
siRNA
7
siRNA
6
siRNA
5
siRNA
4
'area'
'perim
eter'
std
extensi
on'
std
disper
sion'
std
elonga
tion'
mean
velocit
y'
mean
motion
Lineari
ty'
std
compa
ct
factor'
freq
area
freq
perime
ter
freq
extensi
on
freq
disper
sion
freq
elonga
tion
mean
extensi
on
mean
disper
sion
mean
elonga
tion
1
1
0
0
0
0
1
0
-1
-1
-1
-1
-1
1
0
1
0
0
-1
-1
-1
0
0
-1
-1
-1
-1
0
-1
0
0
0
1
1
-1
0
-1
-1
0
0
0
0
0
0
0
0
0
0
1
1
-1
-1
-1
-1
1
-1
-1
0
-1
0
-1
0
-1
0
1
1
0
0
0
-1
-1
0
0
0
0
0
0
0
1
0
1
1
-1
-1
-1
0
1
-1
-1
-1
-1
-1
-1
0
0
0
0
0
-1
-1
-1
-1
-1
-1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
-1
0
0
0
0
0
0
0
0
0
0
1
-1
-1
-1
0
1
-1
-1
-1
-1
-1
-1
0
0
0
1
1
0
0
0
0
0
0
0
0
-1
0
0
1
1
1
0
0
0
0
0
0
0
0
-1
-1
-1
0
-1
0
0
0
Null hypothesis: There are no significant changes in feature measurement compared with control group
Red labeling indicates a right-tail rejection
Yellow labeling indicates a left-tail rejection
Green labeling indicates that null hypothesis cannot be rejected
- 87 -
Appendix VI: Experiment Protocol of siRNA Screening
30-Mar-2008 plate cells in 48 wells plate
plGFP 12 wells a 10000 cells per well
plGFP 12 wells a 20000 cells per well
Erb1 12 wells a 10000 cells per well
Erb1 12 wells a 20000 cells per well
31-Mar-2008 transfection in 48 wells 100nM/Dharmafect 2
plGFP-20000
Erb1-20000
ctrl
siRNA4
ctrl
siRNA4
GFPsiRNA siRNA5
GFPsiRNA siRNA5
ctrsiRNA2 siRNA6
ctrsiRNA2 siRNA6
siRNA1
siRNA7
siRNA1
siRNA7
siRNA2
siRNA8
siRNA2
siRNA8
siRNA3
paxsiRNA siRNA3
paxsiRNA
1-Apr-2008 cell passage in a collagen coated 24 wells glass bottom plate
Coat the glass bottom plate with 30 ug/ml collagen
Trypse the cells with 50 ul trypsin-EDTA
Add to the glass bottom plate where 500 ul complete MEMα μεδωμ ιν εαχη ωελλ ισ
ctrl
GFPsiRNA
ctrlsiRNA2
siRNA1
siRNA2
siRNA3
siRNA4
siRNA5
siRNA6
siRNA7
siRNA8
paxsiRNA
ctrl
GFPsiRNA
ctrlsiRNA2
siRNA1
siRNA2
siRNA3
siRNA4
siRNA5
siRNA6
siRNA7
siRNA8
paxsiRNA
4-Apr-2008 live cell imaging
8h
Starvation
10h I (-EGF)
Starvation II
11h I (+EGF)
12h II (-EGF)
Starvation III
13h II (+EGF)
14h III (-EGF) Starvation IV
15h III (+EGF)
16h IV (-EGF)
17h IV (+EGF)
18h ready
Two positions are observed for each well.
EGF stimulation 1.2 ul per 2 ml
EGF stock solution 10 ul 1.250ml kaal medium
- 88 -
Reference
[1] Changjiang Yang, Ramani Duraiswami, Daniel DeMenthon, Larry S. Davis: Mean-shift analysis using
quasi Newton methods, Perceptual Interfaces & Reality Laboratory, University of Maryland, College Park,
MD 20742, The United State of American
[2] Fernando de la Torre, Shaogang Gong, and Stephen McKenna: View-Based Adaptive Affine Tracking,
Department of Signal Theory, Universitat Ramon Llull, Spain, Department of Computer Science, Queen
Mary and Westfield College, England, Department of Applied Computing, University of Dundee, Scotland.
[3] Rafael C. Gonzalez, Richard E. Woods, Steven L. Eddins, Digital Image Processing using MATLAB,
Person Prentice Hall Publishing, ISBN 0-13-008519-7
[4] Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, Peter Walter, Molecular
Biology of The Cell, Garland Science Publishing, New York, ISBN 0-8153-4072-9, The United State of
American
[5] Ron Shamir, Tel Aviv University School of Computer Science, Fall 2001-2002, Algorithms in Molecular
Biology, Pairwise alignment, Israel
[6] Kai She1, George Bebis, Haisong Gu, and Ronald Miller, Vehicle Tracking Using On-Line Fusion of
Color and Shape Features, Computer Vision Laboratory, University of Nevada, Reno, NV, Vehicle Design R
& A Department, Ford Motor Company, Dearborn, MI, The United State of American
[7] Gunnar Farneb¨ack, Fast and Accurate Motion Estimation using Orientation Tensors and Parametric
Motion Models, Computer Vision Laboratory, Department of Electrical Engineering, Link¨oping Universit, SE581 83 Link¨oping, Sweden
[8] Qi Zang and Reinhard Klette, Object Classification and Tracking in Video Surveillance, Lecture Notes in
Computer Science, CITR, Computer Science Department, The University of Auckland, Tamaki Campus,
Auckland, New Zealand, Publisher Springer Berlin / Heidelberg, ISBN 978-3-540-40730-0
[9] Mei Han, Wei Xu, Hai Tao, Yihong Gong, An Algorithm for Multiple Object Trajectory Tracking, NEC
Laboratories America, Cupertino, CA, USA, University of California at Santa Cruz, Santa Cruz, CA, USA
[10] Kota Miura, Cell Biology and Biophysics Programme, EMBL, Meyerhofstrasse 1, 69117 Heidelberg,
Germany, Advances in Biochemical Engineering/Biotechnology, Publisher Springer Berlin / Heidelberg, ISBN
978- 3- 540- 23698- 6
[11] Reyes-Aldasoro CC, Akerman S, Tozer GM. Measuring the velocity of fluorescently labeled red blood
cells with a keyhole tracking algorithm, Cancer Research UK Tumor Microcirculation Group, Academic Unit
of Surgical Oncology, The University of Sheffield, K Floor, School of Medicine & Biomedical Sciences, 2008
Feb, Beech Hill Road, Sheffield S10 2RX, UK.
[12] http://en.wikipedia.org/wiki/Main_Page the free encyclopedia
[13] http://www.datastructures.info/the-backtracking-algorithmmethod/ a sample of backtracking
[14] J.M.Fuertes Garcia, M. Lucena Lopez, Region Growth Based on Color Gradients, Dept. de Informatica,
Escuela Politecnica Superior, Universidad de Jaen, Spain, N. Perez de la Blanca, J. Fdez-Valdivia, Dept. de
Ciencias de la Computacion e I.A., ETSI Informatica, Universidad de Granada, Spain.
[15] Raziel Álvarez, Erik Millán and Ricardo Swain-Oropezn Multilevel Seed Region Growth Segmentation,
Carretera al Lago de Guadalupe, Atizapán de Zaragoza, Estado de México, 52926, Mexico, Mechatronics
Research Center (CIMe), Tecnológico de Monterrey, Campus Estado de México, Km 3.5, Lecture Notes in
Computer Science, Copyright 2005, Springer Berlin / Heidelberg, ISBN 978-3-540-29896-0
[16] Dr. Ir. Fons. Verbeeks, Dr. Nies Huijsmans, Lecture Notes for Image Analysis in Microscopy, Leiden
University, LIACS, The Netherlands
[17] Virginia Fonte, Brian Hiester, John Yerg, Jmil Ferguson, Susan Csontos, Michael A. Silverman, and
- 89 -
Gretchen H. Conversion of GFP into a toxic, aggregation-prone protein by C-terminal addition of a short
peptide Christopher D. Link, Stein IBG, University of Colorado, Boulder, CO 80309, USA
[18] Drs. William E. Green, Canny Edge Detection Tutorial, Mechanical Engineering and Mechanics, 3141
Chestnut Street, MEM Department, Room 2-115, Drexel University, Philadelphia, PA 19104, 2002
[19] Cox, I.J.; Boie, R.A.; Wallach, D.A.; "Line recognition", Pattern Recognition Proceedings.10th
International Conference on Volume i, Page(s):639 - 645 vol.1, June 1990
[20] Delport, V.; Liesch, D. Fuzzy-c-mean algorithm for codebook design in vector quantization Electronics
Letters Volume 30, Issue 13, 23 Jun 1994 Page(s):1025 - 1026
[21] The Sage Dictionary of Statistics, Duncan Cramer, Dennis Howitt, 2004, ISBN 076194138X
[22] Yael Paran, Micha Ilan, Yoel Kashman, Sofee Goldstein, Yuvalal Liron, Benjamin Geiger, Zvi Kam, Highthroughput screening of cellular features using high-resolution light-microscopy; Application for profiling drug
effects on cell adhesion, Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot
76100, Israel, HTS Center, Department of Zoology, Tel Aviv University, Tel Aviv 69978, Israel, School of
Chemistry, Tel Aviv University, Tel Aviv 69978, Israel, published by Elsevier Inc, 2007,
doi:10.1016/j.jsb.2006.12.013
[23] Suha Naffar-Abu-Amara, Tal Shay, Meirav Galun, Naomi Cohen, Steven J. Isakoff, Zvi Kam, Benjamin
Geiger, Identification of Novel Pro-Migratory, Cancer-Associated Genes Using Quantitative, MicroscopyBased Screeening, Department of Molecular Cell Biology, Weizmann Institute of Science, Rehovot, Israel,
Department of Physics of Complex Systems, Weizmann Institute of Science, Rehovot, Israel, Department of
Computer Science and Applied Mathematics, Weizmann Institute of Science, Rehovot, Israel,
Massachusetts General Hospital Cancer Center, Boston, Massachusetts, United States of America,
Department of Cell Biology, Harvard Medical School, Boston, Massachusetts, United States of America
[24] Neumann B, Held M, Liebel U, Erfle H, Rogers P, Pepperkok R, Ellenberg J, MitoCheck Project Group,
European Molecular Biology Laboratory, Meyerhofstrasse 1, D-69117 Heidelberg, Germany, Highthroughput RNAi screening by time-lapse imaging of live human cells.
[25] Hinton CV, Avraham S, Avraham HK, Division of Experimental Medicine, Beth Israel Deaconess Medical
Center and Harvard Medical School, Boston, MA 02115, USA, Contributions of Integrin-Linked Kinase to
Breast Cancer Metastasis and Tumorigenesis.
[26] Schaller MD (2001). "Biochemical signals and biological responses elicited by the focal adhesion
kinase.". Biochim. Biophys. Acta 1540 (1): 1-21. PMID 11476890.
[27] A comprehensive pathway map of epidermal growth factor receptor signaling. Molecular Systems
Biology doi:10.1038/msb4100014, 2005 May
[28] Nobel Prize in Physiology or Medicine for 1989 jointly to J. Michael Bishop and Harold E. Varmus for
their discovery of "the cellular origin of retroviral oncogenes".
[29] Persad S, Dedhar S (2004). "The role of integrin-linked kinase (ILK) in cancer progression.". Cancer
Metastasis Rev. 22 (4): 375-84. PMID 12884912.
[30] Herdegen T, Leah JD (1999). "Inducible and constitutive transcription factors in the mammalian nervous
system: control of gene expression by Jun, Fos and Krox, and CREB/ATF proteins.". Brain Res. Brain Res.
Rev. 28 (3): 370-490. PMID 9858769.
[31] Martin TA, Harrison G, Mansel RE, Jiang WG (2004). "The role of the CD44/ezrin complex in cancer
metastasis.". Crit. Rev. Oncol. Hematol. 46 (2): 165-86. PMID 12711360.
[32] Hoeflich KP, Ikura M (2005). "Radixin: cytoskeletal adopter and signaling protein.". Int. J. Biochem. Cell
Biol. 36 (11): 2131-6. doi:10.1016/j.biocel.2003.11.018. PMID 15313460.
[33] Goldmann, W.H., Ingber, D.E. (2001). Intact vinculin protein is Require for control of cell shape, cell
mechanics, and rac-Dependent Lamellipodia Formation. Biochemical and Biophysical Research
Communications. 290: 749-755.
- 90 -
- 91 -