Microtubule group movement tracking based on endpoints by optical flow and cluster matching

Abstract

Microtubules usually form groups at high densities and can glide in the same direction on a kinesin-coated glass surface. They tend to move together (snuggling) to avoid collisions and overlaps. Tracking microtubule groups is difficult due to their intrinsically complex and nonlinear behavior, as well as their sudden appearance and disappearance. We developed a microtubule motion analysis workflow incorporating a U-Net-like fully convolutional neural network (FCN) for noise filtering, a template testing method for endpoint detection, Sparse Optical Flow (SOF) for motion detection, SOF clustering for group detection, and cluster matching for group tracking. We fine-tuned the parameters using videos generated by a microtubule gliding simulation system, and then applied this workflow to real experimental videos. This workflow significantly facilitates microtubule motion analysis.

1. Background

1.1 Microtubule and Gliding Assay

Microtubules are polar, dynamic filaments that interact with motor proteins. They play an important role in cell migration and gene regulation [1]. In general, microtubules can be considered tubular structures [2]. These tubular structures can take straight, curved, or circular shapes. They can cross each other, merge into bundles, or diffuse in different directions. As a result, microtubule groups exhibit various motion patterns, such as straight, curved, or wave-like trajectories [3].

The microtubule gliding assay [4] is a biological experiment used to investigate self-organizing microtubule motion patterns powered by molecular motors attached to a glass surface [5]. In the gliding assay, as shown in Figure 1, microtubules are driven by molecular motors fixed on the glass surface when ATP is supplied.

Figure 1. Illustration of the microtubule gliding assay movement in experiments

From bottom to top are the glass surface, molecular motors, and microtubule. On gliding assay, the reaction takes ATP and releases ADP. From a full biological perspective, ADP then turns back into ATP via phosphorylation.

The motility of kinesin [6] or dynein [7] ensembles can be investigated using microtubule gliding assays, where above them are microtubules which are long, hollow cylinders made up of polymerized alpha and beta tubulin dimers [8]. These behaviors in in vitro gliding assays have provided valuable insights into important aspects of biomolecular motor function. Both motor proteins and microtubules are in liquid. Many motor proteins, such as kinesins, are fixed on the glass surface. Once the ATP is infused, motor proteins push microtubules to move. Microtubule-kinesin interactions can be observed through the objective lens and recorded as a video or picture sequence. Since the size of microtubules and kinesins is at the nanometer scale, the video tends to have significant noise [9].

Here we define the group movement as a bundle of microtubules moving in a similar direction, such as nematic, parallel, antiparallel, and ring movements of discrete microtubule filaments. Tracking microtubule groups in a gliding assay is important for understanding the mechanisms of collective motion [10] and elucidating different aspects of microtubule behavior [11]. It is useful for molecular robot design [12]. However, because of their thin, dense, and suddenly appearing or disappearing features, manually tracking is difficult.

1.2 Gliding Assay Simulation

Real-time simulations can help in understanding the movement of microtubules. We mainly used the simulation data from Gutmann et al.’s Particle-based Microtubule Gliding Simulation System [13]. The simulation utilizes a self-propelled microtubule model combined with the Lennard-Jones potential. The attraction is explicitly short-range (LJ-like with a finite cutoff) and is included only as an effective representation of linkers at close separation. No long-range (µm-scale) interaction is assumed. Microtubules are modeled as segmented dynamical systems that evolve under the influence of various forces. So the simulation is mainly about simulating these forces, determining the next positions of each segment, and rendering them on the screen. To achieve real-time performance, balancing the 3D rendering and computational workflows is important. Both CPU and GPU are used in this simulation.

A faithful atomistic model would involve billions of atoms, and even simulations with hundreds of thousands to millions of atoms are computationally demanding in packages like GROMACS [14]. Since a bottom-up simulation is not feasible here, the practical approach is a mesoscopic / experiment-mimicking simulation. We represent each microtubule as a rigid (or semi-rigid) chain constructed by connecting many particles, and we reproduce global dynamics by tuning chain rigidity and related effective parameters. With this representation, physical collisions alone produce aggregation without requiring an explicit “flocking” rule that would otherwise be needed for point-particle models.

In this simulation, specifying absolute length and time scales is not strictly necessary to understand the qualitative global dynamics; it becomes essential mainly if the goal is to compute calibrated energies, forces, or quantitative rates in physical units. The appropriate level of scaling detail depends on the study’s objective. We rendered the simulated gliding assay dynamic picture series using the method above, with the pixel resolutions $3000\times 2000$ . Here are examples in a simulation of 100 frames with 10 seconds interval, shown in Figure 2 and supplement2 ss10.avi. At first, their initial positions are random, as well as the moving direction. Then some move together to form a bundle of microtubules.

Figure 2. Microtubules gliding simulation picture series

1.3 Gliding Assay Experiment

We recorded the microtubule gliding assay video with a kinesin concentration (in feed) of 600nM and a microtubule density of $1.1\times {10}^{-2}MTs/\mu m^2$ (corresponding to about 400 MTs in a $200\mu m\times 200\mu m$ field of view). The frame rate is 5 fps (200 ms per frame) and the pixel resolution is $1860\times 1860$ attached with a $50\mu m$ scale bar. Compared with the simulation video, the main difference is the random noise. See Figure 3 and supplement3 motilitycrop.avi in detail.

Figure 3. Microtubules gliding assay experiment picture series

2. Related Work

Recently, some researchers tried to analyze the microtubule process by tracking them. Jansen, Klara I., et al. used a live-cell marker to visualize stable microtubules throughout the cell cycle [15], and these images include various types of microtubules. Mennona, Nicholas J., et al. made an image analysis tool for dynamic dense microtubule networks and visualized in several views, such as optical flow, difference images and orientation maps [16].

For microtubule tracking, current approaches generally rely on either digital image processing or the construction of geometric models. Saban, M., et al. tracked the entire length of each resolvable microtubule, $\beta$ III-tubulin isoform. Then they used probabilistic models for describing typical activities and an HMM (Hidden Markov Model) for analyzing microtubule dynamics [17]. In microtubule gliding assays, Han, Yuexing, et al. utilized B-spline curves to model microtubules and matched them in neighboring frames based on curvature similarity [18]. To quantify yeast microtubules and spindles, Ansari, Saad, et al. detected linear and curved polymers using a geometrical scanning technique [19].

