Automated characterization of cell shape changes during amoeboid motility by skeletonization
 Yuan Xiong^{1},
 Cathryn Kabacoff^{2},
 Jonathan FrancaKoh^{2},
 Peter N Devreotes^{2},
 Douglas N Robinson^{2} and
 Pablo A Iglesias^{1}Email author
https://doi.org/10.1186/17520509433
© Xiong et al; licensee BioMed Central Ltd. 2010
Received: 22 October 2009
Accepted: 24 March 2010
Published: 24 March 2010
Abstract
Background
The ability of a cell to change shape is crucial for the proper function of many cellular processes, including cell migration. One type of cell migration, referred to as amoeboid motility, involves alternating cycles of morphological expansion and retraction. Traditionally, this process has been characterized by a number of parameters providing global information about shape changes, which are insufficient to distinguish phenotypes based on local pseudopodial activities that typify amoeboid motility.
Results
We developed a method that automatically detects and characterizes pseudopodial behavior of cells. The method uses skeletonization, a technique from morphological image processing to reduce a shape into a series of connected lines. It involves a series of automatic algorithms including image segmentation, boundary smoothing, skeletonization and branch pruning, and takes into account the cell shape changes between successive frames to detect protrusion and retraction activities. In addition, the activities are clustered into different groups, each representing the protruding and retracting history of an individual pseudopod.
Conclusions
We illustrate the algorithms on movies of chemotaxing Dictyostelium cells and show that our method makes it possible to capture the spatial and temporal dynamics as well as the stochastic features of the pseudopodial behavior. Thus, the method provides a powerful tool for investigating amoeboid motility.
Background
The ability of a cell to change shape is crucial for the proper function of many cellular processes, including cell migration. For example, cells of the immune system move in response to pathogen infections by crawling, which involves cycles of protrusions and contractions that deform the entire cell shape [1]. Traditionally, cell motility has been characterized by a number of different parameters [2]. Some, such as velocity, directional persistence and chemotactic index, are determined by the position of the centroid of the cell. Others, including perimeter, area, roundness and body orientation, describe cellular morphology as the cell migrates. These parameters primarily provide global information about motilityinduced cell shape changes. Though they can be used to distinguish between strains, they may be insufficient to distinguish phenotypes based on pseudopodial protrusions and retractions, which typify amoeboid motility.
We demonstrate that the skeleton can be used to identify pseudopods of moving cells in microscopic images. Moreover, the timevarying evolution of the skeleton can be used to capture the dynamic nature of the pseudopodial protrusions and retractions. We emphasize that the primary goal of this article is not to suggest mechanisms of directed cell motility, but to provide a useful computational tool to analyze amoeboid motility behavior, thus to facilitate the understanding of specific mechanisms in different studies.
Results
Our method consists of a number of separate steps which we describe here individually.
Digital image processing
Different imaging techniques are used to generate movies of motile cells: fluorescent, phasecontrast, and differential interference contrast (DIC). Before further image analysis, two steps are performed to process the original cell movies. The first step is segmentation, in which cell areas are extracted from their background and cell shapes in each frame are acquired. Based on the image features, different algorithms have been designed specifically for the three microscopy techniques mentioned (Appendix A). The second step is boundary smoothing, aimed at removing highfrequency spatial fluctuations of the cell boundary caused by imaging noise. Uniform quadratic Bspline curves, which are smooth curves with continuous first order derivatives and are widely applied in computer aided geometric design, are used to fit cell boundaries [13]. The fitting algorithm is based on the shapespace model and recursive leastsquare estimations [14] (Appendix B).
Skeletonization
The skeleton of a shape is defined based on the concept of maximal ball in the context of 2D Euclidean space; straightforward extensions are possible in 3D space [8]. A ball that is inside a closed boundary curve defined by a planar shape is said to be maximal if it cannot be contained in any other larger ball that is also wholly included in the same curve. The skeleton of a planar shape is defined as the medial axis of the shape, which is the locus of centers of maximal balls in that shape. Note that wherever there is a protrusion from the cell membrane, at least one branch occurs in the skeleton (Figure 1AC).
In practice, because shapes are acquired from digital images, computed skeletons are only approximations of the continuous medial axis. Multiple algorithms that determine the skeleton of a shape defined by a discretized curve have been proposed [15–20]. The medial axis described earlier can be determined based on distance maps, "grassfire" simulations, or Voronoi diagrams. In the distance map, a value is assigned to each point inside the shape corresponding to the distance from this point to the boundary [16, 17]. The skeleton points are defined as the ridges in the constructed contour plot. In grassfire simulations, each point on the boundary of the shape serves as the source of a wave that propagates inward with constant speed [18, 19]. The skeleton is determined as the singular points generated as the waves collide. These two classes of algorithms are considered to provide equivalent results. Finally, Voronoi diagrambased algorithms use the fact that the skeleton of a shape with discrete boundary points can be approximated by a subgraph of the Voronoi diagram [20]. The computational load for this class of algorithms is usually greater. Without further processing, skeletons generated by these algorithms are quite sensitive to shape variations since small perturbations in the boundary curve can result in significant changes to the medial axis. To acquire a more meaningful representation of the shape, boundary smoothing (described above) and branch pruning (discussed below) are usually necessary.
Another class of algorithms that create a skeletons which look similar to, but do not necessarily correspond to the medial axis, uses an approach called thinning. In these algorithms, the boundary pixels of a region are progressively peeled away without changing the topology of the region, until only a skeletal structure remains [12]. Although this method is preferred in some applications, as it is more robust to boundary variations, cell shape information of interest could be lost. For this reason, we use medial axis skeletons and acquire them using the function bwmorph in the Matlab Image Processing Toolbox, the underlying algorithm of which is grassfire simulation (Mathworks, Natick, MA).
Branch pruning
First, short skeleton branches correspond to relatively small variations in shape; conversely, long branches arise from more significant shape features. Thus, the skeleton is pruned by removing small branches. The outer branches of a skeleton are defined as the line segments with at least one terminal that is not connected to any other part of the skeleton. The point linking an outer branch with the inner part of the skeleton is denoted as the root of this branch (Figure 3B, red). With a preset length threshold of p_{threshold} pixels and a length ratio r > 1 (empirical values given in Table 1), each outer branch shorter than p_{threshold} is considered:
1. If the branch is independent; that is, it does not share a root with any other outer branch, then it is removed from the skeleton (Figure 3B, green).
2. Otherwise, if the branch shares a root with another branch, then their respective lengths are compared:

