2.1.

Experimental Data

Rabbit skeletal muscle actin, myosin, and heavy meromyosin (HMM) were obtained as previously described.¹ Briefly, actin was labeled with rhodamine phalloidin (R-415, Molecular Probes, Eugene, Oregon) and motility assay flow cells with an approximate volume of $12\mu \mathrm{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\mu \mathrm{g}∕\mathrm{ml}$ . 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.⁸ 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\text{-}\mathrm{W}$ mercury arc lamp, and flow cells were observed using a $1.4\mathrm{nA}$ $100\times$ oil immersion objective (Leica PL APO) on an inverted microscope (DM IRBE, Leica Microsystems, Mannheim, Germany). Fluorescence signals were collected at $10\text{frames}∕\text{second}$ with a CCD camera (PCO sensicam, PCO, Kelheim, Germany), digitized with $12\text{-}\text{bit}$ resolution, and stored as image sequences in TIFF format.

The described algorithms were carried out on image sequences with 100 frames, each $256\times 256\text{pixels}$ (pixel size $dx=0.12\mu \mathrm{m}$ , $dy=0.12\mu \mathrm{m}$ ; frame interval $dt=0.1\mathrm{s}$ ).

2.2.

Image Processing

2.3.

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 $\mu +3\times SD$ with mean value $\mu$ and standard deviation $SD$ . Detected filament structures are stored as a binary-valued image sequence $I_t$ $(t=0\dots N_{\text{frames}}-1)$ for subsequent spatiotemporal segmentation. The signal-to-noise ratio (SNR) of denoised images is defined using these binary images. “Signal pixels” (pixel $\text{value}=1$ in the binary image) in the denoised image yield the average signal pixel value ${\mu }_S$ , and the background (noise) pixels ( $\text{value}=0$ in the binary image) yield the noise parameters ${\mu }_N$ and ${\sigma }_N$ (mean, SD). Finally, the image SNR is defined as $({\mu }_S-{\mu }_N)\times {\sigma }_N^{-1}$ .

2.4.

Segmentation

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 $(\text{space}+\text{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.

2.4.2.

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 $R_1$ and $R_2$ in $I_t$ (each belongs to a different actin filament) merge to form the 2-D isolated region $R_3$ in frame $I_{t+1}$ , while the complete set of regions $R_{1,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 $R_i$ in frame $I_t$ , we check whether it has a connected region, a “successor,” in the subsequent frame $I_{t+1}$ as follows: we perform a 3-D region labeling of a reduced two-frame stack that consists of only region $R_i$ in one frame and $I_{t+1}$ in a second frame. If the result shows no extended objects, $R_i$ 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 $I_t$ $(t=0\dots N_{\text{frames}}-2)$ of the binary-valued image sequence, perform 2-D region labeling to yield a set of 2-D connected regions ${\{R_1,R_2,\dots ,R_n\}}_t$ .

• 2. Check if $R_i$ or $R_j$ 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 ${\{R_i,R_j\}}_t$ , get a new frame $F_0$ that contains only regions $R_i$ and $R_j$ , and compose a reduced two-frame stack $(F_0,F_1)$ , $F_1=I_{t+1}$ .

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

• 5. Determine the number $N_e$ of extended 3-D segments. If $N_e=2$ , regions $R_i$ and $R_j$ do not merge, and thus, no trajectory crossing occurs. If $N_e=1$ , a crossing is detected and the corresponding segment in frame $I_{t+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 $30\text{seconds}$ .

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 $R_1$ and $R_2$ (in frame $I_t$ ) produces a single connected region $R_3$ in frame $I_{t+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).

3.1.

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 ( $256\mathrm{px}\times 256\mathrm{px}\times 100$ frames). The amount of noise reduction was quantified with the SNR as defined above. The average SNR of the raw image sequences was $\sim 3$ . While curvelet denoising increased the SNR to almost 4 ( $3.98\pm 0.14$ , $\text{mean}\pm \mathrm{SEM}$ , $n=10$ ), the wavelet-based algorithm led to an SNR increase to more than 6 ( $6.33\pm 0.26$ , $\text{mean}\pm \mathrm{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$ , $\tau =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$ , $\tau =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 $\tau$ , 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).

3.2.

Segmentation

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 $\sim 1.4\%$ ( $92261\pm 13443\mathrm{px}$ , $\text{mean}\pm \mathrm{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 $\sim 1.1\%$ of the total pixel volume ( $72474\pm 10917\mathrm{px}$ , $\text{mean}\pm \mathrm{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 $(\sim 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.

3.3.

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\text{-}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).

3.4.

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 $\text{mean}\pm \mathrm{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\mu \mathrm{m}∕\mathrm{s}$ , 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\mu \mathrm{m}∕\mathrm{s}$ . 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 $\sim 6\mu \mathrm{m}∕\mathrm{s}$ , 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\mu \mathrm{m}∕\mathrm{s}$ , but also a broad tail toward slower velocities with additional peaks. (b) Selecting trajectories with quasi-constant sliding velocity $(\mathrm{SD}<1.5\mu \mathrm{m}∕\mathrm{s})$ , we derive a narrow velocity distribution with a clear peak at $6\mu \mathrm{m}∕\mathrm{s}$ . 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\mu \mathrm{g}∕\mathrm{ml}$ 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\mu \mathrm{m}∕\mathrm{s}$ , black curve) remains unaltered by the procedure, while the tendency toward slower sliding velocities under the influence of sNC-2 is enhanced. The $6\mu \mathrm{m}∕\mathrm{s}$ 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,¹⁷ the resulting velocity distributions have not been analyzed yet. Here, we show that while isolated actin-HMM interaction produces a preferred sliding velocity of $\sim 6\mu \mathrm{m}∕\mathrm{s}$ , 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 $(\mathrm{SD}<1.5\mu \mathrm{m}∕\mathrm{s})$ . 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 $(\mathrm{SD}<1.5\mu \mathrm{m}∕\mathrm{s})$ , both distributions separate more clearly. In the control case, the main velocity peak at $6\mu \mathrm{m}∕\mathrm{s}$ 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.

3.5.

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 $\pm 12\mu \mathrm{m}∕\mathrm{s}$ . 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.