- Methodology article
Automated characterization of cell shape changes during amoeboid motility by skeletonization
BMC Systems Biologyvolume 4, Article number: 33 (2010)
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.
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.
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.
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 . Traditionally, cell motility has been characterized by a number of different parameters . 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 motility-induced 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.
Recently, there has been much interest in developing means for processing microscopic images of motile cells and acquiring local morphological information automatically or semi-automatically [2–7]. Here we describe a series of automated methods used to characterize both local morphological changes and statistical features during amoeboid locomotion based on the skeleton of a planar shape . Skeletonization, also known as the medial axis transform, is a technique in image processing used to reduce a binary shape into a series of connected lines - the skeleton - that roughly maintains the form of the shape (Figure 1A). This thin-line representation of shape has attracted considerable attention [9, 10]. For example, skeletons have been used to measure the lengths and numbers of junctions of tubule complexes in in-vitro angiogenesis assays and to analyze neuronal structures . Part of the technique's attraction is that skeletons of elongated shape patterns, which are frequently observed in many organisms and biological structures, appear to be close to those perceived by humans . Furthermore, the skeleton facilitates shape analysis and uses less data than the original shape. Though skeletonization has long been used to analyze images in cell biology [9, 10], it has not been applied to track dynamic information about cellular shape.
We demonstrate that the skeleton can be used to identify pseudopods of moving cells in microscopic images. Moreover, the time-varying 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.
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, phase-contrast, 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 high-frequency spatial fluctuations of the cell boundary caused by imaging noise. Uniform quadratic B-spline 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 . The fitting algorithm is based on the shape-space model and recursive least-square estimations  (Appendix B).
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 . 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 1A-C).
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 diagram-based algorithms use the fact that the skeleton of a shape with discrete boundary points can be approximated by a sub-graph of the Voronoi diagram . 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 . 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).
Boundary smoothing helps to eliminate spurious branches that arise from noisy images (Figure 2). Nevertheless, further branch pruning is necessary to map each pseudopod perceived by human eyes to exactly one branch in the skeleton. As the name implies, branch pruning helps to minimize the number of distinct elements in the skeleton description of the cell. Many different pruning methods have been developed [20–24]. All these techniques work by defining a metric of significance to skeleton points. In the case of motility-induced cell shape changes, two criteria are designed to prune spurious branches (Figure 3).
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 pthreshold pixels and a length ratio r > 1 (empirical values given in Table 1), each outer branch shorter than pthreshold 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 pthreshold, then the two branches are combined into a single one arising from the root (Figure 3B, blue).
If the longer branch is longer than pthreshold 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).
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 . 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 (PrDistthreshold; Table 1) is set and all branches with terminals (other than their roots) farther away from the cell membrane are removed.
Detecting pseudopodial activity
As cells move, their shapes deform leading to temporal changes in the corresponding skeletons (Figure 4). We seek to describe this dynamic process and to establish an algorithm to determine whether changes in cell shape correspond to protrusions or retractions. Shape differences refer to the regions where the shape changes during motility from one frame to the next . They can be "positive" or "negative" flow, depending on whether the areas are growing or shrinking. However, owing to image noise and random perturbations of membrane shape, differences are not always well defined. Thus, shape differences are combined with skeletons to detect protrusion and retraction activities.
The process is initiated by segmenting successive images of a moving cell (Figure 4A). Assume that I n and In+1are the binary segmentations at the nth and (n+1)st frames respectively, where the pixels within the cell are set to 1 and those of the background to 0. Denote the pruned skeletons in frames n and (n + 1) by S n and Sn+1respectively (Figure 4B). Let D be the 2D space where all the image frames of a movie are defined, and p be a given pixel in D. Define by
the expanding and retracting regions of the cell, respectively, from the nth 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 nth frame extend out from their root into (Figure 4D).
A protrusion is defined as a region in that includes a branch of Sn+1(Figure 4D1). The starting position of a protrusion on the membrane in frame n is set by the intersection of Sn+1and 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 (L1 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 second-order autoregressive process . 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.
Because pseudopods are believed to split frequently during movement, they are tracked backwards in time . It is necessary to decide whether each protrusion or retraction detected at a given time point is part of an ongoing pseudopod or of a new pseudopod. To achieve this, a weighted distance metric between activities that takes into account differences in the spatial location (x, y), the angle (θ), and the time (t) is defined. More specifically, suppose that two activities p1 = (x1, y1, θ1, t1) and p2 = (x2, y2, θ2, t2), with t2>t1 have been identified and it has been deemed that p2 is associated with activity p3 = (x3, y3, θ3, t3), with t3>t2. The distance score between p1 and p2 is computed as
The terms Δ1,2 ((x1, y1), (x2, y2), (x3, y3)) and Δ1,2(θ1, θ2, θ3) are the spatial and angular distances, respectively, computed from the equation based on the first- and second-order 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, Distthreshold, 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.
Using the methods described above, we calculate a number of quantities for individual pseudopods and single cells (Table 2). Computational results for these parameters and the comparisons between different cell strains are presented in detail below.
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.
For a given cell during a given period of time, an activity variable, A(k, θ), is defined as the length of the activity if an activity is detected at frame number k, representing time, and angle θ of the protrusion or retraction activity on the membrane. If no activity is detected, then A(k, θ) = 0. An intensity variable, I(k, θ), is defined corresponding to the relative intensity value at angle θ and frame k on the membrane. The cross-correlation function between these two signals is given by
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 cross-correlation computed from Equation 1 is a matrix of the size (360/Δθ) × (2N-1), where each row corresponds to an angle value from - (180-Δθ)/Δθ° to 180/Δθ°, and each column corresponds to a time shifting length from -(N-1) to N-1 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 .
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 short-lived, lasting no more than 7 s; the remaining long-lived 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 long-lived 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) long-lived 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].
Despite huge differences in positions and sizes among the long-lived pseudopods, their protrusion/retraction state dynamics fall into three main activity patterns: consistent protrusions, consistent retractions, or alternating periods of protrusions and retractions (Figure 5 and Additional file 4). The average state persistence for this cell is 0.80, which means that pseudopods tends either to protrude or to retract during most of their lifetime. We observed that whether a pseudopod extends or retracts depends on the angle of the pseudopod relative to the chemoattractant gradient. Pseudopods near the front have higher probability of maintaining fast growth (Figures 5D, G). This probability decreases if the pseudopod drifts away from the front (Figures 5F, I). Similarly, the probability of retraction is low at the front of the cell, but increases as the pseudopod approached the back (Figures 5E, H). These observations are consistent with previous results [27, 31].
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 long-lived pseudopods, 78% are generated from splitting, with an average rate of once every 5.6 seconds.
Previously published results in quantifying pseudopod production rates have shown significant differences. Using manual identification, on average, a new pseudopod is produced every 23 s for AX3 cells moving in shallow gradient . Using automated methods based on local curvatures of sampled boundaries, this number drops to around 14 s for AX3 cells in shallow gradient and 13 s in the absence of chemoattractant [28, 32]. Although the experimental conditions might play a role, it is possible that this discrepancy comes from the different ways of identifying pseudopods. To investigate this, the pseudopodia dynamics of the same AX3 cell at different levels are examined (Figure 6). First, the smallest number of pseudopods that cover no less than 50% of all protrusion and retraction activities is selected. Using this criterion, only seven long-lasting pseudopods are identified, of which five contained significant protrusion processes. These pseudopods give a production rate of ~20 s per protruding pseudopod, a number similar to that obtained using manual identification . When the activity percentage level is set to 60%, the number of pseudopods detected with significant protrusions increases to eight. These have an average interval of ~13 s.
Statistics for phenotype characterization
To show how our methods can be used to distinguish different cell strains, movies from two Dictyostelium strains: AX3 and AX3:tsuA mutant are analyzed. AX3:tsuA strain, in which the Dictyostelium homolog of the Fused kinase TsuA is knocked out from the parental line AX3, is constructed as previously described . Undeveloped cells from both strains are placed in the gradient of folic acid released from a micropipette containing 1 M folic acid. Images are captured using a 60× objective at 10 s/frame with a spatial resolution of 0.1852 μm/pixel. The pseudopodial behaviors of individual cells are computed and compared between the two strains (Table 3). AX3:tsuA cells show several chemotaxis defects when responding to folic acid. One of the most significant (p < 0.05, Student's t-test) defect is the ability to split pseudopods: the ratio of newly generated pseudopods arising from splits decreases 27% and the splitting rate decreases 36%. The other discrepancy is the state persistence. The protrusion persistence decreases 15% and the overall state persistence decreases 10% in AX3:tsuA cells. These data suggest that these cells have difficulty generating consistent protrusions. The net speed of pseudopods, which decreases 27%, can also be related to the lower protrusion persistence: pseudopods tend to wander around instead of pulling out the membrane far away.
The chemotactic behaviors of developed Dictyostelium AX2 cells, an axenic laboratory strain, and AX2:dynhp cells, in which dynacortin, an actin cross-linking protein, is depleted by RNAi, are also compared (Table 3). The mutant strain is constructed as previously described . 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 t-test) and the retraction speed decreases 24% (p < 0.1, Student's t-test). The net speed of pseudopods decreases accordingly (p < 0.1, Student's t-test).
To contrast the efficiency of movement during chemotaxis, it is useful to look at the angle distribution of pseudopods. A histogram of protrusion angles for each of the four strains described above is plotted, and a Gaussian curve to characterize the distribution is fitted (Figure 7). The mean angle of activities indicates whether cells are able to sense the external cues and generate different responses at the front and the back, while the variance measures how efficiently the activities are produced to enable the directed movement of the cell body. The data reveal that vegetative AX3:tsuA cells do not align their protrusions with the chemoattractant gradient. On the other hand, developed AX2:dynhp cells orient their protrusions up the gradient, but not as efficiently as their respective control cells.
Correlation between pseudopods and molecular localizations
To investigate the molecular drivers of pseudopod formation, the correlation between the time and location of pseudopod activities with fluorescently-labeled proteins is calculated. An AX3-based cell strain is created where myosin-II is tagged using GFP, and, simultaneously, dynacortin is tagged using mCherry. At the same time, both endogenous expressions are confirmed to be depleted completely. It is commonly accepted that myosin-II enriched at the posterior during chemotaxis, generating contractile force at the back and squeezing the cell body forward [7, 35]. On the other hand, dynacortin is normally concentrated at the leading edge where it is thought to cooperate with actin to influence cortical viscoelasticity . In our study, developed cells are placed in a chemoattractant gradient, and images of chemotaxing cells are obtained using a dual-emission microscope (Figure 8A). The respective fluorescent concentrations are detected around the cellular membrane (Figures 8B, D), and pseudopods are identified and classified as to whether they are expanding or contracting (Figure 8C). The mCherry-dynacortin is found consistently at the front of the cell, aligned with the chemoattractant gradient. In contrast, GFP-myosin-II is found to be mostly at the side.
To determine how these protein localizations aligned with the pseudopodial activities, the cross-correlation functions (Equation 1) between the two fluorescent intensities and the protrusion or retraction activities are computed (Figure 8E). In nearly all cells displaying good mCherry expressions, peaks in dynacortin concentration occurred at the same time and location as protrusions. On the other hand, retractions and regions of high dynacortin concentration are highly negatively correlated. To probe further these correlations we compute the correlation coefficient between dynacortin fluorescence and activity length to be 0.623 with 99% confidence interval of 0.595 to 0.649 (Figure 9A). On the other hand, when compared to the fluorescence at the other side of the cell (180° away) the correlation is negative: -0.577 with 99% confidence interval of -0.606 to -0.577 (Figure 9B). When measured cell to cell, all cells (n = 37) have positive correlations between the length of the activity and fluorescence at the location of the activity, and negative correlation with the activity (180° away) (Figure 9C). The average mCherry intensity at the site of protrusions is significantly higher than at retractions (1.06 vs. -0.72 standard deviations above mean fluorescence intensity at the cell periphery; p < 10-6, Student's t-test).
A similar analysis of GFP-myosin-II reveals that myosin-II concentrations along the cell membrane are negatively correlated with protrusion activities (Figure 9E-G). However, the relationship between myosin-II and retractions is more complex. In some cases (20%, n = 100), myosin-II and retractions are highly correlated. However, in most cases (60%), two peaks of high myosin-II 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 GFP-myosin-II 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 t-test) (Figure 9H). Thus, the relationship between myosin-II concentration and retraction activities may not be simply described as positively or negatively correlated; instead, factors other than myosin-II contraction, such as cortical tension, may contribute to pseudopod retraction . 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 GFP-myosin-II intensity is observed. On the other hand, the lengths of retractions increased with myosin-II intensity, but with a much smaller slope (Figure 8F). Together, these data suggest that myosin-II'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 mCherry-dynacortin did not contribute greatly to the length of the protrusion, but the retraction length did decrease with increasing of dynacortin concentration on the membrane.
Amoeboid motility is a complex cellular process driven by highly-organized 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 semi-automatically in microscopic cell movies have been proposed previously. In the 3D Dynamic Image Analysis System (3D-DIAS, 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 . 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 . 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 bleb-like 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, bleb-like 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 short-lived 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 . 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 F-actin, colocalizes with protrusion activities, consistent with the notion that actin polymerization drives protrusions. On the other hand, myosin-II is depleted from the front of the cells, and is enriched at the sides of cells. As myosin-II contributes ~20-30% to cortical viscoelasticity and to cortical tension , the localized myosin-II 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 myosin-II has a substantial role in enabling of a polarized morphology seen in cells.
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.
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 phase-contrast 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 phase-contrast image, a deblurring step using the Lucy-Richardson method  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 prism-induced shear direction in the image, bright and dark edges coexist on the cell boundary, but little contrast is observed between these edges . Thus, direct methods, either edge detection algorithms or region segmentations, cannot give satisfactory results for cell-shape analysis. A line integration method (line integrated DIC, or LID) has been proposed to transform DIC images to pseudo-fluorescent images . However, this technique can give rise to stripes as a result of the integration of random noise . 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 variation-based 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 . 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 . We note that after LID, the pseudo-fluorescent image of the cell will be slightly shifted in the prism-induced 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 B-spline curve, which defines an initial shape estimate and a template in the shape-space;
3. Search for the feature points along the boundary of the initial shape in the normal directions;
4. Use recursive least-square estimation to compute optimal control points according to the positions of the feature points and construct the resultant B-spline 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
Let d t be the measured data point (spatial position or angle value) at time t. Similarly, and are the measurements at time t- Δt2 and t+ Δt1, respectively. We define two measurements describing how data changes over time. The first
defines the first order difference over two frames. The second
represents the second order difference. These serve as analogues for the first and second time-derivatives, 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.
We define a cost function, used in the tracking algorithm that penalizes large changes in both these measures, weighed by constants α and 1-α:
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:PHcrac-GFP
wt: control Ax3::pJK1:PHcrac-GFP
myoII: GFPmyoII, mCherry dynacortin mhcA (HS1):: pBIG:GFP-myosin-II; pDRH:mCherry-dynacortin
Stossel TP: The E. Donnall Thomas Lecture, 1993. The machinery of blood cell movements. Blood. 1994, 84 (2): 367-379.
Soll DR, Voss E, Wessels D, Kuhl S: Computer-Assisted Systems for Dynamic 3D Reconstruction and Motion Analysis of Living Cells. Imaging Cellular and Molecular Biological Functions. Edited by: Shorte SL, Frischknecht F. 2007, 365-384. full_text. Heidelberg: Springer-Verlag,
Machacek M, Danuser G: Morphodynamic profiling of protrusion phenotypes. Biophysical journal. 2006, 90 (4): 1439-1452. 10.1529/biophysj.105.070383
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): 156-165. 10.1002/cm.20338
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.0003734
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.1000223
Soll DR, Wessels D, Kuhl S, Lusche DF: How a cell crawls and the role of cortical myosin II. Eukaryotic Cell. 2009, 8 (9): 1381-1396. 10.1128/EC.00121-09
Sonka M, Hlavac V, Boyle R: Image processing, analysis, and machine vision. 2008, Toronto: Thompson Learning, 3,
Blum H: Biological Shape and Visual Science 1. Journal of theoretical biology. 1973, 38 (2): 205-287. 10.1016/0022-5193(73)90175-6
Piper J: Interactive Image-Enhancement and Analysis of Prometaphase Chromosomes and Their Band Patterns. Analytical and Quantitative Cytology and Histology. 1982, 4 (3): 233-240.
Wearne SL, Rodriguez A, Ehlenberger DB, Rocher AB, Henderson SC, Hof PR: New techniques for imaging, digitization and analysis of three-dimensional neural morphology on multiple scales. Neuroscience. 2005, 136 (3): 661-680. 10.1016/j.neuroscience.2005.05.053
Lam L, Lee SW, Suen CY: Thinning Methodologies - a Comprehensive Survey. Ieee Transactions on Pattern Analysis and Machine Intelligence. 1992, 14 (9): 869-885. 10.1109/34.161346.
Farin GE: Curves and surfaces for CAGD: a practical guide. San Francisco, CA: Morgan Kaufmann, 5,
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,
Noble PB: Images of cells changing shape: pseudopods, skeletons and motile behavior. Biological Motion. Edited by: Alt W, Hoffmann G. 1990, 42-67. New York: Springer-Verlag,
Carlo A, Gabriella Sanniti di B: Ridge points in Euclidean distance maps. Pattern Recogn Lett. 1992, 13 (4): 237-243. 10.1016/0167-8655(92)90074-A.
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): 1055-1066. 10.1109/34.544075.
Leymarie F, Levine MD: Simulating the Grassfire Transform Using an Active Contour Model. Ieee Transactions on Pattern Analysis and Machine Intelligence. 1992, 14 (1): 56-75. 10.1109/34.107013.
Siddiqi K, Kimia BB, Tannenbaum A, Zucker SW: Shapes, shocks and wiggles. Image and Vision Computing. 1999, 17 (5-6): 365-373. 10.1016/S0262-8856(98)00130-9.
Ogniewicz RL, Kubler O: Hierarchical Voronoi Skeletons. Pattern Recognition. 1995, 28 (3): 343-359. 10.1016/0031-3203(94)00105-U.
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): 449-462. 10.1109/TPAMI.2007.59
Choi WP, Lam KM, Siu WC: Extraction of the Euclidean skeleton based on a connectivity criterion. Pattern Recognition. 2003, 36 (3): 721-729. 10.1016/S0031-3203(02)00098-5.
Siddiqi K, Bouix S, Tannenbaum A, Zucker SW: Hamilton-Jacobi skeletons. International Journal of Computer Vision. 2002, 48 (3): 215-231. 10.1023/A:1016376116653.
Shaked D, Bruckstein AM: The curve axis. Computer Vision and Image Understanding. 1996, 63 (2): 367-379. 10.1006/cviu.1996.0026.
do Carmo MP: Differential geometry of curves and surfaces. 1976, Englewood Cliffs, N.J.: Prentice-Hall,
Soll DR, Voss E: Two- and three-dimensional computer systems for analyzing how animal cells crawl. Motion Analysis of Living Cells. Edited by: Soll DR, Wessels D. 1998, 25-52. New York: Wiley-Liss,
Andrew N, Insall RH: Chemotaxis in shallow gradients is mediated independently of PtdIns 3-kinase by biased choices between random protrusions. Nature Cell Biology. 2007, 9 (2): 193-U191. 10.1038/ncb1536
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.0005253
Box GEP, Jenkins GN, Reinsel GC: Time series analysis: forecasting and control. 1994, Englewood Cliffs, N.J.: Prentice Hall, 3,
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): 1065-1069.
Bosgraaf L, Keizer-Gunnink I, Van Haastert PJ: PI3-kinase signaling contributes to orientation in shallow gradients and enhances speed in steep chemoattractant gradients. Journal of cell science. 2008, 121 (Pt 21): 3589-3597. 10.1242/jcs.031781
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.0006842
Tang L, Franca-Koh 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): 2278-2290.
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/1741-7007-5-53
Iwadate Y, Yumura S: Actin-based propulsive forces and myosin-II-based contractile forces in migrating Dictyostelium cells. Journal of cell science. 2008, 121 (Pt 8): 1314-1324. 10.1242/jcs.021576
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: 409-415. full_text. full_text
Scellier D, Boire J-Y, Thouly C, Maublant J: Application of skeletonization algorithms for myocardial spect quantification. Discrete Geometry for Computer Imagery. 1996, 1176: 227-236.
Xie R, Thompson R, Perucchio R: A topology-preserving parallel 3D thinning algorithm for extracting the curve skeleton. Pattern Recognition. 2003, 36: 1529-1544. 10.1016/S0031-3203(02)00348-5.
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): 471-480. 10.1016/j.cub.2008.02.056
Knecht DA, Shelden E: Three-dimensional localization of wild-type and myosin II mutant cells during morphogenesis of Dictyostelium. Dev Biol. 1995, 170 (2): 434-444. 10.1006/dbio.1995.1227
Shelden E, Knecht DA: Dictyostelium cell shape generation requires myosin II. Cell motility and the cytoskeleton. 1996, 35 (1): 59-67. 10.1002/(SICI)1097-0169(1996)35:1<59::AID-CM5>3.0.CO;2-D
Carasso A: Linear and nonlinear image deblurring: A documented study. SIAM journal on numerical analysis. 1999, 36 (6): 1659-1689. 10.1137/S0036142997320413.
Jähne B: Digital image processing, 6th rev. and ext. edn. 2005, Berlin; New York: Springer,
Serra JP: Image analysis and mathematical morphology. 1982, London; New York: Academic Press,
Soille P: Morphological image analysis: principles and applications. 2003, Berlin; New York: Springer, 2,
Murphy DB: Fundamentals of light microscopy and electronic imaging. 2001, New York: Wiley-Liss,
Kam Z: Microscopic differential interference contrast image processing by line integration (LID) and deconvolution. Bioimaging. 1998, 6 (4): 166-176. 10.1002/1361-6374(199812)6:4<166::AID-BIO2>3.0.CO;2-Y.
Heise B, Sonnleitner A, Klement EP: DIC image reconstruction on large cell scans. Microscopy research and technique. 2005, 66 (6): 312-320. 10.1002/jemt.20172
Chan TF, Shen J: Image processing and analysis: variational, PDE, wavelet, and stochastic methods. Philadelphia: Society for Industrial and Applied Mathematics. 2005,
Vese LA, Osher SJ: Modeling Textures with Total Variation Minimization and Oscillating Patterns in Image Processing. J Sci Comput. 2003, 19 (1-3): 553-572. 10.1023/A:1025384832106.
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).
YX designed and implemented the algorithms. CK and JFK performed experiments under the guidance of DNR and PND, who provided the experimental materials and platforms, participated in discussions, and contributed to the manuscript. PAI conceived the study, and participated in its implementation and coordination. YX and PAI wrote the manuscript which was read and approved by all the authors.