## Abstract

The detection of the precise movement of a vesicle during transport in a live cell provides key information for the intracellular delivery process. Here we report a novel numerical method for analyzing three-dimensional vesicle movement. Since the vesicle moves along a linear cytoskeleton during the active transport, our method first detects the orientation and position of the cytoskeleton as a linear section based on angle correlation and linear regression, after noise reduction. Then, the precise vesicle movement is calculated using vector analysis, in terms of rotation angle and translational displacement. Using this method, various vesicle trajectories obtained via high spatiotemporal resolution microscopy can be understood..

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

When the extracellular molecules are taken up by a live cell as forming vesicles, they are gradually delivered from the membrane into the cytoplasmic area, following the biochemical pathway. This process is called endocytosis [1]. The vesicle is a basic unit of intracellular transport, and the movement of vesicle is closely related to the understanding of the intracellular transport mechanism [2]. Therefore, for various biomedical applications related to the intracellular transport mechanism, such as drug delivery, many researches have been conducted to elucidate the precise movement of vesicle [3–5].

Earlier approaches revealed that the vesicles are actively transported via the interaction with cytoskeletons, such as microtubule and actin filament, mainly based on the purified system where the vesicle is coupled with a single transporting filament recruited by motor protein [6,7]. Those studies explained that a vesicle shows linear movement while walking along the cytoskeletal filament. Recent researches are more focused on the visualization and analysis of vesicle trajectories in a live cell with improved imaging techniques [8–10], since the vesicles in a live cell inevitably interact with multiple cytoskeletons rather than a single filament.

Since the size of vesicles is smaller than one hundred nanometers in diameter in case of common receptor-mediated endocytosis, the trajectory of vesicle is acquired based on fluorescence imaging microscopy. FIONA, Fluorescent Imaging with One-Nanometer Accuracy, theoretically enabled the localization of a single fluorophore with one nanometer precision in two dimensions [11]. Also, recent advances in imaging technique such as dual focus optics now allow us to estimate three-dimensional trajectory of a vesicle in a live cell with high spatiotemporal resolution in practical imaging condition [12–14]. Figure 1 shows an example of single vesicle trajectory in three dimensions acquired via the dual focus optics.

However, the movement trajectories of vesicle acquired via advanced imaging techniques still apparently consist of complex cluster of position data, in that the data not only include measurement error but also the vesicle actually interact with complicated structure of interwound cytoskeletons while being transported. Therefore, for detecting the detailed feature of the vesicle movement in a complex cytoskeleton network, one of the most important tasks is to establish a criterion which distinguishes the active transport of vesicle from random diffusion. For this purpose, several researches suggested mathematical methods, based on mean squared displacement (MSD) [15], angle correlation [16], and further with Hidden Markov Model (HMM) [17], in order to determine the active transport of a vesicle where it shows linear and persistent movement. Those studies have provided insight into building statistical model for recognizing the active transport, yet they only concern the vesicle movement itself, and thus the information between the vesicle and cytoskeleton is not included.

In this research, we present a novel numerical analysis method for vesicle movement in terms of the interaction between vesicle and cytoskeleton. Our method features intuitiveness and robustness in detecting active transport section from three-dimensional vesicle trajectory by estimating the orientation and position of the cytoskeleton, which is defined as a linear movement section. The process of the presented method is as shown in Fig. 2. First, noise filter is applied in order to remove the noise from the acquired position data. Then, the movement trajectory is divided into active transport section and random diffusion section, based on the linearity and persistency in the traveling direction. Next, the orientation and position of cytoskeleton are estimated by finding the linear axis during the active transport, via linear regression using principal component analysis. Lastly, the data points are projected onto the estimated location of cytoskeleton, which can be exploited to reveal the precise movement of the vesicles on the cytoskeleton. With the simulation, the performance of the proposed analysis method is proved. The proposed numerical analysis method is expected to help understand the complex trajectory data acquired from various microscopy imaging, unveiling the detailed shape of the vesicle movement in a complicated cytoskeleton network.

## 2. Numerical method for vesicle movement analysis

Vesicles show linear movement when they are transported by a single cytoskeleton filament, due to the linearity of the cytoskeleton structure [18]. Also, because the active transport of vesicle occurs when the vesicle interacts with cytoskeleton, one of the intuitive methods to distinguish the active transport of a vesicle is to extract the linear section from complex movement data points. To do this, we suggest a numerical approach to detect the linear section from the acquired position data distribution.

