Translator Disclaimer
1 September 2008 Velocity distributions of single F-actin trajectories from a fluorescence image series using trajectory reconstruction and optical flow mapping
Author Affiliations +
We present an approach for the computation of single-object velocity statistics in a noisy fluorescence image series. The algorithm is applied to molecular imaging data from an in vitro actin-myosin motility assay. We compare the relative efficiency of wavelet and curvelet transform denoising in terms of noise reduction and object restoration. It is shown that while both algorithms reduce background noise efficiently, curvelet denoising restores the curved edges of actin filaments more reliably. Noncrossing spatiotemporal actin trajectories are unambiguously identified using a novel segmentation scheme that locally combines the information of 2-D and 3-D segmentation. Finally, the optical flow vector field for the image sequence is computed via the 3-D structure tensor and mapped to the segmented trajectories. Using single-trajectory statistics, the global velocity distribution extracted from an image sequence is decomposed into the contributions of individual trajectories. The technique is further used to analyze the distribution of the x and y components of the velocity vectors separately, and it is shown that directed actin motion is found in myosin extracts from single skeletal muscle fibers. The presented approach may prove helpful to identify actin filament subpopulations and to analyze actin-myosin interaction kinetics under biochemical regulation.



Advances in molecular labeling, microscopy, and image acquisition techniques offer new approaches for the analysis of molecular interactions on micrometer and submicrometer scales. These techniques have been shown to be particularly useful for real-time optical measurements and allow the visualization and quantitative analysis of biomedically relevant molecular processes. For instance, the interaction of the motor proteins actin and myosin ultimately determines the contractile properties of skeletal and cardiac muscle; therefore, a detailed understanding of this process is a crucial component in the study of contractile tissues. The highly complex in-situ situation in the presence of a large number of interacting regulator proteins and small diffusible regulator molecules, e.g., Ca2+ and ATP, possibly at nonconstant concentrations, can be abstracted experimentally. A well-studied example is represented by the actin-myosin in vitro motility assay where actin-myosin interactions are studied in a setting under controlled conditions.1 In the basic assay setup, tetramethylrhodamin-labeled actin filaments slid over a myosin-coated glass surface and the process was visualized with fluorescence microscopy. In further steps, the assay could be extended to include regulator proteins, e.g., troponin and tropomyosin, and physiological conditions could be varied.2 Typically, the quantity of interest is the distribution of filament velocities since these represent the kinetics of the underlying interaction of the actin binding head domain of the myosin molecule with its corresponding binding site on the F-actin filament.

During the image-acquisition process, some well-known artifacts and limitations are introduced in the data sets: (1) blurring by convolution of the object shape with the microscope point spread function (PSF), (2) intensity fluctuations due to the statistical nature of the fluorescence process; and (3) detector noise. Photon counting noise can be partially compensated for by longer pixel dwell times, but in practice, this is limited by the restrictions imposed by the time scale of the observed process. Moreover, computation of the sliding velocities requires the sampling rate to be sufficiently high to minimize the distance traveled by a filament between two frames. Therefore, recorded image sequences tend to be large, and differences between experimental results can be subtle. To process sufficiently large data sets for a reliable statistical analysis, data analysis procedures that are largely automated and reproducible are desirable. Basically, two conceptually different strategies have been presented: tracking involves the detection of actin filaments against the background fluorescence signal and identification of the filament trajectory by one of several methods.3 With this strategy, sliding velocities are calculated from the filament centroid dynamics.4 However, centroid tracking leads to biased velocity estimates, especially in the case of actin filaments moving along curved trajectories. The method can be refined by fitting a model shape, e.g., a Gaussian.5 The velocity estimate is then obtained from the temporal dynamics of the model parameters.

With the second strategy, the structure tensor method, filament velocities are computed from the optical flow vector field. Here, the solution of the eigenproblem for the structure tensor must be solved locally.6 The structure tensor approach avoids the segmentation and tracking problem, and returns the velocity distribution for the entire image sequence. Theoretically, the analytical methods, i.e., model-based fitting and the structure tensor method, can achieve subpixel (nm) accuracy.5

Here, we aim to combine the advantageous properties of both methods. In particular, we first demonstrate an efficient method for data preprocessing to partially invert the image formation process. The method is based on the discrete curvelet transform and emphasizes the role of filament shape restoration and computational efficiency. We then present a simple way to detect actin filaments in denoised image sequences and an approach to uniquely label noncrossing filament trajectories in space and time by combining the information of 2-D (XY) and 3-D (XYT) segmentation. Finally, it is shown that the velocity estimates of the noise-robust structure tensor approach can be mapped to individual trajectories, so, the limitations of the centroid tracking method are avoided. Preliminary results of this work were presented at the European Conference on Biomedical Optics.7




Experimental Data

Rabbit skeletal muscle actin, myosin, and heavy meromyosin (HMM) were obtained as previously described.1 Briefly, actin was labeled with rhodamine phalloidin (R-415, Molecular Probes, Eugene, Oregon) and motility assay flow cells with an approximate volume of 12μl were constructed on top of glass cover slips coated with nitrocellulose. Some of the experiments were carried out in the presence of the sNC-2 fragment of the regulating protein myosin binding protein-C (MyBP-C). In these assays, the MyBP-C peptide was added to the myosin surface at a concentration of 110μgml . In one of the presented experiments, the motility assay was prepared from a single skeletal muscle fiber using the method of Höök and Larsson.8 Myosin was extracted from a glycerinated gastrocnemius muscle of an adult C57/SV129 mouse. Experiments were carried out according to the guidelines of the local animal care committee of the University of Heidelberg.