If either the length of the longer branch is less than r times of the length of the shorter one, or the longer one itself is shorter than p_{threshold}, then the two branches are combined into a single one arising from the root (Figure 3B, blue).

If the longer branch is longer than p_{threshold} and longer than r times of the shorter one, then the shorter branch is removed and the longer one retained.
If more than two outer branches share a common root, they are compared and if necessary, combined sequentially. To combine two branches, the point at the average position of the two terminals other than the root is identified and linked to the root. By doing this, not only is shortening of the branches by simply removing the two short branches avoided, but it also guarantees that the combined branch represents the approximate center line of a pseudopod (Figure 3C).
Values of preset parameters in algorithms^{1}.
Variable  Definition  Value used  

Branch pruning  p _{threshold}  Length threshold to identify a short outer branch  1/10 of the length of the cell body^{2} 
r  Ratio threshold to decide whether to combine two branches or to delete the shorter one  1.5  
PrDist _{threshold}  Distance threshold to delete a branch far away from the boundary  1/6 of the length of the cell body  
Backward tracking  T  The maximal length of time apart by which two activities can be considered as from the same pseudopod^{3}  50 s 
α  Weight on the first order difference in computing the cost function in tracking  0.5  
β  Weight on the spatial distance in computing the distance score  0.5  
R  Cell radius  5 μm  
Dist _{threshold}  Threshold on distance score to group two activities into the same pseudopod  1/10 of the perimeter of the cell boundary 
The second criterion for pruning is based on the local curvature of the boundary curve, which can be determined by the reciprocal of the radius of the osculating circle [25]. Typically, each protrusion or retraction occurs near one point where the local curvature is relatively high compared to other parts of the boundary (Figure 3D, red). Because the terminal point of an outer branch is the center of the circle associated with this boundary point, the distance from the terminal point to the boundary gives a measure of the local curvature. When this distance is big, the boundary has a small curvature locally. Thus, the corresponding branch is unlikely to point to a protruding or retracting region of the cell (Figure 3D, blue). To eliminate these branches, a distance threshold (PrDist_{threshold}; Table 1) is set and all branches with terminals (other than their roots) farther away from the cell membrane are removed.
Detecting pseudopodial activity
the expanding and retracting regions of the cell, respectively, from the n^{th} to the (n+1)^{st} frame (Figure 4C). Note that branches of the skeleton from the (n+1)^{st} frame extend out from their root into and those from the n^{th} frame extend out from their root into (Figure 4D).
A protrusion is defined as a region in that includes a branch of S_{n+1}(Figure 4D1). The starting position of a protrusion on the membrane in frame n is set by the intersection of S_{n+1}and the boundary of I_{ n }, and the length of protrusion between frames n and (n +1) is given by the length of branch residing in (L_{1} in Figure 4D1). Similarly, a retraction is defined as a region in that includes a branch of S_{ n }(Figure 4D2). The starting position on the membrane in frame n and the contraction length are defined in the same way.
There may be situations where a skeleton branch does not touch the cell boundary. In these cases, we extend the branch to the membrane and compute the extended length as the step size of the activity.
If the cell is moving in response to an external cue, for example during chemotaxis, we are also interested in how the protrusion or retraction is aligned to the chemoattractant gradient. To compute this, the 0° line is defined as the ray starting from centroid of the cell along which the chemoattractant concentration increases the fastest. The angle of an activity is calculated as the angle from the previously defined 0° line to the line linking the cell centroid and the starting point of this activity on the membrane (Figure 4E).
The time of occurrences, the spatial positions in the starting and ending frames, and the angles from the cell center are recorded for all the detected protrusions and retractions. These data are used to describe the dynamics of these activities.
Backward tracking for pseudopodia lineage classification
It has been argued that different cell strains exhibit dissimilar motility and morphological changes because they extend, retract or split pseudopods in different ways [27, 28]. To analyze the changing patterns of pseudopodia, the lineage of all detected protrusions and retractions needs to be determined. To this end, we develop an automated tracking algorithm that clusters activities into distinct groups, each representing the protruding and retracting history of a single pseudopod.
The basic idea of the algorithm is to model the pseudopod dynamics using a secondorder autoregressive process [29]. More specifically, the following assumptions are made:
1. The changes in spatial position and angle in activities of the same pseudopod from one frame to the next are much smaller than the distance between two different pseudopods;
2. When either the position or angle of the activity changes, the rate of change tends to be constant.
3. If no other protrusion or retraction is observed in an interval of time after the last recorded activity for an existing pseudopod, a pseudopod has terminated. This takes into account the fact that pseudopods have limited lifetimes.
The terms Δ_{1,2} ((x_{1}, y_{1}), (x_{2}, y_{2}), (x_{3}, y_{3})) and Δ_{1,2}(θ_{1}, θ_{2}, θ_{3}) are the spatial and angular distances, respectively, computed from the equation based on the first and secondorder differences (Appendix C). The values β and 1β weigh these two distances, and R is the radius of the cell. The term for the time difference inside the square root sign considers that the probability of two activities from the same pseudopod decreases with the time distance between them. Values used for β and R are given in Table 1. After this distance score is computed from each unassigned activity to all the assigned ones up to T seconds away, the scores are sorted and the minimal one selected. If it is under a given threshold, Dist_{threshold}, the activities are deemed to come from the same pseudopod; otherwise, the unassigned activity is labeled as the latest detected activity of a new pseudopod. Finally, a split is recorded when two activities find the same earlier activity as their common "ancestor" during backward tracking.
Quantities used to describe morphological behaviors.
Classification  Name  Description 