#### 2.1. Linear section detection

Prior to detecting the linearly distributed data points from complex vesicle movement trajectory, it is required to reduce the noise of the raw position data. For this purpose, Gaussian moving average is applied as a noise filter to each raw *x*, *y*, and *z* coordinates. Applying Gaussian filter is also called as Weierstrass transform in mathematics, which is achieved by transforming input data points via convolution with Gaussian function [19]. The convolution scheme is required in order to prevent the latency involved in the smoothing process. Therefore, the data points in three dimensions are individually filtered by Weierstrass transform function as in Eq. (1).

*i*denotes the time index,

*x*represents the

_{i}*i*-th noise-filtered coordinate while ${x}_{i}^{*}$ indicates the raw coordinate, and

*G*refers to the Gaussian kernel. The operator ⊛ denotes the convolution, and

*m*represents the size of the Gaussian kernel, which is defined as in Eq. (2).

By applying Weierstrass transform to three-dimensional data points according to Eq. (1), it is possible to obtain Gaussian weighted moving average of the raw data points set. Now consider the filtered position as a point P* _{i}*(

*x*,

_{i}*y*,

_{i}*z*) and the corresponding vector

_{i}**r**

*∈ ℝ*

_{i}^{3}. In order to define a linear section with the filtered data point set

**P**= {P

_{1}, P

_{2}, ..., P

*} and corresponding vector set*

_{N}**r**= {

**r**

_{1},

**r**

_{2}, ...,

**r**

*}, the interval where the data points show persistent linearity is required to be distinguished. To do this, we define a linear domain based on the local curvature of the trajectory, which can be investigated from the local angles. As shown in Fig. 3, the local angle*

_{N}*θ*is defined between the two vectors

_{i}**r**

_{i+m}−

**r**

*and*

_{i}**r**

_{i−m}−

**r**

*, the three consecutive data points separated by the size of Gaussian filter*

_{i}*m*, as a criterion for partial coherence in the direction of travel. The angle

*θ*can be calculated as in Eq. (3), by computing the inner product.

_{i}As shown in Fig. 3, if *θ _{i}* =

*π*, the local trajectory is completely linear, while ${\theta}_{i}=\frac{\pi}{2}$ means the perpendicularly changed local trajectory. Therefore, the local movement of a vesicle between P

_{i−m}and P

_{i+m}can be considered as relatively linear and persistent when

*θ*is closed to

_{i}*π*, compared to the case where

*θ*is far smaller than

_{i}*π*. Here we introduce a threshold angle

*θ*, in order to determine and extract the local linear domains from the entire trajectory data. As shown in Fig. 3, a series of data points forming

_{c}*θ*≥

_{i}*θ*is recognized as a persistent movement section, otherwise the data point set is treated as a random movement section.

_{c}In fact, the size parameter *m* for data separation and the threshold angle *θ _{c}* contain a trade-off between noise removal and excessive smoothing. Therefore, it is required to investigate the range of

*m*and

*θ*before we apply this noise reduction method to actual data. In order to provide practical reason for determining proper values of both parameters

_{c}*m*and

*θ*, simulation is conducted as shown in Fig. 4. In the simulation, three different linear segments are prepared in space as shown in Fig. 4(A), mimicking the crossed cytoskeletons located in a cell. Since the size of endocytic vesicle is up to 100 nm in diameter assuming common micropinocytosis [20], vesicle is expected to be detected around the cytoskeleton within the range of its radius. Therefore, the vesicle trajectory along the segments are generated randomly to be located between 0 and 100 nm from the base segment, as shown in Fig. 4(B). In order to determine the appropriate range of the parameter

_{c}*m*, the angle

*θ*which is defined as Fig. 3 is calculated with varying

_{i}*m*, from 1 to 40 data points. The line colored in orange shown in Fig. 4(C) represents the

*θ*calculated on the ground truth. As shown in Fig. 4(D), the proper value of

_{i}*m*can be determined by estimating the ratio of occupation with respect to the area defined by the ground truth. When the threshold angle

*θ*is equal to $\frac{3\pi}{4}$, the optimal value

_{c}*m*is around 16 data points for the simulation.

The range of angle *θ _{i}* to define the persistent movement of a vesicle is suggested as the angle from 0 to $\frac{\pi}{2}$ in the previous research [16], but the local minimum angles shown in Fig. 4(C) imply that $\frac{\pi}{2}$ is too small to recognize two different linear sections. Considering the ideal persistent movement forms