Fluorescence excitation was achieved with a 200-W mercury arc lamp, and flow cells were observed using a 1.4nA 100× oil immersion objective (Leica PL APO) on an inverted microscope (DM IRBE, Leica Microsystems, Mannheim, Germany). Fluorescence signals were collected at 10framessecond with a CCD camera (PCO sensicam, PCO, Kelheim, Germany), digitized with 12-bit resolution, and stored as image sequences in TIFF format.

The described algorithms were carried out on image sequences with 100 frames, each 256×256pixels (pixel size dx=0.12μm , dy=0.12μm ; frame interval dt=0.1s ).


Image Processing


Multiresolution transforms

The multiresolution (MR) transform of an image distributes the complete image information on several resolution levels or scales. A sequence of low-resolution scales is obtained by applying a set of low-pass filters to the original image and storing the complementary information (“details”) separately. If the transform fulfills the perfect reconstruction conditions, the original data can be reconstructed from the lowest resolution level and the detail levels. Moreover, noise estimation and reduction on the detail levels can be carried out before reconstruction, so, an approximation of the original image can be obtained. Here we evaluate two different implementations and their respective abilities to denoise fluorescence images and restore the contained objects in a noniterative manner.


2-D wavelet transform

The 2-D discrete wavelet transform was implemented with the “à trous” algorithm based on the cubic B-spline scaling function.9, 10 In this transform, the low-resolution representation Sk,lj of an image at scale 0<j<J and pixel indices k,l is defined recursively: Sk,lj=m,nhm,nSk,lj1 , where (k,l)=(k+2j1m,l+2j1n) and hm,n represent a 2-D binomial filter. The filter hm,n is obtained from the B-spline-based filter B=116[1,4,6,4,1] as hm,n=BB . The detail or wavelet scales at scale j> 0 are computed as Wk,lj=Sk,lj1Sk,lj . For completeness, Sk,l0 is identified with the original image, and the transform is given by Sk,l0{Wk,l1,,Wk,lJ,,Sk,lJ} . Perfect reconstruction is achieved by Sk,l0=Sk,lJ+j=1JWk,lj .

This wavelet transform is highly efficient for the localization and representation of pointlike singularities. It has been observed, however, that 1-D singularities, e.g., curved edges, are not sufficiently well approximated by 2-D wavelet basis functions.11 As detailed in Sec. 2.2.3 recent advances in the design of wavelet related transforms have succeeded in achieving efficient approximations of these shapes.


Curvelet transform