Besides these traditional methods, some researchers also tried to use machine learning. Jokobs tried using an FCN to predict the traces of dynamic processes in kymographs [20]. Reme, Raphael, et al. proposed an optically enhanced Kalman filtering method for single-particle tracking (SPT)[21]. Tran, Phu N., et al. used deep learning-based optical flow (DLOF) to measure velocity fields [22]. Moreover, Yazdi, R., Khotanlou, H created a cell tracking workflow from time lapse imaging [23].

Owing to the different situations and behaviors of microtubules, there are almost no general public labeled datasets. Also annotating the trace of microtubules is tedious and time-consuming. Furthermore, some overlapped microtubules are ambiguous and not even easy for humans to distinguish. Although these methods mentioned above can perform well in some specific situations, they are not designed for a high density of microtubules or group movement tracking. Moreover the sudden appearance or disappearance [9] which is mainly due to the vertical movement on the glass surface makes it more difficult to track. The novelty of this research is that we proposed a brand new workflow dealing with plenty of microtubules in a gliding assay video. The motion detection unit is a microtubule endpoint instead of the fixed mesh in a frame. We can see not only the motion on each endpoint, but also groups of microtubules with similar motion tendency.

3. Methodology

3.1 Overview of the workflow

The workflow for tracking microtubule gliding assay video contains microtubule segmentation (denoising), endpoint detection, motion detection, group detection and group matching. In the view of individual motion, we can see the motion distribution by motion detection on endpoints. In the view of group motion, we can see microtubules with a similar moving direction by group detection and track them by group matching between frames.

As the gliding assay experiment contains many molecules and noise by Atomic Force Microscope (AFM) [24] or Fluorescent Optical Microscopy [25], microtubule segmentation (denoising) should be applied at first. After preprocessing, the image should be in high contrast gray scale image mode. Then we detect endpoints by the neighborhood template testing method. Next we use Sparse Optical Flow (SOF [26]) to detect the motion on endpoints. Finally we involve group detection and group matching. Microtubules can either move together in a direction or randomly in any direction, or they can even chase each other in a circular shape. Commonly, we define the group according to their positions and moving directions, viewing the group using a bounding box. We use cluster algorithms for group detection and mapping. One group can merge into multiple groups in the next frame, and vice versa.

3.2 Denoising Neural Network

Our idea for denoising is by deep learning, training from the simulation data with random noise. Inspired by the Fully Convolutional Neural Network (FCN) [27] and U-Net [28], we built a U-Net-like neural network with skip layers, which can converge fast with microtubule video denoising. The structure is as Figure 4, in which it diminishes the image (by the diminishing part) and then enlarges the image (by the enlarger part). Unlike the conventional denoising method, our network denoises by instance-wise biological segmentation [29]. The target is to determine whether a pixel belongs to a microtubule. As microtubules are usually small filaments [30], we use smaller pieces for input to decrease the amount of microtubules in one piece. This makes denoising more efficient. The details of the denoising neural network are in the appendix "Denoising Neural Network Structure".

Figure 4. The schematic diagram of denoising neural network

The rectangle represents the convolution layer. The height of each rectangle represents the size of layers. And the connection of the neighbor convolution layers is a convolution calculation.

The input size is $256\times 256\times 3$ , divisible by 32, for 4 times max pooling (diminishing width and height to half) and then 4 times trans-convolution (deconvolution [31], enlarge the width and height twice) in pairs. This technique makes the output size the same as the input, regardless of the input size, as long as it is divisible by 32. In order to control the balance of GPU memory load and feature extraction ability, we use symmetrical structure between diminishing part and enlarger part, and the number of convolution layers and convolution kernels becomes larger when the intermediate output size becomes smaller. Additionally the skip layers are added to merge (by concatenating) the same size intermediate output between the diminishing part (some restored convolution layer output) and the enlarger part. This helps the network remember more details among each stride scale.

We trained the network with simulated images with Gaussian random noise. As the network can recognize microtubule pixel after training, we can infer microtubule pixels on gliding assay experiment image. It is required to split the whole image into pieces matching the network input dimensions in both training and inferring process.

  