Individual pseudopods  
General  Lifetime  The length of time period when the pseudopod has protrusion or retraction activity 
Number of detected activities  The total number of detected activities during the lifetime  
Net speed  The net displacement of local membrane during the lifetime divided by the time length  
State dynamics  Protrusion/retraction ratio  Ratio of protrusion/retraction activities out of all detected activities 
Protrusion/retraction speed  Total length of protruded/retracted displacement divided by the total length of protrusion/retraction time  
Protrusion/retraction persistence  The longest running time of consistent protrusions/retractions divided by the lifetime  
State persistence  The longest running time of consistent state (either protrusion or retraction) divided by the lifetime  
Angle dynamics  Starting angle  The angle at which the first activity of the pseudopod is detected 
Mean angle  The angle averaged over all activities in the pseudopod  
Linear angle changing rate  The slope of the line that fits best to the activity angle changes  
Detrended angle variance  The variance of activity angles after removal of the linear changing rate  
Splitting behavior  Origin  The pseudopod is generated newly or by splitting from existing ones 
Number of splits  How many splits are generated from the pseudopod during its lifetime  
Cell Level  
Group behavior  Average number of pseudopods  The average number of active pseudopods coexisting on the cell membrane at any given time point 
Ratio of pseudopods from splitting  Number of pseudopods generated from splitting divided by the total number of pseudopods  
Average splitting rate  Total number of splits occurred divided by the time length of the period of movement  
Averaged behavior  Average protrusion/retraction persistence  Computed by averaging the corresponding quantities for individual pseudopods over all pseudopods in the cell. 
Average state persistence  
Average protrusion/retraction speed  
Average net speed 
Correlation analysis
Fluorescence microscopy can be used to obtain a quantitative description of the cellular localization of molecular species that are believed to contribute to pseudopod dynamics. The spatial correlation of these data with the locations of pseudopod activities can provide information as to how different proteins contribute to cell motility.
where , , ⟨g⟩_{k, θ}represents averaging over time and angles, and ⟨g⟩_{ θ }represents averaging over angles only.
Note that the above definitions are invalid if multiple membrane points correspond to the same θ. The situation is usually rare for chemotaxing Dictyostelium cells; in our movies no cell exhibited such a shape in any frame. Moreover, even if this situation were present, it would not influence the pseudopod detection and tracking method previously described, since these membrane points are parts of different branches of the skeleton, and each activity is defined not only by the angle from the center, but also by the spatial location (x, y). Thus, if they represent different pseudopodia activities, they will be identified and analyzed separately.
Assuming that the cell membrane is sampled every Δθ degrees, measuring pseudopodial activity and intensity variables in a movie of N frames, the crosscorrelation computed from Equation 1 is a matrix of the size (360/Δθ) × (2N1), where each row corresponds to an angle value from  (180Δθ)/Δθ° to 180/Δθ°, and each column corresponds to a time shifting length from (N1) to N1 frames. A high absolute value at (Δk, Δθ) means that a linear relationship, either positive or negative, exists between the activities and the intensities located (Δk, Δθ) away. This correlation has also been used to analyze the ordered patterns of spontaneous cell migration [5].
Illustrating the method
To demonstrate the ability of our methods to characterize cell shape changes during amoeboid motility we implement the steps described above to images of chemotaxing Dictyostelium cells. Developed AX3 Dictyostelium cells, an axenic laboratory strain, are placed in a shallow gradient of cAMP, and DIC images are acquired at 1 s/frame [Additional file 1]. In this representative cell, 58 pseudopods are detected. Most (60%, n = 58) are shortlived, lasting no more than 7 s; the remaining longlived pseudopods last at least 9 s and account for over 80% of the total detected activities.
From our analysis, it is clear that the behaviors of individual longlived pseudopods vary dramatically [Additional file 2]. Pseudopods emerge throughout the membrane; furthermore, the angle drifting of pseudopods is also noticeable [Additional file 3]. While some (22%, n = 23) longlived pseudopods remain near the front or rear of the cell, most (70%) drift away from the front for at least half their lifetime. When a pseudopod shifts from the front to the side or back of the cell, it usually transitions from protrusion to retraction. The size of pseudopods also varies considerably; moreover, the lengths of individual activities from the same pseudopods changes over time [Additional file 1].
It has been reported that in Dictyostelium cells chemotaxing in shallow gradients or moving in the absence of external cues, newly generated protrusions usually split from existing pseudopods [4, 27]. Using our method, both the origin of a pseudopod and the number of splits during its lifetime are acquired automatically [Additional file 2]. We find that 74% of all the pseudopods are generated from splitting at an average rate of once every 2.3 seconds. For longlived pseudopods, 78% are generated from splitting, with an average rate of once every 5.6 seconds.
Statistics for phenotype characterization
Pseudopod statistics for different cell strains.
Cell strain  (number of cells/frames)  Avg. # of pseudopods  % from split  Avg. splitting rate (#/min)  Avg. persistence  Avg. state persistence  Ave. speed (μm/min)  