To get a good approximation of curved edges or filament-like structures, the basis functions of the corresponding transform should be able to represent arbitrarily oriented ridges. This is achieved by the curvelet transform. First, let ψ() represent a 1-D wavelet. Then the 2-D function given by ψ(xnr) with x,nR2 , n=[cos(θ),sin(θ)] , rR , θ[0,2π[ is called a ridgelet. The ridgelet function is constant along lines with normal vector n=[cos(θ),sin(θ)] and (perpendicular) distance r to the origin. The ridgelet profile along the direction of the normal vector n reproduces the shape of the 1-D wavelet function ψ() . Varying n and r allows the construction of a set of functions for the reconstruction of extended straight lines and edges. To include curves and shapes of limited spatial extensions, the ridgelet transform is applied to a block decomposition of the image, as described elsewhere.11 In practice, the ridgelet transform is computed by applying a 1-D wavelet transform to the discrete radon transform of a single-image block.11 Here, we use the 1-D analogue of the wavelet transform described in Sec. 2.2.2.



For MR transform denoising, the standard soft- and hard-thresholding schemes were applied for the detail coefficients.12, 13 Let Wk,lj denote the detail coefficients at scale j obtained from the 2-D MR transform, and let σ̃j=0.67451×median(Wk,l1)×nj be the robust estimator of the noise standard deviation (SD) at wavelet scale j , where nj is the SD of standard, normally distributed white noise at scale j . The scale-adapted threshold θj=τ×σ̃j is used for (1) hard-thresholding: Wk,lJWk,lj×{Wk,lj> θj} ; or (2) soft-thresholding: Wk,lJ[Wk,ljsgn(Wk,lj)×θj]×{Wk,lj> θj} , where τ is a denoising parameter provided by the user (usually 2 to 5). Boolean expressions indicated by parentheses { } are evaluated at each pixel (k,l) , with a value of 1 if the contained relation is true, and zero otherwise. Finally, the denoised image Dk,l is computed as the sum of the coarsest scale Sk,lJ and the thresholded detail coefficients Wk,lj : Dk,l=Sk,lJ+JWk,lj .


Filament Detection

Actin filaments are detected by thresholding the denoised images. The threshold value is computed from the background intensity statistics of the denoised image as μ+3×SD with mean value μ and standard deviation SD . Detected filament structures are stored as a binary-valued image sequence It (t=0Nframes1) for subsequent spatiotemporal segmentation. The signal-to-noise ratio (SNR) of denoised images is defined using these binary images. “Signal pixels” (pixel value=1 in the binary image) in the denoised image yield the average signal pixel value μS , and the background (noise) pixels ( value=0 in the binary image) yield the noise parameters μN and σN (mean, SD). Finally, the image SNR is defined as (μSμN)×σN1 .



Figure 1 shows representative images of the in vitro actin-myosin motility assay. The segmentation of individual trajectories is complicated by several trajectory crossings, as indicated by the arrows. In three dimensions (space+time) , these trajectories form widely connected structures.

Fig. 1

In vitro actin-myosin motility assay. The image sequence shows snapshots of fluorescently labeled actin filaments moving over a heavy meromyosin (HMM) coated surface. Segmentation of single-filament trajectories is complicated by multiple crossings of actin filaments (indicated by arrows), leading to widely connected spatiotemporal structures. Furthermore, actin filaments often change their shape and direction of motion.


The 3-D segmentation procedure we developed here is applied to the binary-valued image sequences obtained by thresholding, as described above, and consists of three steps: (1) initial labeling of 3-D connected structures, (2) detection of trajectory crossings, entries, and exits, and (3) final relabeling of individual filament trajectories.


Initial labeling

For segmentation, we consider the “six neighborhood” in space and time, i.e., the neighborhood of each pixel (i,j,k) in the image sequence is defined as the set of pixels {(i±1,j,k),(i,j±1,k),(i,j,k±1)} , while border pixels are ignored. In two dimensions, the corresponding “four neighborhood” is used. The initial labeling identifies 3-D connected structures in the binary-valued sequence as sets of positively valued pixels for which each member is connected to at least one other member. Due to multiple crossings, these 3-D structures often contain the trajectories of numerous individual actin filaments.


Detection of trajectory crossings

We developed the following algorithm to detect trajectory crossings, entries, and exits: Consider two subsequent image frames where a putative crossing takes place. It is observed that, within these frames, noncrossing trajectories will yield multiple 2-D segments in each frame and multiple segments when 3-D segmentation is applied. A trajectory crossing, however, connects 2-D isolated structures in time and therefore yields a single 3-D (XYT) connected structure. This situation is illustrated in Fig. 2 , where a typical filament crossing in a reduced two-image stack is shown. Figure 2a shows two consecutive frames of an image sequence. In the lower half of the images, the beginning of a trajectory crossing can be observed. After 2-D segmentation, binary-valued images are obtained, and filaments appear as gray structures in Figs. 2b and 2c. The 2-D isolated regions R1 and R2 in It (each belongs to a different actin filament) merge to form the 2-D isolated region R3 in frame It+1 , while the complete set of regions R1,2,3 forms a single 3-D connected structure. Noncrossing trajectories yield isolated 3-D structures, as in the case of the remaining trajectories in Fig. 2. We will call a 3-D connected structure that extends over more than one frame an “extended segment.” Whether a 3-D segment is extended or not is easily checked by looking at the set of time coordinates of the pixels forming the segment. In the presence of multiple filaments, and given the possibility of crossings of more than two filaments, crossings must be identified by comparison with all possible 2-D segment pairs. Furthermore, new trajectories entering the imaging field and trajectories that exit must be identified. Therefore, for each 2-D region Ri in frame It , we check whether it has a connected region, a “successor,” in the subsequent frame It+1 as follows: we perform a 3-D region labeling of a reduced two-frame stack that consists of only region Ri in one frame and It+1 in a second frame. If the result shows no extended objects, Ri has no successor and the region is classified as a trajectory entry/exit. A formal description of the algorithm is given by:

  • 1. For each of the frames It (t=0Nframes2) of the binary-valued image sequence, perform 2-D region labeling to yield a set of 2-D connected regions {R1,R2,,Rn}t .

  • 2. Check if Ri or Rj represents a trajectory entry or exit. If so, skip the pair and proceed with the next pair of regions.

  • 3. For each pair of regions {Ri,Rj}t , get a new frame F0 that contains only regions Ri and Rj , and compose a reduced two-frame stack (F0,F1) , F1=It+1 .

  • 4. Perform a 3-D region labeling of the reduced stack.

  • 5. Determine the number Ne of extended 3-D segments. If Ne=2 , regions Ri and Rj do not merge, and thus, no trajectory crossing occurs. If Ne=1 , a crossing is detected and the corresponding segment in frame It+1 will be set to zero.

  • 6. Run the algorithm backwards in time to detect “splitting” trajectories.

Given the experimental details described in Sec. 2, this part of the algorithm has an average execution time of less than 30seconds .

Fig. 2

(a, b) In these two consecutive frames of an image sequence, the crossing of two actin filaments is shown (lower half) while all other filaments are moving freely. (c) 2-D segmentation: filament detection in all frames yields a binary-valued sequence. Region labeling of 3-D connected structures provides an initial trajectory segmentation, including all trajectory crossing events. The crossing of filament regions R1 and R2 (in frame It ) produces a single connected region R3 in frame It+1 . All noncrossing trajectories, in contrast, show a one-to-one mapping between the frames. This feature is detected by testing 3-D connectivity in a reduced stack for all pairs of regions (for details, see the text).



Individual filament trajectories

Application of the aforementioned procedure detects single-frame regions where filament trajectories cross and split. Due to the length of the actin filaments used here, and given all possible crossing angles, crossing events often include more than one but usually not more than five frames. Since the kind of molecular interaction in these regions is not well defined, i.e., an isolated actin-myosin interaction cannot be assumed, these regions are excluded from the velocity statistics. Thus, the obtained velocity measures give a more valid representation of the process under observation. So far, all trajectory merging and splitting points have been deleted. To also exclude the remaining short segments to which more than one actin filament contribute, and where no well-defined actin-myosin interaction takes place, we delete all trajectories shorter than five frames. The resulting binary image sequence is relabeled to obtain unique integer labels for individual filament trajectories that can be mapped to the optical flow vector field.


Optical flow and velocity computation

The optical flow vector field and the corresponding velocity measures were computed via the 3-D structure tensor approach as previously described.6 Briefly, at each pixel location, the fluorescence intensity gradient vector is obtained by convolution of the image sequence with a 3-D directional derivative filter. For efficient computation, each of these 3D filters is linearly separated in a set of three 1-D filters. Here, the rotation-invariance-optimized 5-tab Sobel and binomial filters are used.14 Let D and B denote 1-D Sobel and binomial filters, respectively, and let {i,j,k} denote any permutation of the three data coordinates {x,y,t} . The directional derivative tI of the image sequence I in direction i is then given by convolution (denoted by *) of I with three 1-D filters iI=Di*Bj*Bk*I , i.e., derivation in the desired direction i and smoothing in the remaining directions j and k .14 At each pixel, the three directional derivatives are combined to yield the gradient vector I=(xI,yI,tI)T . As the gradient is computed at each pixel, each of the three gradient components xI,yI,tI represents a data set of the same size as the original image sequence. To obtain a local measure of 3-D oriented structures in the image sequence, the structure tensor J is computed by averaging the dyadic product of the gradient vector with itself in a neighborhood of each pixel: J=G*(II) .14 Here, G represents a Gaussian averaging kernel implemented by a discrete 7-tab binomial filter mask. The 3×3 structure tensor contains six independent components Jxx,Jyy,Jtt,Jxy,Jxt,Jyt , where each component Jkl is the locally averaged, pointwise product of two directional derivatives Jkl=G*(kIlI) . Thus, a structure tensor can be constructed for each pixel in the image sequence. To implement Gaussian averaging (the filter G ) for each of these six independent tensor components, a linearly separable convolution with three 1-D binomial filters must be carried out.

Finally, the 3×3 structure tensor J is obtained and diagonalized at each pixel. From the standard eigen-decomposition J=MΛM1 , three real eigenvalues λl=Λll and the corresponding eigenvectors vl=(M1l,M2l,M3l)T are extracted. After reordering the eigenvalues λ1λ2λ3 and their corresponding eigenvectors, the optical flow vector f=(fx,fy) is computed from the eigenvector v3=(M13,M23,M33) to the smallest eigenvalue λ3 as follows: f=(M13M33,M23M33) . Actually, the optical flow vector is only calculated if the rank of the structure tensor equals 2 to avoid the well-known aperture problem.15 In practice, this is evaluated using the quantities “coherency” (coh) and “corner” (cor), derived elsewhere,14 and defined using the ordered eigenvalues, as coh=(λ1λ3λ1+λ3)2 , cor=coh(λ1λ2λ1+λ2)2 . The optical flow is computed and stored for further analysis only if (coh> 0.9)(cor> 0.8) . The use of these measures provides a noise-robust velocity estimation and allows the application of the algorithm to the raw fluorescence sequences. The x and y components of the filament velocity vectors are stored as image sequences to map their locations to the segmented trajectories. Because a structure tensor is defined at each pixel, the number of optical flow vectors obtained could, theoretically, equal the number of pixels in the image sequence.

Due to the aperture problem, which is always present in optical flow estimation, valid estimates are computed only for the subset of pixels that fulfills the conditions defined above. However, the number of velocity estimates for our image sequences (256×256×100) is in the range of 103 to 105 . It is important to note that the number of velocity estimates does not equal the number of filaments, but instead, the number of pixels where the optical flow vector can be estimated reliably. Usually this algorithm is applied to fluorescence image sequences without any knowledge about filament localizations and provides a velocity distribution related to the entire image sequence. Using this conventional algorithm, individual velocities cannot be mapped to filaments—a mapping is achieved only if filaments and their individual trajectories are found by an additional segmentation procedure, as presented above.

Implementation of the discrete curvelet transform is based on previously published algorithms.16 Image processing code was written in Java and IDL6.0 (RSI, Boulder, Colorado). Image sequence analysis was carried out on a PC, Intel Dual Core, 2×2.13GHz , 2GB RAM.


Results and Discussion


Image Sequence Denoising

We evaluated the denoising properties of the 2-D discrete wavelet transform and the 2-D discrete curvelet transform. Both algorithms were applied to each individual frame of 10 different input sequences ( 256px×256px×100 frames). The amount of noise reduction was quantified with the SNR as defined above. The average SNR of the raw image sequences was 3 . While curvelet denoising increased the SNR to almost 4 ( 3.98±0.14 , mean±SEM , n=10 ), the wavelet-based algorithm led to an SNR increase to more than 6 ( 6.33±0.26 , mean±SEM , n=10 ). This result may suggest the conclusion that wavelet denoising is more efficient for noise reduction, but we observed that, although both transforms significantly reduced noise, their respective results in terms of denoised filament shapes also differed in an important way. This difference is illustrated in Fig. 3 . In Fig. 3a, a detail of a representative raw image with several fluorescent actin filaments is shown. The white line indicates a profile along the axis of a single filament, and the measured fluorescence intensity on this line is shown below [Fig. 3d]. The 1-D intensity profile [Fig. 3d] illustrates why the application of a simple threshold operation cannot be carried out and often leads to fragmentation of the filament contour. Figures 3b, 3e, 3h show the results obtained with wavelet denoising (soft-thresholding, J=4 , τ=2.5 ). The background noise level was significantly reduced, and thus, filament detection was facilitated. As observed in Fig. 3e, however, the ragged intensity profile along the filament “ridge” was reproduced by the transform, though with less noise. Yet, due to the marked decrease in background intensity SD, a detection threshold can be chosen with less risk of filament fragmentation compared to the case of raw images. In the right column [Figs. 3c, 3e, 3i], a representative result obtained with the curvelet-based denoising algorithm is shown (hard-thresholding, J=4 , τ=5 ). Although background fluctuations are not so effectively suppressed compared with the wavelet transform, clearly the filament shape is reconstructed more reliably by curvelet basis functions. The 1-D plot [Fig. 3f] shows that the intensity profile on the filament ridge is less ragged and has a larger amplitude than the profile of the wavelet denoised image in Fig. 3e. For the wavelet-based procedure, a further increase in the denoising parameter τ , a change in the number of decomposition levels J , or the use of hard-thresholding did not change the qualitative characteristics of these results (not shown). In total, both transforms facilitated filament detection, but the curvelet-based approach showed a better performance in the case of curved, extended filament shapes. The relative effect of both algorithms on background noise and filament restoration is visualized in Figs. 3g, 3h, 3i. All 3-D surface representations use the same color table as shown in Fig. 3c and represent a region around the white line drawn in Fig. 3a. Based on these results, we subsequently employed the curvelet-based denoising algorithm for image sequence preprocessing prior to the following analysis steps.

Fig. 3

Denoising and object restoration: (a) Detail of an unprocessed image with several actin filaments. Profile plots shown in (d), (e), and (f) are measured along the white line indicated in panel (a). (b) Wavelet denoised image showing a strong reduction of background intensity fluctuations. (c) Curvelet denoising: background noise is not as efficiently reduced as in (b), but the elongated filament shape is restored more reliably. (d)–(f) Intensity profiles measured along the line shown in (a). Wavelet denoising (e) reduces background noise but also enhances the local extrema of the filament shape, while curvelet denoising (f) stabilizes the intensity profile along the filament ridge. Note the different scalings of the y axes (intensity). (g)–(i) Combined effects of background noise reduction and filament shape restoration as a 3-D surface representation [same color table as in (a)–(c)]. (Color online only).




Curvelet denoised image sequences were binarized by applying a threshold, and an initial integer labeling of connected regions was carried out. The volume of these trajectories relative to the total pixel volume was 1.4% ( 92261±13443px , mean±SD ). The initial number of trajectories showed a wide range of values (19 to 390) due to the fact that small regions in single frames were also recognized as potential trajectories. These regions often comprised fewer than five pixels and could have been removed easily by morphological operators (e.g., erosion). Our strategy here, however, was to negatively select these regions by a temporal criterium, i.e., trajectories had to span at least five frames. Since this step is already part of the crossing trajectory detection (as detailed in Sec. 2), no further computational load was imposed, and the detected filament shapes were not further influenced by any operators. After detection and deletion of trajectory crossing points, the average trajectory volume was reduced to 1.1% of the total pixel volume ( 72474±10917px , mean±SD ). The amount of reduction depended on the density of the traveling filaments, because this determines the frequency of trajectory crossings. This analysis also shows that, given the conditions used here, a significant amount of the total trajectory volume (20%) consisted of actin-actin crossings, a situation of rather unclear molecular interactions. After complete segmentation and relabeling, the final number of trajectories (range of 26 to 78) was more uniform than the initial counting because noise-related detections were reliably deleted. The maximum value of this range was produced by a single experiment with a high density of actin filaments. Ignoring this experiment, the range was 26 to 42 trajectories per image sequence. It should be noted that the final number of trajectories generally tends to be lower than the initial counting; however, in the case of a few long trajectories with several crossings, the final trajectory count can exceed the initial number. Likewise, the trajectory count cannot be interpreted as the actual number of different actin filaments present in the experiment—an important fact to remember for statistical analyses. Since filaments also leave and re-enter the imaging field, the trajectory count only provides a coarse estimate of the actual filament count.


Optical Flow Mapping

For each image sequence, the optical flow vector field was computed, and using the appropriate scaling factors, the local velocity vectors were determined. Velocity vectors were obtained only at locations where the optical flow computation was reliable, as detailed in Sec. 2. Taking the intersection of the pixel sets where (1) a uniquely labeled actin filament trajectory existed, and (2) a velocity vector was obtained, we finally obtained a mapping of velocity vectors to trajectories and compared their respective distributions. Figure 4a shows a 3-D representation of the segmented filament trajectories of a single experiment. The velocity of a trajectory segment corresponds to the local slope of the trajectory with respect to the X-Y plane. Different slopes (absolute velocities) and directions of motion can be observed. In particular, straight motion, sharp turns, and even spiralling movements can be seen. Figure 4b shows several single-trajectory velocity distributions. For better visibility, only half of the trajectory distributions (even integer labels) are shown, and all histograms are shown as cubic B-spline interpolated continuous curves. By analyzing the respective contributions of different trajectories, it can be observed that some trajectories display a broad velocity distribution while others have a sharp peak and contribute solely to a narrow velocity band. For a specific distribution of interest, the velocity values can be mapped back to their respective spatial and temporal positions. Previously used methods pooled all velocity values in a global velocity histogram, as shown in Fig. 4c. By comparing Figs. 4b and 4c, it can be seen that the contributions of single trajectories to the global histogram cannot be anticipated.

Fig. 4

(a) 3-D isosurface reconstruction of detected spatiotemporal actin filament trajectories. Filament motions are complex and include straight lines, sharp turns, and spiraling movements. (b) After the mapping of the optical flow vector field derived from the local structure tensor, each trajectory contributes a single velocity distribution. For better visualization, only half of these distributions are shown. Some trajectories contribute only to a single peak in the global velocity distribution (c), while others have a broad distribution of velocities. Distributions are shown as cubic B-spline interpolated histograms. (c) Pooled distribution of absolute velocities and the underlying histogram (bars).



Uniform Motion and Histogram Reconstruction

It is a well-known phenomenon that actin filaments display nonuniform patterns of motion, and eventually they can get stuck due to imperfections in the HMM surface and to the formation of strong binding states. This can lead to either complete or partial immobility of these filaments. Partial immobility can occur in at least two ways: first, a part of the filament may remain tightly bound to the underlying myosin while the remaining part moves erratically around this fixed point, second, an evenly moving filament may suddenly get stuck and either remain immobile or resume its trajectory (a “stop and go” pattern). In any case, these filaments contribute erroneous and mostly low velocities to the histogram. Figure 5 shows velocity time courses of two filament trajectories. Each data point is obtained by averaging the absolute velocities corresponding to a filament in a given image frame. Since there are numerous optical flow vectors for each structure, we present the velocity in a given frame as mean±SD values. The trajectory denoted as “moving filament” (grey curve) represents an evenly moving actin molecule with a quasi-constant sliding velocity. However, the other filament (stop and go, black curve) does not show significant motion for almost 30 frames, except for a small displacement near frame 15, before it eventually starts to move with the typical average sliding velocity of 6μms , interrupted by another short deceleration around frame 43. For evenly moving filaments, perturbations due to stop-and-go phenomena, including several acceleration and deceleration phases, that are caused by geometrical and biochemical inhomogeneities of the preparation can be considered minimal. Thus, to distinguish different experimental conditions, it is especially interesting to identify evenly moving objects.

Fig. 5

By tracking the mean velocity of individual filaments across the image sequence, velocity time courses are obtained. Here, some filaments show a relatively constant sliding velocity (moving filament, grey trace), while others get stuck due to imperfections in the myosin surface and the formation of strong binding states. In that case, filaments can display a “stop and go” motion pattern (black trace). Inclusion of these trajectories can lead to broad velocity distributions, as shown in Fig. 4b and can eventually distort the pooled velocity statistics.


Different ways of identifying evenly moving filaments can be considered. We propose a straightforward way that uses the SD of the time series formed by the mean velocities (filled circles in Fig. 5). Evenly moving objects are identified by velocity values with a low SD, here SD<1.5μms . Figure 6 illustrates the difference between the velocity distribution, including all found trajectories [Fig. 6a], and the distribution computed from the same data set, but includes only evenly moving filaments with a narrow velocity profile as defined by the SD [Fig. 6b]. While both distributions show the same main peak at 6μms , the full distribution in Fig. 6a has a broad tail extending to the left with two smaller peaks at lower velocities. An examination of both distributions shows that these peaks correspond to trajectories with broad velocity profiles and not to evenly moving filaments traveling at a constant but lower velocity. Because the selection criterion is based on the shape of the velocity distribution, there is no bias against high or low velocities.

Fig. 6

Pooled velocity distributions with and without the selection of quasi-constant sliding velocities. The global velocity distribution of the same motility assay experiment computed using two selection criteria. (a) All individual filament trajectories are included in the analysis. The resulting velocity distribution shows a peak around 6μms , but also a broad tail toward slower velocities with additional peaks. (b) Selecting trajectories with quasi-constant sliding velocity (SD<1.5μms) , we derive a narrow velocity distribution with a clear peak at 6μms . The smaller peaks in panel (a) were derived largely from contributions of trajectories with a broad velocity distribution rather than from constantly moving slow filaments.


Motility assays are often used to investigate the influence of regulator proteins on motor protein interactions by reconstituting actin and/or myosin with these regulator molecules. Here, we present experimental data obtained from a motility assay reconstituted with the sNC-2 fragment of the myosin binding protein-C (MyBP-C). Figure 7 shows velocity distributions obtained under both experimental conditions, i.e., from assays in the presence of 110μgml sNC-2 (grey bars and curves), and control assays without the regulator (black bars and curves). Five data sets were analyzed for each condition, each corresponding to a different flow cell. Figure 7a shows the distributions obtained by pooling the velocity values of all detected filament trajectories, and Fig. 7b represents the distributions obtained by selecting only evenly moving filaments. Due to the different number of values obtained under both conditions, area-normalized distributions are shown. The analysis of filaments with a quasi-constant sliding velocity leads to a more pronounced separation of both distributions. The main velocity peak of the control experiment ( 6μms , black curve) remains unaltered by the procedure, while the tendency toward slower sliding velocities under the influence of sNC-2 is enhanced. The 6μms peak disappears under the action of sNC-2 without the appearance of another clear peak velocity, thus, the resulting distribution in Fig. 7b resembles a monotonically decaying curve. Although it is known from other studies that MyBP-C reduces the mean sliding velocity of actin by inhibiting the actomyosin ATPase activity,17 the resulting velocity distributions have not been analyzed yet. Here, we show that while isolated actin-HMM interaction produces a preferred sliding velocity of 6μms , an additional interaction with the sNC-2 fragment of MyBP-C leads to a decaying, broad velocity distribution—which, however, is composed of relatively narrow sliding velocity distributions for individual filaments (SD<1.5μms) . Methodologically, this result can be obtained only by mapping sliding velocities to individual filaments and their trajectories.

Fig. 7

Influence of the sNC-2 fragment of the regulating protein MyBP-C on actin sliding velocity in vitro. Global velocity distributions are derived from 5 flow cells for each condition, i.e., in the presence (grey) or absence (control, grey) of MyBP-C. (a) Including all filament trajectories, two broad and partially overlapping velocity distributions are derived. (b) By selecting quasi-constantly moving filaments (SD<1.5μms) , both distributions separate more clearly. In the control case, the main velocity peak at 6μms is pronounced, while evenly moving filaments in the regulated assay are grouped around lower sliding velocities due to the inhibiting effect of MyBP-C on actomyosin-ATPase kinetics.



Directed Motion

To reliably mimic the intracellular geometry, the orientational coherency of the myosin surface and the oriented motion of actin filaments have become important aspects in advanced motility assay designs. Using the method presented here, we obtained a measure of directionality of motion for a given experimental setup and also for each filament. Actin filament trajectories from our experiments obtained with a modified motility assay are analyzed in terms of their directional components. In this setup, the in vitro actin-myosin motility assay was not carried out using randomly oriented HMM arrays on a surface, but using instead myosin extracted directly from single skeletal muscle fibers. This latter preparation allegedly yields a pattern of oriented myosin heads, and as a consequence, induces directed motion of actin molecules. We used our method to test the directed movement of actin filaments in this assay. The results are presented in Fig. 8 , where the 2-D velocity distributions of two representative experiments are shown as surface plots. Figure 8a shows the cubic B-spline interpolated velocity histogram of actin filaments moving over an HMM-coated surface. Other than the central peak at zero velocity, an almost circular distribution of x and y velocities is observed. Results obtained with the fiber extraction method, as shown in Fig. 8b, yields a highly elongated distribution of velocity vectors that have an almost vanishing x component, and y velocities that peak at approximately ±12μms . These results quantify directed filament motions in this kind of assay.

Fig. 8

2-D velocity distributions showing x and y components of the velocity vector. (a) In the HMM-based motility assay, filaments move in all directions as shown in the almost circular distribution of the velocity vector distribution. (b) When myosin is extracted from single skeletal muscle fibers, filaments move almost exclusively in the y direction.




Molecular assays investigating molecular motion and interaction of motor proteins are important tools in biomedical research. For several reasons, quantitative and semiautomated methods for data analysis are of increasing importance: large data sets can be obtained in shorter periods of time thanks to advances in modern imaging technology; even subtle differences in results obtained under different experimental conditions must be reliably estimated; and new experimental designs must be evaluated.

Image preprocessing for noise reduction and object restoration can significantly increase the performance of these methods. The approach presented here is based on the discrete curvelet transform and offers several advantages over the previously presented algorithm based on anisotropic diffusion filtering.18 The curvelet-based algorithm is a one-step procedure with a faster execution time and has fewer parameters that need adjustment (e.g., number of iterations, integration time interval, parameters for the edge-stopping function). On the other hand, an optimized anisotropic diffusion filtering also provides excellent results for this application. Instead of tracking single particles, we performed object detection in all frames independently, followed by a novel segmentation scheme working on the spatiotemporal 3-D trajectories. Using structure tensor-based optical flow estimation, we avoided problems such as bias against curved and circular motions (centroid tracking) or difficulties with rapidly changing object shapes (cross-correlation tracking). The combination of these methods yielded a highly user-independent algorithm and allowed a more detailed analysis of experimental results, because all contributions of single trajectories to the global velocity distribution were readily obtained.

To analyze functional relationships in the complex molecular assemblies that motor proteins form with numerous regulating proteins, the in vitro motility assay represents a flexible experimental platform where situations of variable complexity can be studied. Thus, in the case of muscle contraction, the basic constituents actin and myosin can successively be reconstituted with their molecular interaction partners, e.g., myosin binding proteins, tropomyosin and troponin, to approximate the in situ situation. We presented data obtained in the presence of the sNC-2 fragment of MyBP-C. MyBP-C binds to myosin and to titin and is an important regulator of sarcomeric structure and contractility in skeletal and cardiac muscle. In human heart muscle, mutations are associated with hypertrophic cardiomyopathy.19 MyBP-C’s inhibiting function on actomyosin-ATPase activity has been described as a molecular brake on the actin-myosin motor.19 With our approach, we analysed the influence of the sNC-2 peptide on cross-bridge kinetics in the case of individual actin filaments. We were able to show that the sNC-2 peptide itself exerts a “brake” function on the skeletal muscle-derived actin-myosin system, and the effect was produced in the presence of only these three proteins. Thus, neither titin nor other regulating proteins are necessary to mediate the observed effect. In particular, the peaked velocity distribution of sliding actin filaments in the absence of sNC-2 was replaced by a monotonously decaying distribution, and we know that in both cases the data were derived from evenly moving filaments with only a quasi-constant sliding velocity. Using this approach, it will be interesting to study other MyBP-C peptides or mutation products to assess their ability to influence actomyosin-ATPase activity in skeletal and cardiac muscle, and to relate these molecular-scale in vitro results to observed phenotypes, e.g., the occurrence of cardiomyopathies. Likewise, the action of the actin-binding proteins tropomyosin and troponin on acto-myosin kinetics offers an exciting perspective to better understand muscle contractility. Finally, it should be noted that in vitro protein reconstitution can lead to a mixed population of reconstituted and nonreconstituted molecules. Only by analyzing individual filament trajectories can the subpopulations be separated for statistical analysis.

The presented algorithm was designed especially for the application to motility assay data, as described. This experimental setup presented some special features not present in all situations where tracking and motion determination are carried out. First, we worked with a 2-D geometry due to the surface-bound myosin molecules. Thus, our image sequences did not represent 3-D to 2-D intensity projections, as is the case with microscopy images taken from whole cells where vesicles or fluorescently labeled proteins are routinely tracked. In the latter case, the deletion of crossing points cannot be carried out using the reasoning presented here, because any trajectory crossing in a given image plane could result from projection. For these cases, interesting new methods that take into account the geometry of the studied cell have been developed.20 Second, actin filaments are extended, flexible structures with a length of up to 8μm . The resulting shape variations provoke computational challenges for the tracking and identification of an individual filament across images. In the future, a model-based analysis of actin motility that is based on structural and biochemical data and takes into account the physical properties of the measurement setup may spark new approaches for image sequence analysis. Similar approaches applied for microtubule dynamics have already shown promising results.21 However, to obtain more information on actin-myosin interactions in vitro, the presented method may prove to be a useful contribution.


Author M. Vogel acknowledges ongoing support from the Center for Nanoscale Systems, Harvard University, Cambridge, Massachusetts, and the National Nanotechnology Infrastructure Network (NNIN; National Science Foundation Cooperative Agreement No. 0335765). Myosin binding protein-C was kindly provided by M. Gautel, King’s College, London, UK.



S. J. Kron and J. A. Spudich, “Fluorescent actin filaments move on myosin fixed to a glass surface,” Proc. Natl. Acad. Sci. U.S.A., 83 6272 –6276 (1986). 0027-8424 Google Scholar


M. Kawai, T. Kido, M. Vogel, R. H. A. Fink, and S. Ishiwata, “Temperature change does not affect force between regulated actin filaments and heavy meromyosin in single-molecule experiments,” J. Physiol. (Paris), 574 (3), 877 –887 (2006). 0021-7948 Google Scholar


M. K. Cheezum, W. F. Walker, and W. H. Guilford, “Quantitative comparison of algorithms for tracking single fluorescent particles,” Biophys. J., 81 2378 –2388 (2001). 0006-3495 Google Scholar


S. S. Work and D. M. Warshaw, “Computer-assisted tracking of actin filament motility,” Anal. Biochem., 202 275 –285 (1992). 0003-2697 Google Scholar


M. I. Wallace, C. Batters, L. M. Coluccio, and J. E. Molloy, “Nanometre resolution tracking of myosin-1b motility,” IEE Proc.: Nanobiotechnol., 150 134 –140 (2003). 1478-1581 Google Scholar


D. Uttenweiler, C. Veigel, R. Steubing, C. Götz, S. Mann, H. Haussecker, B. Jähne, and R. H. A. Fink, “Motion determination in actin filament fluorescence images with a spatio-temporal orientation analysis method,” Biophys. J., 78 2709 –2715 (2000). 0006-3495 Google Scholar


F. von Wegner, T. Ober, C. Weber, O. Friedrich, R. H. A. Fink, and M. Vogel, “Multiresolution transform denoising and segmentation of single molecule motility image series,” Proc. SPIE, 6626 66260Q (2007). 0277-786X Google Scholar


P. Höök and L. Larsson, “Actomyosin interactions in a novel single muscle fiber in vitro motility assay,” J. Muscle Res. Cell Motil., 21 357 –365 (2000). 0142-4319 Google Scholar


M. J. Shensa, “The discrete wavelet transform: wedding the à trous and Mallat algorithms,” IEEE Trans. Signal Process., 40 2464 –2482 (1992). 1053-587X Google Scholar


F. von Wegner, M. Both, and R. H. A. Fink, “Automated detection of elementary calcium release events using the à trous wavelet transform,” Biophys. J., 90 2151 –2163 (2006). 0006-3495 Google Scholar


J. L. Starck, D. L. Donoho, and E. J. Candès, “Astronomical image representation by the curvelet transform,” Astron. Astrophys., 398 785 –800 (2003). 0004-6361 Google Scholar


S. G. Chang, B. Yu, and M. Vetterli, “Adaptive wavelet thresholding for image denoising and compression,” IEEE Trans. Image Process., 9 1532 –1546 (2000). 1057-7149 Google Scholar


D. L. Donoho and I. M. Johnstone, “Adapting to unknown smoothness via wavelet shrinkage,” J. Am. Stat. Assoc., 90 1200 –1224 (1995). 0162-1459 Google Scholar


H. Haussecker and H. Spies, “Motion,” Handbook of Computer Vision and Applications, 2 310 –396 Academic Press, San Diego, California (1999). Google Scholar


H. Scharr and H. Spies, “Accurate optical flow in noisy image sequences using flow adpated anisotropic diffusion,” Signal Process. Image Commun., 20 537 –553 (2005). 0923-5965 Google Scholar


J.-L. Starck, Y. Moudden, P. Abrial, and M. Nguyen, “Wavelets, ridgelets and curvelets on the sphere,” Astron. Astrophys., 446 1191 –1204 (2006). 0004-6361 Google Scholar


M. V. Razumova, J. F. Shaffer, A.-Y. Tu, G. V. Flint, M. Regnier, and S. P. Harris, “Effects of the N-terminal domains of myosin binding protein-C in an in vitro motility assay,” J. Biol. Chem., 281 35846 –35854 (2006). 0021-9258 Google Scholar


D. Uttenweiler, C. Weber, B. Jähne, R. H. A. Fink, and H. Scharr, “Spatiotemporal anisotropic diffusion filtering to improve signal-to-noise ratios and object restoration in fluorescence microscopic image sequences,” J. Biomed. Opt., 8 40 –47 (2003). 1083-3668 Google Scholar


E. Flashman, C. Redwood, J. Moolmon-Smook, and H. Watkins, “Cardiac myosin binding protein C—its role in physiology and disease,” Circ. Res., 94 1279 –1289 (2004). 0009-7330 Google Scholar


S. Bonneau, M. Dahan, and L. Cohen, “Single quantum dot tracking based on perceptual grouping using minimal paths in a spatiotemporal volume,” IEEE Trans. Image Process., 14 1384 –1395 (2005). 1057-7149 Google Scholar


J. Boulanger, C. Kervrann, and P. Bouthemy, “Simulation and estimation of fluorescence microscopy image sequences for intracellular dynamics and trafficking,” 18 –26 (2006). Google Scholar
©(2008) Society of Photo-Optical Instrumentation Engineers (SPIE)
Frederic von Wegner, Tobias Ober, Cornelia Weber, Sebastian Schürmann, Rene Winter, Oliver Friedrich M.D., Rainer H. A. Fink, and Martin Vogel "Velocity distributions of single F-actin trajectories from a fluorescence image series using trajectory reconstruction and optical flow mapping," Journal of Biomedical Optics 13(5), 054018 (1 September 2008).
Published: 1 September 2008

Back to Top