*π*while the returning movement forms 0 between adjacent data points, we conservatively suggest $\frac{3\pi}{4}$ as a threshold angle to determine the linear and persistent movement from complex position data distribution.

Note that the nonlinear movement includes random diffusion and the change in traveling direction. Using the above angle criteria, now we can define the linear sections as the interval where the angle *θ _{i}* is persistently larger than $\frac{3\pi}{4}$, while the random diffusion or the corner of direction change is defined when

*θ*is smaller than $\frac{3\pi}{4}$. Based on this definition and the

_{i}*m*as 16 data points, linear sections in the simulated data point set are detected as shown in Fig. 5(A). The data points colored in pink represent the sections detected as linear and persistent movement while the mint green indicates the groups of data at the corner of the direction change.

#### 2.2. Estimation on the orientation and position of the cytoskeleton

Now it is possible to estimate the orientation and position of cytoskeleton, in the section where vesicle showed linear and persistent movement. For this purpose, Principal Component Analysis (PCA) is applied to locate the axis of traveling direction, which can be interpreted as the traces of the interaction between vesicle and cytoskeleton. As a statistical technique in regression analysis, PCA produces the principal axis which minimizes the orthogonal errors from data points [21]. In order to apply PCA method, the data point set {P_{1}, P_{2}, ..., P* _{n}*}, which is distinguished as a linear section, is rearranged as a

*n*× 3 matrix

**P**

*. Then, the covariance matrix*

_{l}**C**for

**P**

*can be calculated as in Eq. (5).*

_{l}**M**represents the zero mean design matrix of

**P**

*[22]. By utilizing the singular value decomposition (SVD) [23],*

_{l}**C**is decomposed as in Eq. (6). where

**A**comprises the eigenvectors of

**C**and Λ indicates a diagonal matrix of the eigenvalues. When the eigenvector for the largest eigenvalue calculated from Eq. (6) is represented as

**ê**

_{1}= (e

*, e*

_{x}*, e*

_{y}*), this vector indicates the travel direction of the data point set*

_{z}**P**

*, which can be interpreted as the orientation of the cytoskeleton the vesicle interacted with. Since the direction vector*

_{l}**ê**

_{1}alone cannot determine the position of cytoskeleton, we fix this vector at the center position Q

*(*

_{c}*x*,

_{c}*y*,

_{c}*z*), the mean coordinates of

_{c}**P**

*. Therefore, the orientation and position of cytoskeleton can be represented by the addition of two points, (*

_{l}*x*,

_{c}*y*,

_{c}*z*) + (e

_{c}*, e*

_{x}*, e*

_{y}*). Let the end point of this vector be Q′*

_{z}*, as shown in Fig. 5(B). Therefore, the orientation of cytoskeleton can be defined as $\overrightarrow{{\text{Q}}_{c}{\text{Q}}_{c}^{\prime}}$.*

_{c}Additionally, the overall shape of the cytoskeleton can be estimated by perpendicularly projecting the respective initial and final data points in the linear section, P_{1} and P* _{n}*, onto the evaluated position vector of the cytoskeleton. Since the projection of a point onto a position vector can be computed by the inner product scheme of two vectors, we define the position vector of the estimated cytoskeleton as $\overrightarrow{{\text{Q}}_{c}{\text{Q}}_{c}^{\prime}}$, for simplicity. Then, the coordinates of Q

*, the perpendicularly projected point from P*

_{i}*(*

_{i}*x*,

_{i}*y*,

_{i}*z*) ∈

_{i}**P**

*onto $\overrightarrow{{\text{Q}}_{c}{\text{Q}}_{c}^{\prime}}$, can be computed as in Eq. (7).*

_{l}Since the point Q* _{i}* can be expressed as Q

*=*

_{i}*κ*(e

*, e*

_{x}*, e*

_{y}*) + (*

_{z}*x*,

_{c}*y*,

_{c}*z*) using the size factor

_{c}*κ*, the Eq. (7) can be numerically calculated as Eq. (8).

*κ*simply as shown in Eq. (9).

Accordingly, as shown in Fig. 5(B), starting point Q_{1} and ending point Q* _{n}*, respectively projected from P

_{1}(

*x*

_{1},

*y*

_{1},

*z*

_{1}) and P

*(*

_{n}*x*,

_{n}*y*,

_{n}*z*) to the cytoskeleton, can be numerically evaluated as in Eq. (10).