$$\{\begin{array}{c}\hfill ϵ(X)\sim 𝒩(\mu ,\sigma )\hfill \\ \hfill \hat{X}=\frac{60}{255}X+ϵ\left(X\right)\hfill \\ \hfill \mu =60,\sigma =\sqrt{0.5}\hfill \end{array}\left(1\right)$$

We use equation (1) to generate the image with noise. $X$ is the clear image with only microtubules (usually using simulation images for ground truth), while $\hat{X}$ is the generated image with noise; $ϵ\left(X\right)$ is the random noise based on a normal distribution. We want to make the noise signal intensity (brightness) similar to the microtubule signal intensity, and $\hat{X}$ is close to 120 (nearly half of the maximum brightness 255). The generated image by this equation is as Figure 5, showing the comparison of original image and noise-added image.

Figure 5. The comparison of the origin simulation piece and the generated piece with noise

3.3 Endpoint Detection

Typical feature points include corner points, end points, and cross points. Given that microtubules exhibit a slender, smooth, and filamentary morphology, we designate the endpoints as feature points for motion detection. The computation is faster because it only thins microtubules to a 1-pixel-width skeleton, and uses a template to filter the endpoint from the skeleton. Moreover, it can be faster by parallel optimization using GPU devices.

In morphology, commonly used thinning techniques are generally categorized into simple methods, such as repeated erosion and dilation, and more sophisticated algorithms, such as the ZHANGSUEN [32] and GUOHAL [33] thinning algorithms. We have tested these methods on a piece of simulation data and the results in Figure 6. This figure shows the comparison of different thinning methods. In our experiment on that piece, the time for the erosion and dilation method, the ZHANGSUEN method and the GUOHAL method is 0.9ms, 4.8ms, and 4.6ms respectively. It seems that the simple method (erosion and dilation) is much faster, but in the sophisticated method, both of ZHANGSUEN’s and GUOHAL’s computation speed are similar. The defect of thinning methods is gaps and convex points. Gaps break a microtubule into several parts, while convex points add small branch to a microtubule. As for the quality, the gaps and convex points numbers in erosion and dilation, GUOHAL, ZHANGSUEN are roughly 106, 24, 20. The simple method (erosion and dilation) performs not well due to many gaps on a single object. GUOHAL and ZHANGSUEN’s result is similar. We use the GUOHAL method to get the 1-width skeleton of microtubules.

Figure 6. Some thinning methods to get the skeleton of a piece

The next step is to detect the endpoint from the 1-width skeleton. In theory, if 1-width skeletons are extracted perfectly (no gaps or convex points on a microtubule skeleton), an endpoint must have 6 or 7 continuous empty pixels among the 8-neighbor pixels. We use equation (2) and equation (3) to demonstrate the endpoint neighborhood relation features.

  

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}\hfill p_8\hfill & \hfill p_1\hfill & \hfill p_2\hfill \\ \hfill p_7\hfill & \hfill p_0\hfill & \hfill p_3\hfill \\ \hfill p_6\hfill & \hfill p_5\hfill & \hfill p_4\hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill \end{array}\right]\left(2\right)\hfill \\ \hfill \{\begin{array}{c}\hfill \exists p_{i_0},p_{i_1},\cdots p_{i_5}=0\hfill \\ \hfill \exists p_j\ne 0\hfill \\ \hfill i_0,j\in \{1,2,\cdots 8\}\hfill \\ \hfill i_n=(i_0+n-1)mod8+1\hfill \end{array}\left(3\right)\hfill \end{array}$$

In equation (2), $p_0$ is a pixel and $p_1\sim p_8$ are the 8 neighbor pixels of $p_0$ arranged clockwise. And the equation (3) is the necessary condition if $p_0$ is determined as an endpoint. $p_{i0},...p_{i5}$ are $p_0$ ’s continuous neighbors. That means the endpoint should be at least connected to 1 neighbor and at least not connected by 6 continuous neighbors. Under this endpoints detection, there are 4 main different pattern cases (assuming rotated patterns are considered the same), as shown in Figure 7.

Figure 7. Endpoint pattern by 6 or 7 continuous empty method, according to equation (3)

The blue central grid signifies an endpoint. White grids are zero pixels. Red grids are non-zero pixels.

However, only using equation (3) is not enough, as the thinning methods are not perfect. There are still some gaps and convex points left. In Figure 8 (b), gaps and convex points by thinning method creates false endpoints. So that we develop a template pattern filter method to reduce these defects. This method detects more pixels around the point so that they can observe the gaps and convex points, allowing us to remove fake endpoints. It requires a $7\times 7$ judgement matrix and each neighborhood pixel matrix. As all of the gaps or convex points in Figure 8 (b) are below 3 pixels length, $7\times 7$ judgement matrix is suitable. And the neighbor should be rotated to a certain angle, transforming it to the basic cases in Figure 7, to adapt to the fixed judgement matrix. Then use matrix element-wise calculation after removing some special cases, and see if all the values are larger than zero. If they satisfy, the center point is a true endpoint.

We use equation (3) to roughly test every pixel and find candidates of endpoints. Then we apply equation (4) on these candidates to judge whether they are true endpoints. (x, y) is the endpoint candidate’s coordinate by equation (3). $J$ is judgement matrix containing values $\{-1,0,1\}$ to test the target value by pixel. $A$ is a $7\times 7$ area matrix of (x, y) in the 1-width thinned image I with zero padding. $J*A$ is the element-wise matrix multiplication, thus the result $R$ is a matrix. Note that $A$ should be rotated to the pattern case as Figure 7. The details of these patterns are in appendix "Endpoint detection template pattern". Figure 8 (c) shows the image with a $7\times 7$ area matrix filtered.

Figure 8. As a defect of the thinning method, there are gaps and convex points mistaken as endpoints

It shows the comparison of endpoints without or with a filter.

  

$$\begin{array}{c}\hfill \{\begin{array}{c}\hfill R=J*AJ,A\in {\mathbb{Z}}^{7\times 7}\hfill \\ \hfill \forall a_{ij}\in A,a_{ij}=I(x+i-3,y+j-3)\hfill \\ \hfill \forall r_{ij}\in R,r_{ij}\ge 0\hfill \\ \hfill \hfill \end{array}\Rightarrow (x,y)\mathrm{\text{is a true endpoint}}\hfill \end{array}\left(4\right)$$

3.4 Optical Flow

Optical Flow is a pattern of apparent motion in objects, surfaces and edges, in a visual scene caused by the relative motion between an observer and a scene. In computer vision, they can also be defined as the distribution of apparent speed of movement with brightness pattern changes between frames [34]. In other words, finding the vector field which describes how the image is changing with time.

Several strategies exist for computing optical flow, including feature-based [35], correlation-based [36], and gradient-based approaches [37]. Considering the drawbacks that the feature-based approach needs detection features (such as image edge, corner), correlation-based approach does not perform data reduction and is computationally expensive. In contrast, the gradient-based method uses spatial and temporal partial derivatives to estimate image flow at every position in the image. As microtubules are small filaments and their motion makes the surrounding pixels change significantly, we chose the gradient-based optical flow method. The details of this method are in appendix "Gradient-based Optical Flow Calculation".

In the microtubule gliding assay process, microtubules can be considered as moving on the surface of the gliding assay and with few occlusions above them. Also, the brightness has no significant fluctuation in the video. Therefore optical flow is effective to detect the movement of microtubules. There are mainly two types of optical flow in implementation, Dense Optical Flow (DOF) [38] and Sparse Optical Flow (SOF) [39]. The DOF calculates motion on all pixels, while the SOF calculates motion only on some pixels in which the user is interested. We firstly used the DOF to view the movement intuitively in all pixels, and then used the SOF on the endpoints, as a basis for defining microtubule group in a next step.

3.5 SOF Cluster

We have made motion detection on endpoints by SOF, then we will detect microtubule group with similar moving tendency by a clustering method. As there are plenty of microtubules in a gliding assay experiment, the group motion can simplify the moving information by group average motion, and also it can reduce some error caused by the individual movement detection. In this part, we will show our group definition and the cluster algorithm to find groups.

We define the microtubule group by the position and orientation similarity of each microtubule in that group. The position of similarity can be measured by the distance of the two endpoints, and the orientation similarity can be measured by the angular distance of the two endpoints’ optical flow, as shown in equation (5). We use $dis$ operator to calculate the distance between objects, and the specific equation depends on the object’s type (like function overriding in programming). For example, $dis(e{p}_i,e{p}_j)$ is the distance between endpoints, while $dis(e{p}_i,g{p}_j)$ is the distance from an endpoint to a group. $p_i,p_j$ is the endpoint position, and $o_i,o_j$ is the endpoint orientation (calculated by SOF). $\underset{\_}{p_j},\underset{\_}{o_j}$ are the average microtubule position and orientation in a group. $max\left(d_p\right)$ is the max value of endpoint displacement detected by optical between two frames. $mix(d_p,d_o)$ is the mixing function of position distance and orientation distance. As the fundamental unit of the distance and angle is different, adding together is meaningless in theory. So that we need to normalize them, by using a max range value to divide and eliminate the fundamental unit. Then we can use the result by the linear combination of distance and angle (In programming, to make the angle range and distance range simple, we still do not normalize, as all the frames in a video follow the same range criteria).

  

$$\{\begin{array}{c}\hfill dis(e{p}_i,e{p}_j)=mix(\parallel p_i-p_j{\parallel }_2,\parallel o_i-o_j\parallel )\hfill \\ \hfill dis(e{p}_i,g{p}_j)=mix(\parallel p_i-\overline{p_j}{\parallel }_2,\parallel o_i-\overline{o_j}\parallel )\hfill \\ \hfill mix(d_p,d_o)=\frac{d_p}{U_p}+3\frac{d_o}{U_o}\hfill \\ \hfill U_p=max\left(d_p\right),U_o=2\pi \hfill \end{array}\left(5\right)$$

Under this definition, an endpoint is in a group, if and only if both its position distance and orientation distance between this endpoint and group are smaller than the threshold ${\theta }_p$ for position, and ${\theta }_o$ for orientation. This group definition can be used in many situations, such as parallel movement, anti-parallel movement, cross direction movement, and ring movement.

Finding the microtubule group can also be described as finding the endpoints’ sparse optical flow cluster (SOF Cluster). Inspired by the Mean-shift method [40], we developed a cluster method for optical flow vector. It does not restrict the number of clusters, and generate the cluster by the threshold provided. Furthermore, comparing each endpoint with the existing clusters, the complexity is no more than $O\left(n^2\right)$ . We need to know the similarity of optical flow vector, and find a proper way to determine if an optical flow vector belongs to either a previous cluster or form a new cluster. We implemented algorithm (1) (the similarity on SOF vector) and algorithm (2) (SOF cluster method to find groups) for the SOF cluster.

The amount of groups (how many groups in a picture) with the variation of angular and distance threshold is shown in Table 1. This table shows the group number according to different pixel and angle parameters. Groups with a single item are excluded because they are trivial and meaningless. In general, the more groups, the lower average similarity movement in a group. The group number becomes smaller as the threshold increases. But when ${\theta }_d$ is small, adding ${\theta }_o$ even creates more groups, as some endpoints join into the group with only 1 item. We chose the middle point (300, 30) as the threshold, so that the group size and average similarity is balanced.

Table 1. Group number table (excluding the group with only 1 item) when the endpoint number is 647

The row is the threshold of orientation (degree), and the column is the threshold of displacement (pixel).

3.6 SOF Cluster Matching

It is challenging for microtubule group matching, because not only can microtubules suddenly appear or disappear [41], but many groups also join together or a group can separate into many groups, as shown in Figure 9. The points aligned in the same column represent the groups in the same frame.

Figure 9. An example of an SOF cluster group matching graph, illustrating cases of sudden appearance, disappearance, splitting, and merging

The vertex means microtubule SOF cluster, the edge means matching of neighbor frames.

Dealing with these problems, we take the relation of the current frame and the previous frame as a complete bipartite graph, in which the vertex refers to the SOF cluster, the weight of edge refers to the similarity of SOF clusters, and the SOF cluster of current frame and previous frame are in the two parts, respectively. We define that the similarity value between each SOF cluster is the similarity value between the center mass of each SOF clusters (Similarity definition is declared by the equation (5)). The center mass of a SOF cluster, also called average optical flow, is from the equation (6). In that equation, $C$ represents a microtubule group (cluster); $p_{start}$ represents the current position of an endpoint, while $p_{end}$ represents the next position of this endpoint. This definition takes advantage of the group feature in a cluster and computationally cheap. Then, the matching steps are as algorithm (3), after that we can get the SOF cluster matching graph like the example in Figure 9.

  

$$C_{center}=\frac{1}{\left|C\right|}(\Sigma p_{start},\Sigma p_{end})\left(6\right)$$

The relationships of frame, SOF Cluster, and SOF Cluster matching are illustrated as a graph in Figure 10. This figure shows the calculation sequence from frame to SOF cluster matching.

Figure 10. The structure of frames, SOF Clusters and SOF Cluster matching graph

4. Results

4.1 Denoising

We use the simulation pieces and the counterpart generated pieces by (1) as output (labeled data for learning) and input respectively, in Figure 5. The dataset includes thousands of pieces in total, with the ratio 3:1 between the training set and validation. Furthermore, we use categorical cross-entropy as the loss function and the Adam algorithm with a learning rate ${10}^{-4}$ as optimizer. It converges fast (in Figure 11) and the results are shown in Figure 12. As the accuracy is high, though the edge is kind of coarse, the predicted result is near the ground truth.

Figure 11. The training curves of the denoise neural network

Figure 12. The processed image by the denoising neural network in the simulation data with added noise

Table 2 shows the comparison of the classical denoising method (by the default NL-means [42] method) and our denoising neural network method on microtubules gliding simulation with noise generation. We use Mean Squared Error (MSE), Peak Signal-to-Noise Ratio (PSNR) [43] and Structural Similarity Index Measure (SSIM) [44] metrics to show the similarity between target image and origin image. MSE calculates the error on each pixel, while PSNR and SSIM calculate at the signal level. As our denoising method considers the problem at the segmentation level using a neural network, the result is better than classical denoise method.

Table 2. The comparison of a noised image, an NL-means image and our denoised image on a microtubule gliding simulation

MSE, the higher the worse. PSNR, the higher the better. SSIM, range $[0,1]$ , the higher the better.

image	MSE	PSNR/dB	SSIM
noised	95.88	28.31	0.012
NL-means denoise	87.45	28.71	0.013
our denoise	5.81	40.48	0.638

Figure 13 shows the denoising result on microtubule gliding assay experiment. In the view of defects (imperfect parts of denoising), the classical method creates some of the haze around microtubules, while our method leaves some single points in the empty area. Although we can use some threshold to somehow eliminate the haze, the threshold changes according to the microtubule and the environment brightness. In contrast, our method can remove the noise around microtubules, and the single points left can be removed by erosion and dilation⁴⁵ no matter how the environment brightness changes.

Figure 13. The comparison of NL-means and our denoising neural network microtubule assay experiment

4.2 Endpoint Detection

We select some frames (the first frame, middle frame, the last frame, correspond to Figure 2) in the simulation video and then test our endpoint template testing method on these frames. We count all of the detected endpoints programmatically and count failed endpoints manually. There are two types of failed endpoints, false-positive (FP) for non-endpoints detected as endpoints and false-negative (FN) for endpoints not detected. As for the performance, we tested them on both CPU and GPU.

The endpoint detection result is in Table 3. The "MT Density" column is the microtubule density at pixel level (the percent of pixel microtubules occupied). Assuming that each microtubule has 2 endpoints (ignore microtubules at the boundary), we roughly estimate the amount of microtubule ("Est. MT" column) by $(N_{detected}+N_{FN}-N_{FP})/2$ . Our method performs well if the microtubules are separate, in Figure 14 (a) as all the endpoints are correctly detected. Sometimes, when microtubules go across each other, some wrong endpoints might be detected in Figure 14 (b). But when microtubules go together to form a large bundle (in Figure 14 (c)), the endpoints inside the bundle will be hard to detect, because the microtubules seem to be connected.

Table 3. Endpoint detection performance in microtubule simulation video, MT for microtubule, EP for endpoint, FP for false-positive, FN for false-negative, Tcpu for executed time using CPU, Tgpu for executed time using GPU

Frame	MT Density	Est. MT	Detected EP	FP EP	FN EP	$T_{cpu}$ /s	$T_{gpu}$ /s
1	3.6%	540	1018	19	82	13.04	1.36
50	3.1%	441	794	5	93	13.25	1.33
103	2.4%	366	643	6	95	13.08	1.37

Figure 14. The Endpoint Detection examples. Blue marks denote correct endpoints detected (true-positive, TP)

Yellow marks denote endpoints not detected (false-negative, FN). Red marks denote wrong endpoints detected (false-positive, FP).

4.3 Individual Motion Detection

In this part, we will at first show the Dense Optical Flow (DOF) image to view the velocity of each pixel intuitively, and then Sparse Optical Flow (SOF) image only on endpoints. We use Gunnar Farneback’s algorithm [38], one of the efficient ways to calculate every pixel’s optical flow. The results are shown in Figure 15 (a). We find that the endpoint is more sensitive in optical flow, as the gradient change of them is large, but most of the center parts of a microtubule are detected wrongly as not moving (black represents zero velocity) because of the aperture problem [46]. The aperture problem is a visual perception and computer vision issue where the direction and speed of a moving object are ambiguous when viewed through a small window (aperture) because only a small part of its contour is visible, leading to local motion signals that don’t reveal the true global motion.

Figure 15. The dense and sparse optical flow in a microtubule gliding simulation

Owing to the computationally expensive and aperture problem, it is better to use SOF than DOF, in which we only need to detect the feature points’ movement. The demonstration of SOF is shown in Figure 15 (b). In this figure, color or arrow direction represents the velocity direction, and brightness or arrow length represents the velocity intensity.

4.4 Group Motion Detection

Each node in the graph (in Figure 9) is calculated by the two neighbor nodes in the lower layer. The SOF cluster matching is as Figure 16, showing the groups and how groups move by matching groups between frames. In this figure, rectangles represent the boundaries of each group, and numbers above each group are marked in $Grou{p}_{current}[Grou{p}_{previous}]$ format for showing matching. For example, we can observe that the "group with index 71" in (a) will move to "group with index 77" in (b). As for the group matching accuracy, we use the criteria $\frac{N_{correct\_matched}}{N_{allgroup}}$ to evaluate the accuracy in a frame. Considering group separation and merge, we defined the correct situation by formula 7 in most cases except two special cases mentioned later. In this formula, $C_i$ is the SOF cluster in the current frame and $C_j$ is the matched SOF clusters in the next frame. In addition, the group motion detection video is supplement4 ss10groupmatch.mp4.

Figure 16. SOF cluster matching neighbor frames example, in simulation data

The number on the top-left of a bounding box is the group index, number in the ‘[]‘ is the matching group index in the previous frame.

  

$$\{\begin{array}{c}\hfill C_i=C_{i1}\cup C_{i2}...\cup C_{in}\hfill \\ \hfill C_j=C_{j1}\cup C_{j2}...\cup C_{jn}\hfill \\ \hfill C_i,C_j\notin ⌀\hfill \\ \hfill |C_i\cap C_j|\ge 0.8\left|C_i\right|\hfill \\ \hfill |C_i\cap C_j|\ge 0.8\left|C_j\right|\hfill \end{array}\left(7\right)$$

There are two special cases of matching, one is that the group in the current frame is merged from multiple groups in the previous frame (multi_matched), the other is that the group in the current frame is from no group in the previous frame (no_matched). If one group is from multi_matched, the matching is correct when all the matched groups are correct. If one group is from no_matched , the matching is correct when this group is newly appeared (for example, moving into the observed area). We tested some frames tracking in the simulation, counted the number of no_matched, multi_matched and correct_matched groups, and the results are shown in Table 4.

Table 4. SOF cluster matching accuracy in simulation video, including first frame, middle frame, last frame

With cluster threshold θd=300, θang=30, cluster matching threshold θ=330. In the table, no_matched means it hasn't found the similar SOF cluster in the previous frame, multi_matched means it has found multi similar SOF cluster in the previous frame.

frame	endpoint	group	no_matched	multi_matched	correct	acc
0	721	104	22	12	77	74.03%
3	713	104	31	9	73	70.19%
6	668	101	23	11	81	80.20%
7	697	113	42	13	74	65.48%

As for the performance evaluation, we firstly used only single thread on CPU (python code only) and then used GPU to accelerate (python in company with cuda kernel code). The result about a processing time per frame is in Table 5 (Note that it does not include the denoising process, because the simulation does not need denoising. The denoise process requires about 2 seconds when the network is loaded) based on Intel I7-6900K and NVIDA Geforce GTX TITAN X. In this table, $N_{endpoint}$ and $N_{cluster}$ are the number of endpoints and clusters; $T_{endpoint}$ and $T_{sof}$ are the time for each step in the individual movement detection; $T_{cluster}$ and $T_{match}$ are the time for each step in the group movement detection. In CPU mode, the bottleneck is in the endpoint detection, because every point and its neighborhood are required to be calculated. Therefore, using GPU parallelism can significantly improve the performance, making it nearly 8x faster. At that time, the main time cost is about CPU and GPU memory transfer. Due to the relatively light calculation in other processes, using only the CPU is tolerable. The overall tracking process performance is about 4s/frame.

Table 5. The performance of this microtubules tracking system about both CPU and GPU on simulation and microtubules gliding assay experiment. Resolution 1860 × 1860 is the experiment, 3000 × 2000 is the simulation

5. Discussion

In microtubules gliding assay, microtubules come together in a bunch, move in a direction under some conditions. Tracking microtubules with similar motion can not only be useful for microtubule behavior analysis but also for reducing calculation and detection errors. It seems that the most ambiguous thing is the definition of a microtubule group, as there are many situations. Inspired by the K-means cluster methods, we considered the microtubule groups as clusters related to the moving direction and position. Then the definition becomes similar to the cluster method, and we can use a quasi-cluster algorithm to generate the group, whose size and number can be adjusted by the threshold.

The difference between this research and normal object tracking is that normal object tracking deals with the object’s projection taken by a camera from 3D world, but our research deals with plenty of microtubules observed by Fluorescent Optical Microscopy. The normal object tracking method takes more interest in rotation, scale, illumination, occlusion etc. However microtubules are in thin, linear, curvature structure, usually occurring in high density, also with noise. Under the consideration of these features, we use the bottom-up method, that detects endpoints as a base, finds the moving trend by SOF, and make use of the statistical information (cluster generation and matching) of the base for tracking. This behavior is independent of the specific curvature of the microtubule and concentrates on the motion information for tracking. Avoiding making the curvature model of microtubules can reduce a lot of calculations, also the motion information directly shows the changes of each frame. So that our method shows the balance of both performance and accuracy, up to 4s/frames and 75% accuracy with an old type CPU and GPU. Furthermore, the workflow can be overlapped and more processes can be optimized, then it may be able to be used for the real time analysis.

We made this workflow from the view of computer science (image processing). It can be a general solution for small filament motion detection based on endpoints. We also tested the public dataset of Microtubule Gliding Assay [47]. The dataset contains videos from gliding motility assays using kinesin-1 and 29% rhodamine tagged microtubules. Here’s the example of "Kinesin-1 gliding motility assay, whole casein passivation.avi". Because microtubules in this movie move slowly, we sample every 5 frames. As the image resolution is not so high ( $478\times 360$ ), we use the threshold ${\theta }_p=100$ and ${\theta }_o={60}^{\circ }$ . Figure 17 contains the 30th and 35th frames in that video. Figure 18 shows the individual endpoint motion by SOF while Figure 19 shows the group motion of these endpoints.

Figure 17. The original frames from "Kinesin-1 gliding motility assay, whole casein passivation.avi"

Figure 18. The individual motion detection from "Kinesin-1 gliding motility assay, whole casein passivation.avi"

Figure 19. The group motion detection from "Kinesin-1 gliding motility assay, whole casein passivation.avi"

Our method is in a new perspective of tracking microtubules, which concentrates on groups. And it is also significant for the microtubule analysis in biological behavior. The pattern based endpoint detection can roughly count the amount of microtubules. The individual motion detection can show the overall tendency of microtubules and the group motion detection unveils the distribution of microtubules and how they converge or diverge. We can infer what situations they were in or what reactions happened by the distribution of the movement. And then predict which biological process will probably happen. The speed, strength, and direction can be easily and clearly seen by the wind rose [48] illustration, as in Figure 20, which shows the velocity distribution in each direction by the radius and color.

Figure 20. The motion WindRose illustration in a frame

The radius length represents distribution ratio (unit percent), the color represents different speed range (unit pixel/frame).

6. Conclusion

In this research, we proposed a microtubule tracking workflow in a new perspective focusing on both individual and group motion based on endpoint. We showed the definition of microtubule groups and explained each part of the workflow in detail, including microtubule segmentation (denoising), endpoint detection, motion detection, group detection and group matching. We also used CUDA to accelerate the endpoint detection process, the most computationally heavy part. Finally we showed the accuracy and efficiency of our method, and made a brief discussion about why we used it in that way and the comparison between normal object tracking and microtubule tracking. In addition, the statistical motion features can be seen easily and clearly by the wind rose illustration, and these features can be helpful for inferring the situation of microtubules. This reveals a new promising way of tracking, which is significant for biological behavior analysis.

The authors express special thanks to the advisors in the co-creation environment (CCE) project for molecular robotics with VR and AI technologies, especially for Prof. Akinori Kuzuya, Dr. Hidefumi Sawai (former NICT), Prof. Shigenori Tanaka (Kobe University), Prof. Takefumi Yamashita and Prof. Taro Toyota (University of Tokyo) for their valuable advice in promoting this work.

Supplement file 1 — Appendix for this research

Appendix for showing the details of this workflow including the denoise unet structure, endpoint detection template pattern and SOF cluster algorithm.

Supplement file 2 — The microtubule simulation video

The microtubule gliding assay simulation video

Supplement file 3 — The microtubule experiment video

The microtubule gliding assay experiment video

Supplement file 4 – The microtubule group tracking video

The microtubule gliding assay simulation group matching in each frame and WindRose of the speed intensity and direction distributions.

References

• [1]  Dustin, P. Microtubules; Springer Berlin Heidelberg, 2012.
• [2]  Bicek, A. D.; Tüzel, E.; Kroll, D. M.; Odde, D. J. Analysis of Microtubule Curvature. Methods Cell Biol. 2007, 83, 237–268. doi: 10.1016/S0091-679X(07)83010-X.
• [3]  Inoue, D.; Gutmann, G.; Nitta, T.; Kabir, A. M. R.; Konagaya, A.; et al. Adaptation of Patterns of Motile Filaments Under Dynamic Boundary Conditions. ACS nano. 2019, 13 (11), 12452–12460. doi: 10.1021/acsnano.9b01450.
• [4]  Duke, T.; Holy, T. E.; Leibler, S. " Gliding Assays" for Motor Proteins: A Theoretical Analysis. Physical review letters 1995, 74 (2), 330. doi:10.1103/PhysRevLett.74.330.
• [5]  Inoue, D.; Mahmot, B.; Kabir, A. M. R.; Farhana, T. I.; Tokuraku, K.;et al. Depletion Force Induced Collective Motion of Microtubules Driven by Kinesin. Nanoscale 2015, 7 (43), 18054–18061. doi:10.1039/C5NR02213D.
• [6]  Cross, R. A. The Kinetic Mechanism of Kinesin. Trends in biochemical sciences 2004, 29 (6), 301–309. doi: 10.1016/j.tibs.2004.04.010.
• [7]  Roberts, A. J.; Kon, T.; Knight, P. J.; Sutoh, K.; Burgess, S. A. Functions and Mechanics of Dynein Motor Proteins. Nat Rev Mol Cell . 2013, 14 (11), 713–726. doi: 10.1038/nrm3667. Epub 2013 Sep 25.
• [8]  Ludueńa, R.; Shooter, E.; Wilson, L. Structure of the Tubulin Dimer. J Biol Chem.1977, 252 (20), 7006–7014.
• [9]  Mahemuti, B.; Inoue, D.; Kakugo, A.; Konagaya, A. Investigation of the Microtubule Dynamics with Probabilistic Data Association Filter. In 2016 IEEE 11th annual international conference on nano/micro engineered and molecular systems (NEMS); IEEE, 2016; pp 101–106.
• [10]  Araki, S.; Beppu, K.; Kabir, A. M. R.; Kakugo, A.; Maeda, Y. T. Controlling Collective Motion of Kinesin-Driven Microtubules via Patterning of Topographic Landscapes. Nano Lett. 2021, 21 (24), 10478–10485. doi: 10.1021/acs.nanolett.1c03952.
• [11]  Zhang, D.; Hu, Y.; Chen, Y. MTrack: Tracking Multiperson Moving Trajectories and Vital Signs with Radio Signals. IEEE Internet of Things Journal 2020, 8 (5), 3904-3914. doi:3914.10.1109/jiot.2020.3025820.
• [12]  Konagaya, A.; Gutmann, G.; Zhang, Y. Co-Creation Environment with Cloud Virtual Reality and Real-Time Artificial Intelligence Toward the Design of Molecular Robots. J Integr Bioinform.2022. doi: 10.1515/jib-2022-0017.
• [13]  Gutmann, G.; Inoue, D.; Kakugo, A.; Konagaya, A. Parallel Interaction Detection Algorithms for a Particle-Based Live Controlled Real-Time Microtubule Gliding Simulation System Accelerated by GPGPU. New Generation Computing 2017, 35 (2), 157–180. doi:10.1007/s00354-017-001.
• [14]  Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. GROMACS: Fast, Flexible, and Free. J Comput Chem. 2005, 26 (16), 1701–1718. doi: 10.1002/jcc.20291.
• [15]  Jansen, K. I.; Iwanski, M. K.; Burute, M.; Kapitein, L. C. A Live-Cell Marker to Visualize the Dynamics of Stable Microtubules Throughout the Cell Cycle. J Cell Biol.2023, 222 (5), e202106105. doi: 10.1083/jcb.202106105.
• [16]  Mennona, N. J.; Sedelnikova, A.; Echchgadda, I.; Losert, W. Filament Displacement Image Analytics Tool for Use in Investigating Dynamics of Dense Microtubule Networks. Phys Rev E. 2023, 108 (3), 034411. doi: 10.1103/PhysRevE.108.034411.
• [17]  Saban, M.; Altinok, A.; Peck, A.; Kenney, C.; Feinstein, S.; Wilson, L.; Rose, K.; Manjunath, B. Automated Tracking and Modeling of Microtubule Dynamics. In 3rd IEEE international symposium on biomedical imaging: Nano to macro, 2006.; IEEE, 2006; pp 1032–1035.
• [18]  Han, Y.; Daisuke, I.; Kakugo, A.; Konagaya, A. Tracking Single Microtubules by Using b-Spline Curves and Hausdorff Distance. In 2014 4th international conference on image processing theory, tools and applications (IPTA); IEEE, 2014; pp 1–6.
• [19]  Ansari, S.; Gergely, Z. R.; Flynn, P.; Li, G.; Moore, J. K.; et al. Quantifying Yeast Microtubules and Spindles Using the Toolkit for Automated Microtubule Tracking (TAMiT). Biomolecules 2023, 13 (6), 939. doi: 10.3390/biom13060939.
• [20]  Jakobs, M.A.H.; Dimitracopoulos, A.; Franze, K. KymoButler: A Deep Learning Software for Automated Kymograph Tracing and Analysis. bioRxiv 2019, 405183. doi:10.1101/405183.
• [21]  Reme, R.; Newson, A.; Angelini, E.; Olivo-Marin, J.-C.; Lagache, T. Particle Tracking in Biological Images with Optical-Flow Enhanced Kalman Filtering. In 2024 IEEE international symposium on biomedical imaging (ISBI); IEEE, 2024; pp 1–5.
• [22]  Tran, P. N.; Ray, S.; Lemma, L.; Li, Y.; Sweeney, R.; et al. Deep-Learning Optical Flow for Measuring Velocity Fields from Experimental Data. Soft matter 2024, 20 (36), 7246–7257. doi: 10.1039/D4SM00483C.
• [23]  Yazdi, R.; Khotanlou, H. A Survey on Automated Cell Tracking: Challenges and Solutions. Multimedia Tools and Applications 2024, 83 (34), 81511–81547. doi: 10.1007/s11042-024-18697-9.
• [24]  Binnig, G.; Quate, C. F.; Gerber, C. Atomic Force Microscope. Phys Rev Lett. 1986, 56 (9), 930. doi: 10.1103/PhysRevLett.56.930.
• [25]  Sanderson, M. J.; Smith, I.; Parker, I.; Bootman, M. D. Fluorescence Microscopy. Cold Spring Harb Protoc.2014, 2014 (10), pdb–top071795. doi: 10.1101/pdb.top071795.
• [26]  Beauchemin, S. S.; Barron, J. L. The Computation of Optical Flow. ACM computing surveys (CSUR) 1995, 27 (3), 433–466. doi:10.1145/212094.212141.
• [27]  Long, J.; Shelhamer, E.; Darrell, T. Fully Convolutional Networks for Semantic Segmentation. In Proceedings of the IEEE conference on computer vision and pattern recognition; 2015; pp 3431–3440.
• [28]  Ronneberger, O.; Fischer, P.; Brox, T. U-Net: Convolutional Networks for Biomedical Image Segmentation. In International conference on medical image computing and computer-assisted intervention; Springer, 2015; pp 234–241.
• [29]  Böhm, A.; Ücker, A.; Jäger, T.; Ronneberger, O.; Falk, T. Isoo Dl: Instance Segmentation of Overlapping Biological Objects Using Deep Learning. In 2018 IEEE 15th international symposium on biomedical imaging (ISBI 2018); IEEE, 2018; pp 1225–1229.
• [30]  Wade, R. H. On and Around Microtubules: An Overview. Mol Biotechnol.2009, 43 (2), 177–191. doi: 10.1007/s12033-009-9193-5.
• [31]  Noh, H.; Hong, S.; Han, B. Learning Deconvolution Network for Semantic Segmentation. In Proceedings of the IEEE international conference on computer vision; 2015; pp 1520–1528.
• [32]  Zhang, T.; Suen, C. Y. A Fast Parallel Algorithm for Thinning Digital Patterns. Commun ACM. 1984, 27 (3), 236–239. doi:10.1145/357994.358023.
• [33]  Guo, Z.; Hall, R. W. Parallel Thinning with Two-Subiteration Algorithms. Commun ACM. 1989, 32 (3), 359–373. doi:10.1145/62065.62074.
• [34]  Horn, B. K.; Schunck, B. G. Determining Optical Flow. Artificial intelligence 1981, 17 (1-3), 185–203. doi:10.1016/0004-3702(81)90024-2.
• [35]  Castelow, D. A.; Murray, D. W.; Scott, G. L.; Buxton, B. F. Matching Canny Edgels to Compute the Principal Components of Optic Flow. Image Vision Comput.1988, 6 (2), 129–136. doi:10.1016/0262-8856(88)90008-X.
• [36]  Bergen, J. R.; Burt, P. J.; Hingorani, R.; Peleg, S. Computing Two Motions from Three Frames. In [1990] proceedings third international conference on computer vision; IEEE, 1990; pp 27–32.
• [37]  Enkelmann, W. Obstacle Detection by Evaluation of Optical Flow Fields from Image Sequences. Image and Vision Computing 1991, 9 (3), 160–168. doi:10.1016/0262-8856(91)90010-M.
• [38]  Farnebäck, G. Two-Frame Motion Estimation Based on Polynomial Expansion. In Scandinavian conference on image analysis; Springer, 2003; pp 363–370.
• [39]  Bouguet, J.-Y.; others. Pyramidal Implementation of the Affine Lucas Kanade Feature Tracker Description of the Algorithm. Intel corporation 2001, 5 (1-10), 4.
• [40]  Cheng, Y. Mean Shift, Mode Seeking, and Clustering. IEEE transactions on pattern analysis and machine intelligence 1995, 17 (8), 790–799. doi:10.1109/34.400568.
• [41]  Kleele, T.; Marinković, P.; Williams, P. R.; Stern, S.; Weigand, E. E.; et al. An Assay to Image Neuronal Microtubule Dynamics in Mice. Nature communications 2014, 5 (1), 4827. doi: 10.1038/ncomms5827.
• [42]  Buades, A.; Coll, B.; Morel, J.-M. A Non-Local Algorithm for Image Denoising. In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05); IEEE, 2005; Vol. 2, pp 60–65.
• [43]  Tanchenko, A. Visual-PSNR Measure of Image Quality. J. Vis. Commun. Image Represent. 2014, 25 (5), 874–878. doi:10.1016/j.jvcir.2014.01.008.
• [44]  Ndajah, P.; Kikuchi, H.; Yukawa, M.; Watanabe, H.; Muramatsu, S. SSIM Image Quality Metric for Denoised Images. In Proc. 3rd WSEAS int. Conf. On visualization, imaging and simulation; 2010; Vol. 53.
• [45]  Jankowski, M. Erosion, Dilation and Related Operators. In 8th international mathematica symposium; 2006; pp 1–10.
• [46]  Hess, R.; Baker, C.; Zihl, J. The" Motion-Blind" Patient: Low-Level Spatial and Temporal Filters. J. Neurosci. 1989, 9 (5), 1628–1640. doi: 10.3389/fnint.2015.00006.
• [47]  Maloney, A.; Koch, S. Microtubule Gliding Assay, 2011. http://hdl.handle.net/1928/12559 (accessed 2024-07-21).
• [48]  Crutcher, H. L. On the Standard Vector-Deviation Wind Rose. Journal of Atmospheric Sciences 1957, 14 (1), 28–33. doi:10.1175/0095-96.