protr.  retr.  protr.  retr.  net  
AX3 (14/1075)  mean  4.2  42  1.54  0.52  0.39  0.82  10.8  9.3  7.8 
STD  0.9  10  0.67  0.06  0.04  0.06  2.2  2.5  2.2  
AX3:tsuA (8/756)  mean  3.6  31  0.99  0.44  0.39  0.74  9.3  7.4  5.7 
STD  0.3  10  0.31  0.05  0.04  0.08  1.8  1.3  1.1  
pvalue  0.07  0.02  0.04  0.01  0.94  0.01  0.12  0.05  0.02  
AX2 (17/1043)  mean  2.9  34  1.67  0.40  0.36  0.67  16.8  11.3  8.5 
STD  0.7  14  1.09  0.09  0.08  0.13  8.4  4.2  3.1  
AX2:dynhp (9/869)  mean  2.8  31  1.66  0.42  0.36  0.72  10.2  8.6  6.6 
STD  0.9  15  1.07  0.07  0.05  0.10  2.2  1.5  1.5  
pvalue  0.67  0.56  0.98  0.54  0.91  0.35  0.03  0.08  0.095 
The chemotactic behaviors of developed Dictyostelium AX2 cells, an axenic laboratory strain, and AX2:dynhp cells, in which dynacortin, an actin crosslinking protein, is depleted by RNAi, are also compared (Table 3). The mutant strain is constructed as previously described [34]. Cells are placed in the gradient of cAMP released from a micropipette containing 1 μM cAMP. Images are captured using a 40× objective at 5 s/frame with spatial resolution of 0.3077 μm/pixel. The AX2:dynhp cells do not show obvious splitting deficiency. The state persistence is also comparable to that of wild type cells. However, this strain produces pseudopods that protrude and retract more slowly: the protrusion speed decreases 39% (p < 0.05, Student's ttest) and the retraction speed decreases 24% (p < 0.1, Student's ttest). The net speed of pseudopods decreases accordingly (p < 0.1, Student's ttest).
Correlation between pseudopods and molecular localizations
A similar analysis of GFPmyosinII reveals that myosinII concentrations along the cell membrane are negatively correlated with protrusion activities (Figure 9EG). However, the relationship between myosinII and retractions is more complex. In some cases (20%, n = 100), myosinII and retractions are highly correlated. However, in most cases (60%), two peaks of high myosinII correlation appear almost symmetrically with respect to the 0° line. In the remaining cells, irregular patterns are observed, not showing obvious maxima or minima. The average GFPmyosinII is significantly higher 90° away from a retraction than at the actual retraction (0.283 vs. 0.034 standard deviations above mean fluorescence intensity at the cell periphery; p < 10^{6}, Student's ttest) (Figure 9H). Thus, the relationship between myosinII concentration and retraction activities may not be simply described as positively or negatively correlated; instead, factors other than myosinII contraction, such as cortical tension, may contribute to pseudopod retraction [7]. To probe this connection further the length changes for all the detected protrusion or retraction activities are plotted as a function of the local fluorescence intensities (Figure 8F). A range of intensities where the length of the protrusions decreases almost linearly with increasing GFPmyosinII intensity is observed. On the other hand, the lengths of retractions increased with myosinII intensity, but with a much smaller slope (Figure 8F). Together, these data suggest that myosinII's role is not just to create retractions, but to suppress lateral protrusions. Interestingly, although the protrusion speed in AX2:dynhp cells is considerably smaller than that of control cells (Table 3), the localized presence of mCherrydynacortin did not contribute greatly to the length of the protrusion, but the retraction length did decrease with increasing of dynacortin concentration on the membrane.
Discussion
Amoeboid motility is a complex cellular process driven by highlyorganized cytoskeletal dynamics involving alternating cycles of pseudopodial protrusions and retractions.
In some cases, the subtle differences in motility between different strains may only be revealed by comparison of pseudopod dynamics. However, because of the highly stochastic fluctuations in the behavior of pseudopods, large data sets may be necessary, making manual detection, identification and quantification impractical. In this study we presented an automated method based on skeletonization that detects and characterizes pseudopodial behavior of cells. These algorithms have been illustrated on movies of chemotactic Dictyostelium cells.
Several methods for detecting pseudopods automatically or semiautomatically in microscopic cell movies have been proposed previously. In the 3D Dynamic Image Analysis System (3DDIAS, Soll Technologies), pseudopods are identified based on manual identification of the "nonparticulate" cytoplasm using DIC images [2, 7]. In addition to being completely automatic, our proposed method has the advantage that only shape information is needed, thus allowing one to apply this to fluorescent images, making correlation studies like those of Figure 8 possible. Moreover, though all of our examples are illustrated in 2D images, the skeletonization technique can be easily expanded to 3D shapes. Other applications for 3D skeletons have been proposed [36–38].
Pseudopodial activities have also been identified by approximating the local curvature of the boundary using a closed polygon formed from a chain of nodes [4]. This form of detection may depend on the magnification, image quality, as well as sampling density. In contrast, the skeletonization technique uses topological and geometrical information from the whole shape, and takes into account not only the local curvature, but also its relationship to the complete path of the membrane, thus allowing great flexibility in applications.
Related research, based on level set methods, tracks the evolution of virtual markers to determine the spatial and temporal dynamics of a region of interest on the cell membrane [3]. In general, this method requires high resolution imaging in both space and time to deliver topologically consistent solutions. In contrast, our method does not need special boundary markers and can be applied to a relatively broader range of resolutions.
For accurate pseudopod detection, the time resolution of the movies should not be too low, so that a considerable part of the cell body overlaps from frame to frame. The valid range should be decided by the size and speed of the specific cell types of interest. In our case, cells move at ~10 μm/min, and the average length of the cell is 10 μm. The minimal frame rate we use is 10 s/frame, which is equivalent to a distance of ~1/6 of their body length between successive frames. At the same time, the spatial resolution must be selected so that desired pseudopodia activities can be identified at least by eye. The spatial resolution for most of our movies is either ~0.3 μm/pixel when using a 40× objective, or ~0.2 μm/pixel when using a 60× objective. The lower the resolution, the more shape features are lost during imaging, and fewer pseudopodial activities, especially subtle ones, can be captured. This will cause a general problem for all detection method, including manual counting. Other parameters in Table 1 are chosen to minimize the discrepancy between results obtained using the automated method and manual observations of pseudopodial activity. These work well with different cell strains, frame rates and magnifications. However, when tracking significantly different cell movements, other parameters may need to be selected.
One possible drawback of relying on topological information is that skeletonization sometimes cannot detect protrusions with bleblike structures, since the local curvature at these structures is relatively low compared to typical activities. This problem may also exist in other automated methods based on local curvature calculations. However, in most of our studies, bleblike structures occur rarely and do not influence the results statistically. For example, there is only one bleb in 101 s in the movie [Additional file 1] at time 59 s, which seems to be a shortlived pseudopod but is missed by the algorithm. Other than the blebs, we do not find shape structures that consistently cause the method to fail in the movies of Dictyostelium cells we analyzed. Attempts have been made to characterize the extension of pseudopodia by amoeboid cells in the absence of external cues as well as in shallow gradient of chemoattractant, based on the assumption that spatial differences in chemoattractant receptor occupancy gives rise to biases in the direction of pseudopod extensions [4, 28].
However, other studies suggest that the bias might come from pseudopod retractions [27]. Here, we view pseudopodial behavior as a dynamic process that includes both protrusions and retractions, and modeled this dynamic behavior using an autoregressive model  in which the state at a given time depends on the historical activities  for each pseudopod. Our results are consistent with models that suggest that, for a Dictyostelium control cell moving in a shallow gradient, a considerable fraction of pseudopods experience both protrusion and retraction, and they tend to retract back when shifting far away from the right direction (Figure 5I). At the same time, the consistent protruding pseudopods may play a more important role in leading the movement of the cell centroid compared to other pseudopods (Figure 5G).
Published results show a discrepancy regarding the production rate of pseudopods, with manual counting reporting less pseudopods than methods that record pseudopods automatically. Our data illustrates that these differences can be attributed to a counting bias. In manual techniques, only the most prominent and persistent protruding activities are identified as pseudopods, whereas automatic methods are able to detect smaller protrusions and retractions. Thus, our method can be applied to quantify both the dominant deformations as well as subtler dynamic perturbations of shape. Our results suggest a higher ratio of "de novo" pseudopods relative to those from splitting during chemotaxis in shallow gradients, compared to previous published results (Additional file 2, [27, 32]). We conjecture that this also comes from the fact that not all meaningful boundary activities are captured in these analyses. We note that for cells chemotaxing to a micropipette, the fraction of pseudopods arising from splitting is even lower (Table 3).
By using the cross correlation method to analyze the molecular drivers of protrusions and retractions, we find that dynacortin, as a marker of Factin, colocalizes with protrusion activities, consistent with the notion that actin polymerization drives protrusions. On the other hand, myosinII is depleted from the front of the cells, and is enriched at the sides of cells. As myosinII contributes ~2030% to cortical viscoelasticity and to cortical tension [39], the localized myosinII modulated increase in cortical tension at the side of cells may thus help to inhibit lateral pseudopods, as previously suggested [40, 41]. This argues that myosinII has a substantial role in enabling of a polarized morphology seen in cells.
Conclusions
In this article we propose an automated method to characterize cell shape changes during amoeboid motility. Based on the skeletonization technique, this method makes full use of both the global and local information of cell shape to detect and track pseudopodial protrusions and retractions, and correlates them with molecular localizations. Using the proposed method, the pseudopodial behavior for single cells can be described, the discrepancies among different strains can be disclosed, and the distinct roles of molecules in driving membrane deformation can be revealed. Thus, it provides a powerful tool to investigate amoeboid motility.
Appendix
A. Image segmentation
We design three different approaches to segment cell areas, according to the imaging techniques used to acquire a specific microscopic movie.
Fluorescence microscopy allows for good contrast between the bright, but of possibly heterogeneous intensity, cell and the dark background. Initial segmentation is achieved by selecting an adaptive threshold intensity that lies between those of the background and cell.
The intensity differences in phasecontrast images come from the anisotropic properties of the medium through which the light travels. Pixel intensities inside the cell are noisy because of the numerous compartments distributed there. Outside the cell, the intensities are considerably smoother. However, a phase ring can be usually observed around the cell periphery, which blurs the boundary and makes the segmentation difficult. To obtain an accurate cell boundary in a phasecontrast image, a deblurring step using the LucyRichardson method [42] is first used to remove the phase ring. Next, gradient mappings of the pixel intensity are generated. High gradients are seen in regions where the cell has considerable topological changes and low gradients are seen at the background. Segmentation is achieved by connecting disparate areas using morphological operations [43–45].
In DIC images, because of the existence of a prisminduced shear direction in the image, bright and dark edges coexist on the cell boundary, but little contrast is observed between these edges [46]. Thus, direct methods, either edge detection algorithms or region segmentations, cannot give satisfactory results for cellshape analysis. A line integration method (line integrated DIC, or LID) has been proposed to transform DIC images to pseudofluorescent images [47]. However, this technique can give rise to stripes as a result of the integration of random noise [48]. Deconvolution and/or the Hilbert transform can be used to correct these images, but neither technique avoids introducing distortions in the shape of the cell. Instead, we use a total variationbased approach, in which an optimal image is found by minimizing an energy function defined by the observed image and the properties of the restored image [49]. Two algorithms: one using texture extraction after LID and the other using denoising before LID, are found to be effective and to give similar results [50]. We note that after LID, the pseudofluorescent image of the cell will be slightly shifted in the prisminduced shear direction when compared to the cell boundary perceived by human eyes (Figure 1C). However, since this is a constant offset everywhere, it does not affect cell shape or the other variables computed during skeletonization.
B. Boundary smoothing
The algorithm contains four steps:
1. Segment the image to get the initial boundary;
2. Sample from the initial boundary to get a group of control points for the construction of a Bspline curve, which defines an initial shape estimate and a template in the shapespace;
3. Search for the feature points along the boundary of the initial shape in the normal directions;
4. Use recursive leastsquare estimation to compute optimal control points according to the positions of the feature points and construct the resultant Bspline curve that fits the boundary best.
As illustrated in Figure 2D, boundary smoothing can eliminate part of spurious branches from the skeleton. We note that splines are appropriate for smoothing as long as the prominent features of a shape are retained. In Figure 2C, the spline faithfully represents the two spots with high local curvatures at the bottom, thus it would not influence the skeletonization results in spite of the missing of small details. However, if a cell has a lot of spiky structures, such as the shape of a fibroblast, and these structures represent the main features of the shape, the smoothing process might be inappropriate as well as unnecessary.
C. First and second order differences
represents the second order difference. These serve as analogues for the first and second timederivatives, respectively, of the function d_{ t }. Thus, the first order difference measures how much the data value changes, and the second order difference measures how much the rate of changes varies over time.
Minimizing the first order difference corresponds to the assumption that changes in spatial position and angle in activities of the same pseudopod from one frame to the next are much smaller than the distance between two different pseudopods. On the other hand, minimizing the second order difference corresponds to the assumption that whenever the position or angle of the activity changes, the rate of change tends to be constant.
Because computing the second order difference for an activity requires three time points, if only two are available, the cost function uses only the first order difference.
D. Strains used
wt: dynhp Ax2:Rep orf+ (HS1000):: pLD1A15SN:dynacortin RNAi
wt: control Ax2:Rep orf+ (HS1000):: pLD1A15SN
wt: tsuA Ax3 ΔtsuA::pJK1:PHcracGFP
wt: control Ax3::pJK1:PHcracGFP
myoII: GFPmyoII, mCherry dynacortin mhcA (HS1):: pBIG:GFPmyosinII; pDRH:mCherrydynacortin
Declarations
Acknowledgements
We are grateful to N. Andrew and R. H. Insall for providing us the movie in Additional file 1. We also thank members of the Devreotes, Robinson and Iglesias Labs for many helpful discussions. This work was supported in part by grants from the NIH, GM71920 (PAI), GM66817 (DNR), GM28007 (PND), GM34933 (PND), and the NSF 0621740 (PAI).
Authors’ Affiliations
References
 Stossel TP: The E. Donnall Thomas Lecture, 1993. The machinery of blood cell movements. Blood. 1994, 84 (2): 367379.PubMedGoogle Scholar
 Soll DR, Voss E, Wessels D, Kuhl S: ComputerAssisted Systems for Dynamic 3D Reconstruction and Motion Analysis of Living Cells. Imaging Cellular and Molecular Biological Functions. Edited by: Shorte SL, Frischknecht F. 2007, 365384. full_text. Heidelberg: SpringerVerlag,View ArticleGoogle Scholar
 Machacek M, Danuser G: Morphodynamic profiling of protrusion phenotypes. Biophysical journal. 2006, 90 (4): 14391452. 10.1529/biophysj.105.070383PubMed CentralView ArticlePubMedGoogle Scholar
 Bosgraaf L, van Haastert PJ, Bretschneider T: Analysis of cell movement by simultaneous quantification of local membrane displacement and fluorescent intensities using Quimp2. Cell motility and the cytoskeleton. 2009, 66 (3): 156165. 10.1002/cm.20338View ArticlePubMedGoogle Scholar
 Maeda YT, Inose J, Matsuo MY, Iwaya S, Sano M: Ordered patterns of cell shape and orientational correlation during spontaneous cell migration. PloS one. 2008, 3 (11): e3734 10.1371/journal.pone.0003734PubMed CentralView ArticlePubMedGoogle Scholar
 Tsukada Y, Aoki K, Nakamura T, Sakumura Y, Matsuda M, Ishii S: Quantification of local morphodynamics and local GTPase activity by edge evolution tracking. PLoS computational biology. 2008, 4 (11): e1000223 10.1371/journal.pcbi.1000223PubMed CentralView ArticlePubMedGoogle Scholar
 Soll DR, Wessels D, Kuhl S, Lusche DF: How a cell crawls and the role of cortical myosin II. Eukaryotic Cell. 2009, 8 (9): 13811396. 10.1128/EC.0012109PubMed CentralView ArticlePubMedGoogle Scholar
 Sonka M, Hlavac V, Boyle R: Image processing, analysis, and machine vision. 2008, Toronto: Thompson Learning, 3,Google Scholar
 Blum H: Biological Shape and Visual Science 1. Journal of theoretical biology. 1973, 38 (2): 205287. 10.1016/00225193(73)901756View ArticlePubMedGoogle Scholar
 Piper J: Interactive ImageEnhancement and Analysis of Prometaphase Chromosomes and Their Band Patterns. Analytical and Quantitative Cytology and Histology. 1982, 4 (3): 233240.Google Scholar
 Wearne SL, Rodriguez A, Ehlenberger DB, Rocher AB, Henderson SC, Hof PR: New techniques for imaging, digitization and analysis of threedimensional neural morphology on multiple scales. Neuroscience. 2005, 136 (3): 661680. 10.1016/j.neuroscience.2005.05.053View ArticlePubMedGoogle Scholar
 Lam L, Lee SW, Suen CY: Thinning Methodologies  a Comprehensive Survey. Ieee Transactions on Pattern Analysis and Machine Intelligence. 1992, 14 (9): 869885. 10.1109/34.161346.View ArticleGoogle Scholar
 Farin GE: Curves and surfaces for CAGD: a practical guide. San Francisco, CA: Morgan Kaufmann, 5,Google Scholar
 Blake A, Isard M: Active contours: the application of techniques from graphics, vision, control theory and statistics to visual tracking of shapes in motion. 1998, London; New York: Springer,View ArticleGoogle Scholar
 Noble PB: Images of cells changing shape: pseudopods, skeletons and motile behavior. Biological Motion. Edited by: Alt W, Hoffmann G. 1990, 4267. New York: SpringerVerlag,View ArticleGoogle Scholar
 Carlo A, Gabriella Sanniti di B: Ridge points in Euclidean distance maps. Pattern Recogn Lett. 1992, 13 (4): 237243. 10.1016/01678655(92)90074A.View ArticleGoogle Scholar
 Ge YR, Fitzpatrick JM: On the generation of skeletons from discrete Euclidean distance maps. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1996, 18 (11): 10551066. 10.1109/34.544075.View ArticleGoogle Scholar
 Leymarie F, Levine MD: Simulating the Grassfire Transform Using an Active Contour Model. Ieee Transactions on Pattern Analysis and Machine Intelligence. 1992, 14 (1): 5675. 10.1109/34.107013.View ArticleGoogle Scholar
 Siddiqi K, Kimia BB, Tannenbaum A, Zucker SW: Shapes, shocks and wiggles. Image and Vision Computing. 1999, 17 (56): 365373. 10.1016/S02628856(98)001309.View ArticleGoogle Scholar
 Ogniewicz RL, Kubler O: Hierarchical Voronoi Skeletons. Pattern Recognition. 1995, 28 (3): 343359. 10.1016/00313203(94)00105U.View ArticleGoogle Scholar
 Bai X, Latecki LJ, Liu WY: Skeleton pruning by contour partitioning with discrete curve evolution. IEEE Transactions on Pattern Analysis and Machine Intelligence. 2007, 29 (3): 449462. 10.1109/TPAMI.2007.59View ArticlePubMedGoogle Scholar
 Choi WP, Lam KM, Siu WC: Extraction of the Euclidean skeleton based on a connectivity criterion. Pattern Recognition. 2003, 36 (3): 721729. 10.1016/S00313203(02)000985.View ArticleGoogle Scholar
 Siddiqi K, Bouix S, Tannenbaum A, Zucker SW: HamiltonJacobi skeletons. International Journal of Computer Vision. 2002, 48 (3): 215231. 10.1023/A:1016376116653.View ArticleGoogle Scholar
 Shaked D, Bruckstein AM: The curve axis. Computer Vision and Image Understanding. 1996, 63 (2): 367379. 10.1006/cviu.1996.0026.View ArticleGoogle Scholar
 do Carmo MP: Differential geometry of curves and surfaces. 1976, Englewood Cliffs, N.J.: PrenticeHall,Google Scholar
 Soll DR, Voss E: Two and threedimensional computer systems for analyzing how animal cells crawl. Motion Analysis of Living Cells. Edited by: Soll DR, Wessels D. 1998, 2552. New York: WileyLiss,Google Scholar
 Andrew N, Insall RH: Chemotaxis in shallow gradients is mediated independently of PtdIns 3kinase by biased choices between random protrusions. Nature Cell Biology. 2007, 9 (2): 193U191. 10.1038/ncb1536View ArticlePubMedGoogle Scholar
 Bosgraaf L, Van Haastert PJ: The ordered extension of pseudopodia by amoeboid cells in the absence of external cues. PloS one. 2009, 4 (4): e5253 10.1371/journal.pone.0005253PubMed CentralView ArticlePubMedGoogle Scholar
 Box GEP, Jenkins GN, Reinsel GC: Time series analysis: forecasting and control. 1994, Englewood Cliffs, N.J.: Prentice Hall, 3,Google Scholar
 Keller PJ, Schmidt AD, Wittbrodt J, Stelzer EH: Reconstruction of zebrafish early embryonic development by scanned light sheet microscopy. Science (New York, NY). 2008, 322 (5904): 10651069.View ArticleGoogle Scholar
 Bosgraaf L, KeizerGunnink I, Van Haastert PJ: PI3kinase signaling contributes to orientation in shallow gradients and enhances speed in steep chemoattractant gradients. Journal of cell science. 2008, 121 (Pt 21): 35893597. 10.1242/jcs.031781View ArticlePubMedGoogle Scholar
 Bosgraaf L, Van Haastert PJ: Navigation of chemotactic cells by parallel signaling to pseudopod persistence and orientation. PloS one. 2009, 4 (8): e6842 10.1371/journal.pone.0006842PubMed CentralView ArticlePubMedGoogle Scholar
 Tang L, FrancaKoh J, Xiong Y, Chen MY, Long Y, Bickford RM, Knecht DA, Iglesias PA, Devreotes PN: tsunami, the Dictyostelium homolog of the Fused kinase, is required for polarization and chemotaxis. Genes & development. 2008, 22 (16): 22782290.View ArticleGoogle Scholar
 Kabacoff C, Xiong Y, Musib R, Reichl EM, Kim J, Iglesias PA, Robinson DN: Dynacortin facilitates polarization of chemotaxing cells. BMC Biol. 2007, 5 (1): 53 10.1186/17417007553PubMed CentralView ArticlePubMedGoogle Scholar
 Iwadate Y, Yumura S: Actinbased propulsive forces and myosinIIbased contractile forces in migrating Dictyostelium cells. Journal of cell science. 2008, 121 (Pt 8): 13141324. 10.1242/jcs.021576View ArticlePubMedGoogle Scholar
 Palágyi K, Balogh E, Kuba A, Halmai C, Erdőhelyi B, Sorantin E, Hausegger K: Sequential 3D Thinning Algorithm and Its Medical Applications. Information Processing in Medical Imaging. Lecture Notes in Computer Science. 2001, 2082: 409415. full_text. full_textView ArticleGoogle Scholar
 Scellier D, Boire JY, Thouly C, Maublant J: Application of skeletonization algorithms for myocardial spect quantification. Discrete Geometry for Computer Imagery. 1996, 1176: 227236.View ArticleGoogle Scholar
 Xie R, Thompson R, Perucchio R: A topologypreserving parallel 3D thinning algorithm for extracting the curve skeleton. Pattern Recognition. 2003, 36: 15291544. 10.1016/S00313203(02)003485.View ArticleGoogle Scholar
 Reichl EM, Ren Y, Morphew MK, Delannoy M, Effler JC, Girard KD, Divi S, Iglesias PA, Kuo SC, Robinson DN: Interactions between myosin and actin crosslinkers control cytokinesis contractility dynamics and mechanics. Curr Biol. 2008, 18 (7): 471480. 10.1016/j.cub.2008.02.056PubMed CentralView ArticlePubMedGoogle Scholar
 Knecht DA, Shelden E: Threedimensional localization of wildtype and myosin II mutant cells during morphogenesis of Dictyostelium. Dev Biol. 1995, 170 (2): 434444. 10.1006/dbio.1995.1227View ArticlePubMedGoogle Scholar
 Shelden E, Knecht DA: Dictyostelium cell shape generation requires myosin II. Cell motility and the cytoskeleton. 1996, 35 (1): 5967. 10.1002/(SICI)10970169(1996)35:1<59::AIDCM5>3.0.CO;2DView ArticlePubMedGoogle Scholar
 Carasso A: Linear and nonlinear image deblurring: A documented study. SIAM journal on numerical analysis. 1999, 36 (6): 16591689. 10.1137/S0036142997320413.View ArticleGoogle Scholar
 Jähne B: Digital image processing, 6th rev. and ext. edn. 2005, Berlin; New York: Springer,Google Scholar
 Serra JP: Image analysis and mathematical morphology. 1982, London; New York: Academic Press,Google Scholar
 Soille P: Morphological image analysis: principles and applications. 2003, Berlin; New York: Springer, 2,Google Scholar
 Murphy DB: Fundamentals of light microscopy and electronic imaging. 2001, New York: WileyLiss,Google Scholar
 Kam Z: Microscopic differential interference contrast image processing by line integration (LID) and deconvolution. Bioimaging. 1998, 6 (4): 166176. 10.1002/13616374(199812)6:4<166::AIDBIO2>3.0.CO;2Y.View ArticleGoogle Scholar
 Heise B, Sonnleitner A, Klement EP: DIC image reconstruction on large cell scans. Microscopy research and technique. 2005, 66 (6): 312320. 10.1002/jemt.20172View ArticlePubMedGoogle Scholar
 Chan TF, Shen J: Image processing and analysis: variational, PDE, wavelet, and stochastic methods. Philadelphia: Society for Industrial and Applied Mathematics. 2005,Google Scholar
 Vese LA, Osher SJ: Modeling Textures with Total Variation Minimization and Oscillating Patterns in Image Processing. J Sci Comput. 2003, 19 (13): 553572. 10.1023/A:1025384832106.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.