_{n}Therefore, the line $\overline{{\text{Q}}_{1}{\text{Q}}_{n}}$ is the estimated cytoskeleton in terms of the 3D orientation and position for the linear section **P*** _{l}*.

In order to evaluate the robustness of this estimation model, especially against the noise level and the diameter of vesicle, the discrepancy between the estimated cytoskeleton position and the ground truth was investigated. As shown in Fig. 5(C), for the estimated cytoskeleton $\overline{{\text{Q}}_{1}{\text{Q}}_{n}}$, the corresponding ground truth is $\overline{{\text{G}}_{1}{\text{G}}_{n}}$, and the angle between $\overline{{\text{Q}}_{1}{\text{Q}}_{n}}$ and $\overline{{\text{G}}_{1}{\text{G}}_{n}}$ is defined as *δ* which represents the estimation error. The distribution of *δ* according to the change in noise level and vesicle diameter are shown in Fig. 5(D). Noise level includes the system noise which influences the detection stability, and is defined as the standard deviation of the measurement value. From one hundred independent simulations with changing noise level, the median of *δ* was smaller than 1° up to 10 nm of noise level. Also, the distribution of *δ* increased according to the size of the vesicle, but the median value of *δ* was smaller than 1° when the diameter of vesicle is around 100 nm or smaller. Since the recent advances in imaging system such as super-resolution microscopy enabled high spatial resolution imaging and the diameter of vesicle is around 100 nm assuming common endosomes except for phagosomes [24], it is expected that the proposed model can be applied to various vesicle movement trajectory analysis acquired via high spatial resolution microscopy imaging.

Based on the estimated orientation and position of cytoskeleton, the angles between crossing cytoskeletons in an interwound structure also can be calculated. Assuming that a vesicle moves from one cytoskeleton to another as shown in Fig. 5(C), and the estimated start and end coordinates are Q_{1} and Q* _{n}* for the former and Q′

_{1}and Q′

*for the latter, the angle*

_{n}*ϕ*between two cytoskeletons can be computed by using outer and inner product of $\overrightarrow{{\text{Q}}_{1}{\text{Q}}_{n}}$ and $\overrightarrow{{\text{Q}}_{1}^{\prime}{\text{Q}}_{n}^{\prime}}$, as in Eq. (11).

Note that the angle *ϕ* is calculated by four-quadrant inverse tangent function, also known as atan2, in order to prevent any confusion in determining the rotation direction [25,26].

#### 2.3. Rotational and translational movement analysis of vesicle on the cytoskeleton

In order to analyze the movement of vesicle in terms of the interaction with cytoskeleton, the data points in the linear section are projected to corresponding estimated location of cytoskeleton, as shown in Fig. 6. The projection is conducted in the same manner described in Eq. (7). Figure 6(A) shows the projected points Q* _{i}* and Q

_{i+1}from P

*and P*

_{i}_{i+1}, respectively. The instant movement of vesicle around the estimated cytoskeleton at time

*t*can be explained based on the relation between two right-handed coordinate systems, Σ

_{i}*and Σ*

_{i}_{i+1}, which are defined as shown in Fig. 6(A). Therefore, the local movement from P

*to P*

_{i}_{i+1}can be described by the coordinate system transformation from Σ

*to Σ*

_{i}_{i+1}, with rotation matrix

**R**

*and translation vector*

_{i}**t**

*, as shown in Eq. (12).*

_{i}**R**

*and*

_{i}**t**

*are defined as in Eq. (13).*

_{i}Note that **R*** _{i}* concerns the vesicle rotation around the axis of cytoskeleton as shown in Fig. 6(B), while

**t**

*indicates the translational displacement along the cytoskeleton orientation. The angle*

_{i}*ρ*denotes the rotation angle of Σ

_{i}_{i+1}with respect to Σ

*around*

_{i}**ẑ**axis, and

*ψ*measures the angle between Σ

_{i}_{1}and Σ

*in a cumulative manner, as described in Fig. 6(B). The cumulative angle*

_{i}*ψ*is a useful term when vesicle shows continuous rotation around cytoskeleton, and it can be calculated as in Eq. (14).

_{i}Preparing these angle criteria is important in that vesicles are known to move spirally along cytoskeletons from single-molecule experiments, such as microtubule and actin filament [27,28]. Therefore, it is expected that we can detect any angular motion of the target vesicle on cytoskeleton, using these definitions for *ρ _{i}* and

*ψ*. Also, since Q

_{i}*can be interpreted as the actual position of a vesicle while interacting the cytoskeleton, instant translational velocity*

_{i}*v*along the axis of cytoskeleton between the time

_{i}*t*and

_{i}*t*

_{i+1}can be calculated as in Eq. (15).

Note that *v _{i}* is also an essential factor for identifying the type of motor protein involved in the interaction between vesicle and cytoskeleton.

## 3. Application to actual vesicle movement data

The proposed numerical analysis method described is now applied to the actual vesicle movement data. As an example case, the movement of an endocytic vesicle in a live cancer cell was tracked using three-dimensional imaging microscopy, and its movement trajectory was analyzed following the processes presented above.

#### 3.1. Three-dimensional position data acquisition

The cell line exploited in the experiment is KPL-4 human breast cancer cell [29], which was kindly provided by Dr. Kurebayashi (Kawasaki Medical school, Kurashiki, Japan). The endocytic vesicles were labeled by adding 4 nM of carboxylate-functionalized quantum dots (Thermo Fisher Scientific, Inc.) to the cells. Quantum dot is a semiconductor nanocrystal, which is widely applied in biological imaging as a fluorescent probe [30]. Since the cell engulfs surrounding molecules including quantum dots, it is possible to detect the position of vesicles by tracking the light emitted from the quantum dot [31,32]. The movement of vesicle in a live cell was imaged through fluorescent microscopy and reconstructed in three dimensions based on the principle of dual focus imaging optics [13]. The *x*, *y* coordinates of a vesicle were detected as a position of the Gaussian peak of point source intensity, and the *z* coordinates were evaluated by comparing the intensities of the point source on dual image planes. In the live-cell imaging, the initial position of the specimen was kept stable by using the axial position locking system [33]. The detection precision in the imaging system is within 5 nm in three dimensions, as shown in Fig. 7(A). Total imaging time was 10 s, and the images were taken at 100 Hz frame rate, which produced one thousand position data points. Note that fluorescent microscopy is not mandatory for applying the proposed analysis method, but the vesicle position data acquired via various imaging system can be exploited as long as the detection noise level and the size of vesicle are within the range of 10 nm and 100 nm, respectively, as explained in Fig. 5(D).

#### 3.2. Vesicle movement analysis

The result of the presented method applied to a representative trajectory of vesicle movement acquired in a live cell is shown in Fig. 7(B). The red dots represent the location of vesicle labeled by quantum dot, and the shape of the target cell is colored in gray scale, which was taken by phase-contrast microscopy. In order to compare the shape of interaction in a cytoskeleton network, we chose two vesicles in different location in a cell : one vesicle marked as **P*** _{c}* near the cell edge where cytoskeletons form a dense network structure, and the other vesicle

**P**

*found at perinuclear area where a few microtubules run due to the existence of a nucleus. The reconstructed trajectory of*

_{d}**P**

*and*

_{c}**P**

*in three dimensions are shown in Fig. 7(C) and Fig. 7(D), respectively. When our linear section finding algorithm is applied to the acquired trajectories, three individual linear sections were recognized for the vesicle*

_{d}**P**

*and four linear sections were detected for*

_{c}**P**

*. The linear sections in*

_{d}**P**

*were labeled as*

_{c}*c*

_{1},

*c*

_{2}, and

*c*

_{3}, and the linear sections in

**P**

*were marked as*

_{d}*d*

_{1},

*d*

_{2},

*d*

_{3}, and

*d*

_{4}. Also, the orientation and position of cytoskeleton for each linear section was evaluated as a red segment for both trajectories, as respectively shown in Fig. 7(C) and Fig. 7(D).

When the angles between the estimated cytoskeletons, which is notated as *ϕ* in Fig. 5, were computed via Eq. (11), the estimated cytoskeletons in **P*** _{c}* formed 158° and 19°, while the angles between the estimated cytoskeletons in

**P**

*were 156°, 159°, and 98°. Although we only investigated two example cases, very acute angle (19°) was detected in the case of*

_{d}**P**

*that was tracked near the cell edge, which implies that the vesicle can transfer between the cytoskeletons crossed with quite acute angles in a complex network structure. It is expected that identification of the cytoskeleton in this type of transfer and the precise relation between the angle*

_{c}*ϕ*and the intracellular location of the vesicle can be further studied in the future researches.

Also, for each detected linear section, the angular motion of the vesicle on the estimated cytoskeleton was calculated in terms of cumulative angle *ψ _{i}* according to Eq. (14). Figure 7(E) represents the change of cumulative angle

*ψ*and the corresponding translational distance of the vesicle in the linear section

_{i}*c*

_{1},

*c*

_{2}, and

*c*

_{3}, for the case of

**P**

*. Note that the*

_{c}*ψ*at the initial position of the vesicle on each linear section was set at 0°, and the translational distance of the vesicle on the estimated cytoskeleton was calculated by adding up the $\Vert \overrightarrow{{\text{Q}}_{i}{\text{Q}}_{i+1}}\Vert $ in Eq. (13). In the same manner, the change of cumulative angle

_{i}*ψ*and the translational distance for

_{i}*d*

_{1},

*d*

_{2},

*d*

_{3}, and

*d*

_{4}in the trajectory

**P**

*were computed as shown in Fig. 7(F).*

_{d}As a result, as shown in Fig. 7(E) and Fig. 7(F), the cumulative angle *ψ _{i}* in the linear sections seemed to be either increasing or decreasing for a certain extent but not until the end of linear section, in both cases

**P**

*and*

_{c}**P**

*alike. Since the changes in cumulative angle indicates a rotational movement around the cytoskeleton, the detected motion can be interpreted as a composition of both clock-wise and anti clock-wise rotation. This non-uniformity in the direction of rotation implies the possibility that the interaction between vesicle and cytoskeleton in a live cell might not be directly governed by the structure of the cytoskeleton, which induces a regular helical motion around the cytoskeleton as observed in the case of*

_{d}*in vitro*experiment [27,28]. In addition, while the vesicles were fluctuating in terms of angular motion around the axis of estimated cytoskeleton, they walked along the cytoskeleton showing linearly increasing translational distance as shown in Fig. 7(E) and Fig. 7(F). The mean speeds in the linear sections were distributed from 100 nm/s to 360 nm/s for the case of

**P**

*, while the linear sections in*

_{c}**P**

*showed values between 250 nm/s and 380 nm/s. Particularly, the translational speed between the adjacent linear sections of*

_{d}**P**

*changed considerably, in contrast to the small variance of the mean speeds in the case of*

_{c}**P**

*. Since the speed of vesicle on the cytoskeleton is determined by the interaction between the vesicle and cytoskeleton, which is mediated by motor proteins, the large variance in mean speed found in*

_{d}**P**

*can imply the difference in the type of involved motor proteins or cytoskeletons. Although there can exist any external influences to the vesicle movement from the cytoskeletons, such as the movement of cytoskeleton itself [34], the movement in the linear section occurs in a sub-second scale, thus the speed here can be interpreted independently from such influences.*

_{c}Therefore, using our numerical analysis method, complicated movement of vesicle in a cytoskeleton network can be investigated and understood in quantitative terms. If substantial data are collected, it is expected that the unsolved questions regarding the vesicle delivery in a cell, such as the identification of the types of cytoskeleton involved in the vesicle transfer and the relation between the transfer angle and the frequency of transfer, can be answered with a reliable evidence.

## 4. Conclusion

In this research, we proposed a novel numerical analysis method for detecting precise movement of vesicle in a complex trajectory, focusing on the interaction between the vesicle and cytoskeleton. Since the vesicles are actively transported along cytoskeleton as showing linear movement, our computational method started with locating the linear section in three-dimensional trajectory and estimated the orientation and position of cytoskeleton based on PCA and vector analysis. In the simulation, the proposed method showed high accuracy in detecting the linear sections, proving the effectiveness of the algorithm. Also, we suggested criteria to show the detailed feature of vesicle movement on cytoskeleton, by projecting the position of the vesicle in the linear section onto the estimated cytoskeleton using vector computation. Particularly, we showed that rotational and translational movement of the vesicle on the axis of the cytoskeleton can be explained by the transformation of coordinate system. When the presented method was applied to the actual vesicle trajectory acquired from a live cell, several linear sections were recognized from apparently complex trajectory data. The movement of vesicle on the estimated cytoskeleton was also analyzed in terms of the angular and translational movement. Since the presented numerical analysis method features intuitiveness and simplicity, it is expected to be practically applied to various trajectory data analyses which suffer from noise and complexity, and to open a door for understanding the precise movement of vesicle in terms of interaction with cytoskeleton.

## Funding

Grants-in-Aid for Scientific Research (A and B) (H. H. 23247022 and 16H04773) and Challenging Exploratory Research (H. H. 17K19343) from Japan MEXT.

## References and links

**1. **Bruce Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, *Molecular Biology of the Cell* (Garland Sciences, 2002).

**2. **A. Sorkin and M. von Zastrow, “Signal transduction and endocytosis: close encounters of many kinds,” Nat. Rev. Mol. Cell. Biol. **3**(8), 600–614 (2002). [CrossRef] [PubMed]

**3. **L. Y. T. Chou, K. Ming, and W. C. W. Chan, “Strategies for the intracellular delivery of nanoparticles,” Chem. Soc. Rev. **40**, 233–245 (2011). [CrossRef]

**4. **S. Matveev, X. Li, W. Everson, and E. Smart, “The role of caveolae and caveolin in vesicle-dependent and vesicle-independent trafficking,” Advanced Drug Delivery Reviews **49**, 237–250 (2001). [CrossRef] [PubMed]

**5. **J. Seo, K. Cho, S. Lee, and S. Joo, “Concentration-dependent fluorescence live-cell imaging and tracking of intracellular nanoparticles,” Nanotechnology **22**, 235101 (2011). [CrossRef] [PubMed]

**6. **H. Hess and V. Vogel, “Molecular shuttles based on motor proteins: active transport in synthetic environments,” Rev. Mol. Biotechnol. **82**, 67–85 (2001). [CrossRef]

**7. **R. Dixit, J. L. Ross, Y. E. Goldman, and E. L. F. Holzbaur, “Differential regulation of dynein and kinesin motor proteins by tau,” Science **319**, 1086–1089 (2008). [CrossRef] [PubMed]

**8. **M. Liu, Q. Li, L. Liang, J. Li, K. Wang, J. Li, M. Lv, N. Chen, H. Song, J. Lee, J. Shi, L. Wang, R. Lal, and C. Fan, “Real-time visualization of clustering and intracellular transport of gold nanoparticles by correlative imaging,” Nat. Commun. **8**:15646 (2017). [CrossRef] [PubMed]

**9. **M. Nozumi, F. Nakatsu, K. Katoh, and M. Igarashi, “Coordinated movement of vesicles and actin bundles during nerve growth revealed by superresolution microscopy,” Cell Rep. **18**, 2203–2216 (2017). [CrossRef] [PubMed]

**10. **Š. Bálint, I. V. Vilanova, Á. S. Álvarez, and M. Lakadamyali, “Correlative live-cell and superresolution microscopy reveals cargo transport dynamics at microtubule intersections,” Proc. Natl. Acad. Sci. U.S.A. **110**(9), 3375–3380 (2013). [CrossRef] [PubMed]

**11. **A. Yildiz and P. Selvin, “Fluorescence Imaging with One Nanometer Accuracy: Application to Molecular Motors,” Accounts of Chemical Research **38**, 574–582 (2005). [CrossRef] [PubMed]

**12. **E. Toprak, H. Balci, B. Blehm, and P. Selvin, “Three-Dimensional Particle Tracking via Bifocal Imaging,” Nano Letters **7**, 2043–2045 (2007). [CrossRef] [PubMed]

**13. **T. M. Watanabe, T. Sato, K. Gonda, and H. Higuchi, “Three-dimensional nanometry of vesicle transport in living cells using dual-focus imaging optics,” Biochem. Biophys. Res. Commun. **359**(1), 1–7 (2007). [CrossRef] [PubMed]

**14. **R. Velmurugan, J. Chao, S. Ram, E. S. Ward, and R. J. Ober, “Intensity-based axial localization approaches for multifocal plane microscopy,” Opt. Express **25**(4), 3394–3410 (2017). [CrossRef] [PubMed]

**15. **D. Arcizet, B. Meier, E. Sackmann, J. O. Rädler, and D. Heinrich, “Temporal analysis of active and passive transport in living cells,” Phys. Rev. Lett. **101**(24), 248103 (2008). [CrossRef] [PubMed]

**16. **A. W. Harrison, D. A. Kenwright, T. A. Waigh, P. G. Woodman, and V. J. Allan, “Modes of correlated angular motion in live cells across three distinct time scales,” Phys. Biol. **10**, 036002 (2013). [CrossRef] [PubMed]

**17. **M. Röding, M. Guo, D. A. Weitz, M. Rudemo, and A. Särkkä, “Identifying directional persistence in intracellular particle motion using Hidden Markov Models,” Math. Biosci. **248**, 140–145 (2014). [CrossRef] [PubMed]

**18. **M. Smith, H. Li, T. Shen, X. Huang, E. Yusuf, and D. Vavylonis, “Segmentation and Tracking of Cytoskeletal Filaments Using Open Active Contours,” Biophysical Journal **100**, 445a (2011). [CrossRef]

**19. **A. Jayachandran and R. Dhanasekaran, “Automatic detection of brain tumor in magnetic resonance images using multi-texton histogram and support vector machine,” Int. J. Imaging Syst. Technol. **23**(2), 1098 (2013). [CrossRef]

**20. **H. S. Kruth, N. L. Jones, W. Huang, B. Zhao, I. Ishii, J. Chang, C. A. Combs, D. Malide, and W. Zhang, “Macropinocytosis is the endocytic pathway that mediates macrophage foam cell formation with native low density lipoprotein,” J. Biol. Chem. **280**(3), 2352–2360 (2004). [CrossRef] [PubMed]

**21. **H. Abdi and L. J. Williams, “Principal component analysis,” Comput. Stat. **2**(4), 433–459 (2010). [CrossRef]

**22. **F. Wood, K. Esbensen, and P. Geladi, “Principal component analysis,” Chemometr. Intel. Lab. Syst. **2**, 37–52 (1987) [CrossRef]

**23. **A. Ilin and T. Raiko, “Practical approaches to principal component analysis in the presence of missing values,” J. Mach. Learn. Res **11**, 1957–2000 (2010).

**24. **S. Conner and S. Schmid, “Regulated portals of entry into the cell,” Nature **422**, 37–44 (2003). [CrossRef] [PubMed]

**25. **A. Ukil, V. H. Shah, and B. Deck, “Fast computation of arctangent functions for embedded applications: A comparative analysis,” Proc. IEEE Int. Symp. Ind. Electron. **2011**1206–1211 (2011).

**26. **J. Diebel, “Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors,” Matrix (2006).

**27. **S. Can, M. A. Dewitt, and A. Yildiz, “Bidirectional helical motility of cytoplasmic dynein around microtubules,” eLife **3**:e03205 (2014). [CrossRef] [PubMed]

**28. **M. Y. Ali, S. Uemura, K. Adachi, H. Itoh, K. Kinosita Jr., and S. Ishiwata, “Myosin V is a left-handed spriral motor on the right-handed actin helix,” Nat. Str. Biol. **9**(6), 464–467 (2002). [CrossRef]

**29. **J. Kurebayashi, T. Otsuki, C. Tang, M. Kurosumi, S. Yamamoto, K. Tanaka, M. Mochizuki, H. Nakamura, and H. Soono, “Isolation and characterization of a new human breast cancer cell line, KPL-4, expressing the Erb B family receptors and interleukin-6,” British Journal of Cancer **79**, 707–717 (1999). [CrossRef] [PubMed]

**30. **T. Kakizuka, K. Ikezaki, J. Kaneshiro, H. Fujita, T. M. Watanabe, and T. Ichimura, “Simultaneous nano-tracking of multiple motor proteins via spectral discrimination of quantum dots,” Biomed. Opt. Express **7**(7), 2475–2493 (2016). [CrossRef] [PubMed]

**31. **L. W. Zhang and N. A. Monteiro-Riviere, “Mechanisms of quantum dot nanoparticle cellular uptake,” Toxicol. Sci. **110**(1), 138–155 (2009). [CrossRef] [PubMed]

**32. **Y. Xiao, S. P. Forry, X. Gao, R. D. Holbrook, W. G. Telford, and A. Tona, “Dynamics and mechanisms of quantum dot nanoparticle cellular uptake,” J. Nanobiotechnol. **8**:13 (2010). [CrossRef]

**33. **S. Lee, H. Kim, and H. Higuchi,“Focus stabilization by axial position feedback in biomedical imaging microscopy,” Proc. IEEE Sensors Applications Symposium (SAS). 309–314 (2018).

**34. **I. Kulic, A. Brown, H. Kim, C. Kural, B. Blehm, P. Selvin, P. Nelson, and V. Gelfand, “The role of microtubule movement in bidirectional organelle transport,” Proceedings of the National Academy of Sciences **105**, 10011–10016 (2008). [CrossRef]