AMI  Vol.9 No.4 , October 2019
A Method for Astral Microtubule Tracking in Fluorescence Images of Cells Doped with Taxol and Nocodazole
Abstract: In this paper, we describe an algorithm that performs automatic detection and tracking of astral microtubules in fluorescence confocal images. This sub-population of microtubules only exists during and immediately before mitosis and aids in the spindle orientation by connecting it to the cell cortex. Anomalies in their dynamic behaviour play a causal role in many diseases, such as development disorders and cancer. The main novelty of the proposed algorithm lies in the fact it provides a fully automated estimation of parameters related to microtubule dynamic instability (growth velocity, track length and track lifetime), and helps in understanding the effects of intermediate drug concentrations. Its performance has been objectively assessed using publicly available synthetic data and largely employed metrics. Moreover, we present experiments addressing cell cultures doped with different concentrations of taxol and nocodazole. Such drugs are known to suppress the microtubule dynamic instability, but their effects at intermediate concentrations are not completely assessed. The algorithm has been compared with other state-of-the-art approaches, tested on consistent real datasets. The results are encouraging in terms of performance, robustness and simplicity of use, and the algorithm is now routinely employed in our Department of Molecular Biotechnology.

1. Introduction

This paper is focused on the estimation of dynamic instability of astral microtubules (MT) in in vivo fluorescence microscopy images. MTs are highly dynamic cytoskeleton polymers playing a pivotal regulatory role in several biological functions: intracellular trafficking in interphase cells, formation of the mitotic spindle, establishment and maintenance of cell morphology and motility [1] .

The structural elements of MTs are heterodimers composed of two kinds of globular polypeptides, α- and β-tubulin. Dimers polymerize into linear proto-filaments; 13 of them, arranged around a hollow core into head-to-tail arrays, make up the MT.

The intracellular pool of heterodimer subunits and the MT polymers are in a complex, dynamic equilibrium. Polymerization (or growing) occurs by a nucleation-elongation mechanism; free heterodimers are incorporated at the ends of a MT nucleus, thanks to non-covalent bonds enabled by the GTP (Guanosine-5’-triphosphate) hydrolysis. On the other hand, the polymerized structure releases heterodimers into the soluble tubulin pool (shrinking). As MTs are polarized elements, a plus-end and a minus-end can be recognized. The plus-end is characterized by a faster growth speed, hence the MT growth mainly happen at plus-ends.

Two main processes describe the MT dynamics: treadmilling, i.e. simultaneous growth at one MT end and shortening at the opposite end, and dynamic instability. This latter represents the spontaneous switching between sustained growth and rapid shortening (catastrophe). In a third possible state, the pause, the MT stops growing but does not depolymerize; the factors that regulate this state are still not fully clear. The MT dynamic behaviour is regulated by the concentration of free tubulin and chemical mediators such as calcium and magnesium ions. Moreover, microtubule-associated proteins (MAPs) play a regulatory role, as well as different tubulin isotypes (e.g. γ-tubulin [2] ), post-translational modifications, and tubulin mutations, responsible for several disorders.

Three types of MTs are recognized [1] . Kinetochore MTs contribute to the mitotic spindle assembly by linking to the chromosomes via a particular protein, the kinetochore. Astral MTs interface with the cellular cortex and are involved in several functions, including the spindle orientation. They exist during mitosis and in interphase cells about to enter mitosis. Their MT minus-end is linked to the cell centrosome (also defined Microtubule Organizing Center), a membranous structure located near the nucleus in interphase cells, while the plus-ends extend towards the cell cortex. Non-kinetochore MTs provide stability to the spindle.

The frequency of catastrophe events and the MT growth rate affect the effectiveness of the spindle assembly [1] . When the mitotic spindle is not correctly oriented, abnormal chromosome segregation can occur, and pathologies can arise, such as the human primary microcephaly [3] , a disorder of the neuro-development in which the patients are affected by a reduced head circumference with different degrees of intellectual disability.

MTs represent the target for many cancer chemotherapy drugs, the Microtubule Targeting Agents (MTA). In 1990s, paclitaxel (taxol) was a first-line drug for the treatment of many cancers, even though its application was hampered by heavy toxicity and resistance phenomena. Since then, several other MTAs, characterized by better toxicity profiles, were introduced in the clinical practice.

The general MTA mechanism is to perturb MT dynamics, so interfering with the mitotic spindle formation, arresting the cell cycle in mitosis, and possibly promoting apoptosis. They can be broadly grouped in MT stabilizing and destabilizing agents. Generally speaking, both these actions lead to a reduction of the MT dynamics. However, recent in vitro experiments have shown that, in the presence of given End-Binding (EB) proteins, low doses of MTAs can increase the MT dynamic instability instead. In any case, the MTA effects are highly concentration-dependent. In [4] , a mathematical model describing MT dynamics is proposed, and applied to estimate the catastrophe frequencies in the presence of several molecules. The model shows that MTAs and EBs are likely to interact in modifying the MT dynamic instability; however, the MT response to intermediate MTA concentrations is not completely clarified.

Subcellular components such as MTs, as well as their dynamic behaviour, can be analyzed in vivo using confocal laser scanning microscopy (CLSM). Thanks to a point illumination source and a pinhole in an optically conjugate plane in front of the detector, it is able to suppress the out-of-focus signal, achieving an optical resolution much better than in wide-field microscopy. This comes at the expenses of decreased image intensity, as a large percentage of light is blocked at the pinhole; actually, long exposure times are often required. As CLSM focuses a narrow light beam at a specific depth level, it achieves an extremely precise depth of focus. Multiple images at different depths in a sample can be captured, so enabling the reconstruction of three-dimensional structures (optical sectioning). The focal plane thickness is directly proportional to the light wavelength, and inversely proportional to the numerical aperture of the objective lens, but also depends on the optical properties of the specimen.

A major step forward in CLSM is related to the discovery of a naturally fluorescent protein in living organisms, the green fluorescent protein (GFP) [5] . Numerous other markers, with different spectral properties, have been engineered for labelling various types of proteins and cellular structures. Moreover, transgenic techniques can create organisms that produce their own fluorescent chimeric molecules, allowing biologists to detect specific genes, evaluate their kinetic parameters, and quantify the interactions among molecules [6] .

Nevertheless, CLSM is affected by limitations due to both instrumentation and samples [5] . The resolution of the system is limited to about 100 nm, even though it can be improved using an immersion layer (e.g. oil) between the lens and the sample [5] . Other phenomena that degrade images, namely low signal-to-noise ratio, variability of the biological samples, photo-bleaching, auto-fluorescence and photo-toxicity, will be discussed in Sect. 3.

MT dynamic features are usually studied in time-lapse images employing tracers built up with tubulin covalently linked to fluorophores: the fluorescently-tagged End Binding Proteins (EB-EGFP) [7] . In many experiments, Type 1 and Type 3 EGFP (EB1-EGFP and EB3-EGFP) are employed. Since the available binding sites for free tubulin heterodimers decrease exponentially along the MT, the fluorescence profile appears in the images as a comet-shaped object [7] . The cell cultures addressed in this paper are treated with EB3-tdTomato, a microtubule plus-end tracking protein (+TIPs) selective for the MTs plus-ends in the assembly phase. Such a marker allows one to visualize only the growing MT phase. Other phenomena such as shrinkages and pauses cannot be directly observed but only be inferred.

Although, in the last years, several tools have been proposed for MT tracking, due to the variability in the experimental conditions [5] the biologists still review the samples manually in many cases. Since the number of particles to be detected can be as large as several hundred, such a manual analysis is extremely time-consuming, hardly reproducible, strongly affected by inter- and intra-observer variability.

The objective of our method is to achieve a fully automated characterization of the dynamic behaviour of astral MTs, taking as input fluorescence confocal microscopy image stacks, in terms of MT growth velocity, track length and track lifetime. A major novelty aspect is that the algorithm is conceived to be an easy-to-use, practical tool to be employed in the daily activities of a biotechnology lab, yielding a reliability comparable to that of manual stack analysis, and requiring little or no intervention by the end-user. Moreover, we want to assess the effects of different concentrations of MTAs on cell cultures. In the experiments presented in this paper, we address two specific MTAs, namely taxol and nocodazole. The goal of this selection is twofold. From one hand, it allows to assess the algorithm performance, as the effects of taxol and nocodazole are well known at high concentrations. On the other hand, the results obtained at intermediate drug concentrations can help in the interpretation of the biological effects of such drugs. The effects of other MTAs can be similarly assessed.

This paper is organized as follows. In Section 2, an overview of the available approaches for detecting and tracking MTs in time-lapse images is discussed. In Section 3, the dataset used in this work is presented. Section 4 describes the developed algorithm. In Section 5, the obtained results are presented, and in Section 6 conclusions are drawn.

2. MTs Detection and Tracking: State of the Art

The detection and tracking of different kinds of particles (including MT) in time-lapse fluorescence image sequences have been addressed in many studies. Nevertheless, at present, there is no standard protocol to follow [8] [9] . The main reasons can be found in the extreme variability of the biological processes and the equipment used to acquire the image sequences. Moreover, even though many proposals aim at achieving flexibility as for the kind of particles to be tracked (e.g., MTs, vesicles, viruses), the extremely different nature of both the particles themselves and the motion they are expected to exhibit, makes this task even more challenging.

Most proposed algorithms are divided into four steps [6] :

1) Image data pre-processing, mainly devoted to noise reduction.

2) Particle detection, which consists in recognizing and sealing off the objects of interest from the background on a frame-by-frame basis.

3) Particle linkage for time-tracking of the previously identified objects.

4) Post-processing of the results, in order to provide quantitative information about the biological phenomenon at hand.

Regardless of the adopted solution, the performance of any algorithm in terms of accuracy, robustness and precision gets dramatically impaired when the signal to noise ratio (SNR) drops at very low levels [8] [9] . Moreover, the dominant noise that corrupts images is not additive. Other problems to tackle are the low contrast of the images, the auto-fluorescent background [10] , and the fact that objects might exit the focal plane during the experiments.

Some methods exploit a Bayesian approach by designing a filter aiming to predict the particle positions from a series of measurements. The filter design embeds both the dynamic (representing the spatial-temporal particle behaviour) and the actual measurement model. The Kalman filter is addressed [10] , which is the optimal estimator of the state of a linear system when the noise and the error affecting the models are zero-mean, normally distributed, statistically independent random variables. However, it achieves good performance even if these conditions are not exactly fulfilled. The piecewise-stationary motion model smoother (PMMS) algorithm [11] aims at tracking several kinds of molecules subject to rapid motion changes in high-density scenarios. A stochastic smoothing stage detects the particles of interest on a frame-by-frame basis using an iterative approach and assuming a Gaussian intensity model. As for tracking and trajectory reconstruction, an update of the publicly available u-track software [11] is addressed, including further motion models besides the linear and Brownian one. PMMS is suitable for the detection of objects characterized by heterogeneous or jerky motion, in images acquired with a reduced frame rate. However, the MT dynamic behaviour is well described as a linear motion, and embedding such information in a tracking algorithm seems more efficient.

In [12] , an automatic tracker employing a Bayesian probabilistic framework is proposed and evaluated using simulated sequences, for which ground truth was available. However, its usefulness in practical situations cannot be assessed.

In [13] , a non-parametric regression method for denoising fluorescence microscopy image sequences corrupted by Poisson-Gaussian noise is proposed. A global energy functional is minimized, involving spatial-temporal image characteristics. The algorithm performance, evaluated on both synthetic and real image sequences, is heavily dependent on a number of design parameters.

Other authors choose to locate particles through enhancing techniques, and reconstruct trajectories by means of a nearest-neighbour criterion [10] [14] . Most approaches implement a search strategy exploiting peak intensities [15] .

Balzarini et al. [14] use thresholding to detect particle positions, using little a priori knowledge of the motion regime; this makes the procedure less expensive from a computational standpoint. As the approach is affected by a large false positive rate, a position refinement strategy is implemented. Once the particle positions are estimated, they are linked to build up the trajectories with a nearest-neighbour criterion. Let us assume that p is a particle belonging to the i-th frame of a stack, and q a particle belonging to the successive j-th frame of the stack, with j = i + 1 . Let ( x p , y p ) and ( x q , y q ) be the spatial coordinates of particles p and q. The linking of particles through two consecutive frames is based on the minimization of a cost function that takes into account both spatial displacement and intensity information:

Φ = ( x p x q ) 2 + ( y p y q ) 2 + [ m 0 ( p ) m 0 ( q ) ] 2 + [ m 2 ( p ) m 2 ( q ) ] 2 (1)

where m 0 ( ) and m 2 ( ) represent the moments of zero-th and second order respectively of intensity of pixel p and q (see [14] for further details on the employed cost function and its rationale). This method provides good results and is quite efficient from the computational point of view. This is the reason why it has been taken as a starting point in several works, such as [16] and also in this paper. However, its performance is suboptimal when the images exhibit low quality and high particle density.

Mahemuti et al. [17] investigate the MT dynamics using morphological information for detection and a probabilistic data association (PDA) filter for tracking. Moreover, the authors perform an object decomposition in order to detect single particles making up compound structures. Once the MT positions are estimated, the tracks are identified via a probabilistic approach. Elements into consecutive frames are considered as belonging to the same MT if the measured and estimated positions are similar in direction and movement. Again, the PDA algorithm is based on the Kalman filter for data association and track updating. The technique exhibits good accuracy in environments with low density of particles, indeed the performance impairs in high-density video.

In [18] , a multiple-stage algorithm is described, whose main feature is the use of morphological transformations in the detection step in order to limit the effects of noise on the subsequent segmentation. This technique has been compared with a preliminary version of our algorithm, using the same set of real data. The two algorithms yielded comparable results as for mean velocity and track length, but [18] is subject to a larger number of false positives.

Recently, and specifically for the MT scenario, some authors have addressed the issue of very low SNR. In [19] an approach based on a Gaussian Process Regression (GPR) is addressed to perform an estimation of MT dynamic parameters (namely, speed, track length, lifetime). The GPR model heavily relies on prior knowledge on the MT motion model, so trading precision with generality. The method has been tested on both synthetic and real data.

In [20] , a robust MT tracking method is proposed, whose main feature is that the particle coordinates over the frames are treated jointly and not as independent items. This reduces the likelihood that tracks showing an unexpected jagged pattern are selected. The task is performed implementing an adaptive hierarchical energy-based trajectory smoothing approach.

In [21] , an improved robust MT detection method is presented, based on few assumptions at the object-level. As MTs appear as filaments in microscopy fluorescence images, the authors model each particle with three connected points standing for the two spindle pole bodies and the plus-end. Moreover, they focus on the effects of photo-bleaching, and provide a model incorporated into a particle filter employed to track spots.

As a general comment, we point out that, even though several methods have been proposed, they have often been validated on different datasets, and/or using different metrics. As a consequence, it is difficult to objectively compare their performance. This is the rationale behind the International Competition described in [9] , where a common (publicly available) dataset was provided to the participants. Attendees proposed their own algorithm to detect and track different kinds of particles (including MTs) in a number of simulated scenarios. Performance was evaluated using commonly defined metrics. The competition confirmed that a single best method for multi-particle tracking does not exist, as each algorithm has crucial parameters to be tuned specifically on the available dataset [9] [16] [22] . It is worth noticing that the complexity of real experimental data is so high that simulated images cannot be assumed to be fully representative; hence, the algorithm ranking obtained with such data cannot be expected to be exactly reproducible in a real experimental set up. Nevertheless, the data employed in the International Competition still remain the best ground truth to compare with, in order to assess any algorithm performance.

Finally, we briefly discuss machine learning-based approaches. Machine learning is a powerful tool that requires little or no a priori knowledge on the particles to analyze. Nevertheless, an efficient method for MT detection and tracking has not been validated yet. The main reason can be found in the fact that getting an appropriate (real-data) training set is difficult, because of the diverse nature of the objects of interest [22] . Cellular imaging is affected by many aspects that are difficult to finely control (e.g. temperature), so that each experiment exhibits features that are difficult to generalize. Moreover, a ground truth to use for the training stage is seldom available in case of real data. On the other hand, in order to build a simulated data set representative of real experiments, a comprehensive study of the physical and motion model of the particles at hand would be needed, using detailed information about the experiment, with high cost and little generalization potential (see also [23] [24] ). This explains why, despite the huge potential of machine learning, at present approaches that exploit available a priori information are still adopted in most cases. Nevertheless, some authors have recently proposed deep learning-based approaches. In [22] a convolutional neural network is used to detect sub-micro-scale particles, in order to optimize tracking procedures. The method is based on predefined tuned parameters and was tested on both synthetic and real data. Unfortunately, the approach is not suitable for the detection of filaments, such as MTs. Yao et al. [25] aimed at improving the accuracy of the tracking phase via a proper tuning of the initial algorithm parameters. Their approach is based on a recurrent neural network that learns and models the object behaviour, given a training set. However, this network has been validated on synthetic data only, and the performance in real scenarios cannot be assessed.

As a conclusion, despite a large number of proposals, at present, there is no detection and tracking algorithm whose performance fits every scenario, in terms of particle class and density, noise levels and types. From these considerations stems the novel approach of this present proposal, i.e. to abandon any claim of generality and to focus on a very sensible method, making use of all the available a priori pieces of information on MT dynamics, and set up on real experimental data acquired with the instrumentation available in loco.

3. Dataset Description

In our experiments, a dataset of 40 time-lapse sequences has been produced. The images have been acquired with a Leica TCS SP5-AOBS 5-channel confocal system, equipped with a 561 nm DPSS laser.

A HeLa-K (HeLa Kyoto) cell line expressing EB3-td Tomato, was chosen to carry out the experiments. This cell line is largely employed in the scientific research, and is the first human cancer cell line immortalized in tissue culture. It is named after Henrietta Lacks, a woman to which the scientific community owes a lot, as cells were extracted with a biopsy of the adenocarcinoma of the cervix she was affected by [26] . The cell culture has been maintained in DMEM-GlutaMAX (Invitrogen) medium supplemented with 10% fetal bovine serum, 100 U∙ml1 penicillin, 100 μg∙ml−1 streptomycin, 200 μg∙ml1 geneticin (Sigma) and 0.5 μg∙ml1 puromycin.

Our experiments were included in a larger study involving astral MTs, and aiming at clarifying their role in the mitotic spindle orientation. Hence, we have addressed interphase cells, treated with different concentrations of nocodazole and taxol: 0 nM (control), 0.1 nM (taxol only), 1 nM (nocodazole only), 10 nM and 100 nM. Nocodazole is a MT destabilizer, whereas taxol acts as a MT stabilizer. Despite their actions are opposite, at high concentrations the neat effect of both drugs is an inhibition of MT dynamic instability. Hence, in these scenarios, we expect to detect a drastic suppression of MT dynamics, and this can be exploited to validate the effectiveness of our tool. On the other hand, these experiments may help to get some insight on the drug effects at intermediate concentrations, and replicated using different MTAs.

After one-hour incubation, videos of astral MTs were acquired using the already mentioned Leica TCS SP5-AOBS confocal system. During the acquisition, cells were stored in the microscope incubator at 37˚C with CO2 5%. For each dosage, in both cases, five stacks have been acquired and saved in TIFF format. The main characteristics of the image stacks are summarized in Table 1.

Noise Characterization

Fluorescence confocal microscopy images are affected by numerous noise sources. First of all, photon shot noise, caused by the random emission of photons [5] , becomes relevant when the number of photons is so small that the uncertainty related to the Poisson distribution cannot longer be neglected [27] . In the case of fluorescence images, the source intensity has to be kept very low for several reasons. First of all, excessive light intensity can affect the living cell behaviour (if not the life itself). Then, photo-bleaching must be avoided, i.e. the fact that markers lose their capability to fluoresce over time to an extent related to light intensity and exposure time [5] . Finally, the achieved spatial resolution is related to the pinhole detector as follows [6] :

ϵ = σ N

where σ is the standard deviation of the instrument point spread function and N is the average number of photons detected in the exposure time [27] . Hence, it could be improved by choosing a small pinhole diameter detector, but this further limits the signal intensity and impairs the shot noise. Unfortunately, photon shot noise can be limited only augmenting the light intensity, hence it is unavoidable in practical applications on living cells.

Speckle noise is a multiplicative noise process that degrades images making them look grainy. It becomes relevant when coherent imaging systems are employed, such as laser in confocal microscopy, and is caused by random interferences between the coherent returns. The effect on grayscale images is an increase of mean intensity in a local area [27] . Other noise sources are autofluorescence, i.e. the property of some molecules to naturally fluoresce at wavelengths in the range of visible spectrum, overlapping with the fluorophore; background noise, caused by the ambient radiations; dark current, due to the thermal agitation of particles at high temperature inside the detector, which leads to spontaneous emissions; quantization noise of the digital output; scattering of

Table 1. Main characteristics of employed image stacks.

light, which occurs when the object dimensions are comparable with wavelength size [5] .

It is generally agreed that, if the SNR (defined as in [10] ) drops below 4 dB, the performance of virtually any algorithm is drastically impaired [9] [10] . Table 2 reports SNR values evaluated for the image stacks considered in this paper. It can been noticed that, even though the datasets exhibit some variability, most stacks are affected by noise levels below or very close to this critical threshold. Hence, it is crucial to implement an effective denoising, focusing on those noise sources (e.g. speckle) that can be effectively faced.

4. AMicro: The Proposed Algorithm

In this section, we describe our algorithm for MT detection and tracking, which will be labelled AMicro in the following. The main objective is to achieve reliable estimates of a few parameters of interest, namely average and standard deviation of:

· MT growth velocity.

· MT track length.

· MT track lifetime.

These are the same parameters that the expert biologists evaluate manually. In

Table 2. SNR values (dB) of the stacks belonging to the addressed dataset.

the control stacks analyzed to monitor astral MTs, they are able to effectively identify some dozens of tracks, used to work out the average metrics of interest. As a consequence, we have set up our algorithm so as to achieve a number of tracked MTs comparable with that of manual analysis. We trade a possibly higher False Negative Rate (FNR) with a lower False Positive Rate (FPR). In fact, as we are interested to measure average parameters, we want the selected tracks to be very reliable, even at the expenses of disregarding a number of true tracks.

With respect to state-of-the art methods, we have devoted particular attention to the following aspects.

· Robustness. The experimental conditions yield very noisy images. Hence, we focus on effective denoising, taking into account the statistical properties of the main noise sources, and focusing on those that can be more effectively limited.

· Ease of use. As it is designed to be routinely used by the biologists, we have privileged a solution with few or no parameter to be manually set.

· Computational efficiency. The method is expected to manage huge amounts of data every day. Hence, the computational complexity should not be excessive.

These features have been traded off with generality. In fact, our algorithm is only suitable for astral MT detection and is not efficient in tracking other types of particles. We have exploited all the possible a priori information of MT dynamic behaviour. Also, in order to enable a simple and efficient use of the tool, we have selected parameters to match typical experimental conditions, i.e.: detection of MTs in interphase cells, and medium MT density (about 100 MTs per field).

The main feature that differentiates between our algorithm and the previously proposed ones, is the fact that it is tuned on a specific application, conceived for easy everyday use by biologists, aiming at achieving performance comparable with those obtained by manual analysis.

AMicro has been developed in Matlab 2017a, and, following a typical approach, is divided into three main steps: enhancement, detection and tracking.

4.1. Enhancement

Our experimental data are affected by high levels of Poisson and speckle noise. As discussed, Poisson noise can only be limited by augmenting the light intensity, which is not feasible in in vivo experiments. In order to reduce the speckle noise, we have devised a simple heuristic procedure called the LOG-Wiener Transform. First, we apply a logarithmic operator to the image in order to map multiplicative noise into an additive one. Then, we process frames with a Wiener filter with a neighbourhood 3 × 3 wide, based on the working assumptions that noise and signal are not correlated, and the noise process is additive in the transformed domain. Finally, the filtered image is subject to inverse logarithmic transform. This process is very effective to limit speckle noise and has proven to provide smooth background, making the subsequent particle identification easier. However, we stress the fact that the selection of this procedure is driven by heuristic considerations, and the assumption of additive Gaussian noise in the transformed domain is not theoretically guaranteed.

4.2. Detection

This step is devoted to the detection of comets in single images of the stack.

Calibration: In order to limit the FPR, and to automatically match proper parameters with the image stack at hand, the algorithm encompasses a calibration phase. We know that astral MTs stem from the centrosome towards the cell cortex. Hence, it is relatively easy to select a portion of a sample frame not containing MTs, representative of the background, and a portion centred around the centrosome, hence representative of the signal-containing area. In the calibration phase, for each stack, the user is asked to select two regions of a sample image, respectively including and not including the centrosome. Then, the algorithm estimates the sample distribution of pixel intensities in this area, and, in particular:

· the mean value I ¯ b and the standard deviation σ b of the background;

· the mean value I ¯ o and the standard deviation of the intensity of the objects of interest (i.e., the astral MT plus-ends).

It is worth noticing that the same information is used to work out the SNR [9] :

SNR = 10 log I ¯ o I ¯ b σ b (2)

Once these parameters are estimated, we set a threshold T h on the amplitude of an object to be considered a MT plus-end. It must exceed the average amplitude level by a factor that depends on the object standard deviation. From an experimental evaluation on several images (not reported for brevity), the sample distribution of comet plus-end amplitudes exhibits a rather small standard deviation. Hence, we have set T h = I ¯ o . Variations of the threshold with respect to this value provide different trade-off between FPR and FNR.

We assume that the parameters estimated on a sample image of the stack (typically the first one) hold valid for the whole stack at hand. This is not exactly true, due to photo-bleaching. However, this choice is dictated by simplicity issues, and has proven to be rather robust.

Comet Detection: Once the threshold T h is defined, the actual detection of comets is based on a local maxima search over the frame. The search is carried out locally, employing a squared scrolling window whose side is about 400 nm (after [16] ) applied to the enhanced image. Within the k-th window in the t-th frame, a local maximum M k , t at spatial coordinates ( x , y ) is considered as a comet plus-end if and only if its intensity exceeds T h :

M k , t ( x , y ) > T h (3)

The comet positions are then refined as in [16] , by centering the squared window on the local maxima previously detected, and recalculating the peak intensity exploring the newly selected neighbourhood; this limits the problem of recognizing as split two objects actually belonging to the same MT.

A visual representation of the detection results is shown in Figure 1.

4.3. Tracking

Tracking is based on the assumption that MTs exhibit a uniform linear motion during their growth phase (i.e., the only one directly detectable). First, the coordinates of the particles identified in the previous step are linked in order to build up partial trajectories. To this purpose, plus-end positions are connected frame-by-frame minimizing a simple cost functional, where only the contribution due to particle displacement has been kept (as the intensity contributions have revealed to be little significant in our experiments due to the narrow distribution of the comet intensity):

Φ = ( x p x q ) 2 + ( y p y q ) 2 (4)

The maximum displacement allowed between two particles in order to link them is determined as follows. Let us assume that F is the frame rate (in frames/s) and V max is the maximum expected velocity of growing plus-ends (in nm/s). Hence, the maximum displacement (in nm) that can be expected from a growing plus-end between two subsequent frames is

S max = V max F (5)

In our experiments, we have F = 2 frames / s . From studies in literature and thanks to the expertise of the expert biologists, we set V max = 55 μ m / min (this

Figure 1. Cumulative MTs detected at increasing time instants, corresponding to different frames of the same stack.

value has also been validated a posteriori—see Sect. 5). Hence, we have worked out an upper bound S max 450 nm / frame , corresponding to a window 7 × 7 pixel-wide in the present experiments.

This step yields partial tracks, because a track may be temporarily lost due to the MT entering a pause state. Hence, as a common practice, once the partial tracks are available, the algorithm provides for their linking. Two partial tracks T n , T m are connected if and only if:

d ( T n , T m ) S max (6a)

t ( T n , T m ) T max (6b)

where d ( , ) is an operator that measures the spatial displacement between the last pixel of T n and the first pixel of T m , and t ( , ) is an operator measuring the time displacement (in seconds) between the two partial tracks. According to Equation (6a), the linking occurs if and only if the end of the first track and the beginning of the second one have a maximum displacement of S max as evaluated in Equation (5). The rationale behind this choice is that, in case a MT has entered a pause and then it recovers, it is reasonable to search for its plus-end within a window dictated by the maximum expected displacement between adjacent frames. This choice does not take into account the case of the track being fragmented because its plus-end has temporarily got out of the focus plane. However, the expert biologists deem the pause far more likely than this latter event. As for the parameter T max in Equation (6b), it has been set to 2.5 s, as this is considered representative of typical pause events in the experiments at hand.

As a refinement, tracks are fitted with a second-degree polynomial. Finally, in accordance with biologists, we have decided to discard tracks shorter than 2.5 s. This allows one to remove non-reliable tracks exhibiting Brownian motion, which does not meet the assumed linear motion model.

5. Results

In this section, we present the results achieved by AMicro. First of all, we provide an assessment of its performance, referring to synthetic, publicly available data used for the International Competition described in [8] . Then, we test our method on real data, and evaluate the mean MT velocity v m , track length λ m , track lifetime τ m , and their standard deviations σ v , σ λ , σ τ , for both taxol- and nocodazole-doped cell cultures. The achieved results are compared with those yielded by similar algorithms. In case the related software is released, we have run it on the same data stacks. Otherwise, we refer to published results, worked out on data comparable to that addressed in our experiments. Finally, AMicro has been compared with the results achieved by the manual processing performed by the biologists.

5.1. Algorithm Assessment with Standard Synthetic Data

In order to objectively assess the AMicro performance, we have run the algorithm on the same synthetic data and employing the same metrics as in the International Competition [8] . As the tool is not devised for tracking miscellaneous particles (e.g. also viruses, vesicles), its validation has been carried out on data related to the microtubule scenario. We have considered four increasing SNR levels (1, 2, 4, 7 dB) and the mid-density case. The synthetic data for this scenario foresee about 500 tracks per video, hence can be considered as rather demanding for AMicro, which has been designed to manage about a hundred tracks per video.

The performance has been expressed in terms of the average α and β measures and Jaccard similarity coefficient JSC. These metrics evaluate the closeness between the selected tracks and the ground truth (i.e., the tracks actually simulated); the reader is referred to [8] for more details. The obtained results are reported in Figure 2; the average values are reported to enable comparison with different algorithms. The performance of AMicro turns out to be comparable with other state-of-art algorithms [8] .

We stress that our tool, by construction, is heavily dependent on some assumptions specific of the experimental conditions at hand:

· Astral MT detection in interphase cells, hence linear motion of MT stemming from the centrosome and directed towards the cell cortex.

· Moderate MT concentration.

· Presence of speckle noise, besides Poisson and additive Gaussian one.

The simulated data are not fully representative of these assumptions, as both the motion model and the multiplicative noise are sub-optimally represented in these data.

It is worth also discussing our technique performance related to false-positive and false-negative rates. The performance of virtually any algorithm in terms of α , β and JSC are known to be more sensitive to FNR than FPR [15] . On the other hand, by construction, we have privileged low FPR, even at the expenses of a higher FNR, in order to measure parameters related to very reliable tracks. When tested on the simulated data, AMicro achieves an average FPR of about 15% and FNR in excess of 22%. These values, considered jointly with the α , β and JSC metrics reported in Figure 2, allow us to conclude that the performance of our algorithm is in line with other state-of-the-art methods (refer to [15] for more details).

5.2. Results: Nocodazole-Doped Cells

In a first set of experiments, AMicro has been tested on cell cultures doped with nocodazole 0, 1, 10, 100 nM. In this preliminary validation phase, a limited number of drug concentrations have been considered, although spamming a large range, due to limited resources. The average and standard deviations of growth velocity (μm/min), track length (μm) and duration (s) of the detected astral MT tracks are reported in Table 3. Zero values denote that no track matching the linear motion model has been identified in the stack at hand, due to the

Figure 2. Performance of AMicro in terms of average α , β and Jaccard similarity coefficient.

Table 3. Nocodazole-doped cells: results.

suppression of MT dynamics. Actually, as the drug concentration increases, fewer and fewer tracks are detected: from an average of about 100 in controls to a dozen at 100 nM. At concentration as high as 100 nM, in most movies all MTs are disassembled.

The average growth velocity increases significantly at 10 nM, and then decreases again, reaching its minimum at 100 nM concentration. This behaviour represents an interesting outcome of our work. A finer analysis of the velocity trend in cell cultures doped with more nocodazole concentrations is left to future developments. For the sake of clarity, the growth velocity boxplots are displayed in Figure 3.

Due to the decision to discard particles not matching the uniform linear motion model, most tracks are longer than 1 μm. Mean lifetime exhibits a trend similar to length, namely it decreases at 1 nM, although a small increase occurs at 100 nM. However, at that concentration only one stack is evaluated, hence this result is not statistically sound.

Figure 4 reports the sample distributions of velocity, length and lifetime in

Figure 3. Boxplots of growth velocity at different nocodazole concentrations.

Figure 4. Nocodazole-doped cells: sample distributions of velocity, speed and lifetime. Median and mean values numerically evaluated on the sample distribution.

control stacks (similar results can be obtained for other drug concentrations). The sample distribution of velocity is approximatively normal, whereas both length and lifetime exhibit an exponential decay. The mean and median values exhibit little difference, meaning that outliers have not a significant impact on the algorithm performance. It can be noticed that the velocity values are well below V max = 55 μ m / min addressed in Equation (5).

5.3. Results: Taxol-Doped Cells

Taxol is a MTs stabilizer, so it increases the polymer mass and suppresses MT dynamic instability. In Table 4, the average and standard deviation of growth velocity (μm/min), track length (μm) and duration (s) of the detected astral MT tracks are reported.

The number of detected tracks ranges from about 100 (control stacks) to 19 at the highest drug concentration addressed. Comparing the nocodazole and taxol results in control stacks (where actually no drug is employed), we can notice that

Table 4. Taxol-doped cells: results.

in the second case molecules exhibit higher mean velocity. This cannot be explained in terms of drug effects, since cells have not been doped in either case. This points into evidence the extreme variability and complexity of the problem, since cell functions are altered not only by drugs, but also by environmental factors (e.g. temperature) very difficult to finely control.

Velocity mean values exhibit a uniformly decreasing trend; this is coherent with the theoretical knowledge of the taxol effect on MT dynamic behaviour. The growth velocity boxplots are reported in Figure 5.

Figure 6 shows the sample distributions of speed, length and lifetime for a taxol concentration of 0.1 nM. Considerations similar to those related to Figure 4 still hold true in this case.

5.4. Statistical Data Analysis

In order to assess the statistical reliability of the obtained results, the standard error of the mean (SEM) has been worked out for velocity, length and lifetime. SEM is an indicator of the value variability among different experiments, and it is defined as:

SEM = σ M

where σ is the standard deviation of the distribution of the parameter at hand, and M is the sample size. In this work, for each drug concentration, the sample standard deviation has been employed, whereas M is the cumulative number of tracks detected in each stack. Table 5 reports the SEM for the three parameters taken into account for both taxol and nocodazole-doped cells. The SEM values are reasonably low, hence we can conclude that the obtained values at each concentration, are sound from the statistical point of view.

Finally, even though the three features of interest, namely velocity, length and

Figure 5. Boxplots of growth velocity at different taxol concentrations.

Figure 6. Taxol-doped cells: sample distributions of velocity, speed and lifetime.

Table 5. Standard error of the mean (SEM) for MT velocity, speed and lifetime.

lifetime, have been estimated independently, their mean values are clearly correlated. In the following, we assume that two variables (velocity and lifetime) are independently estimated, and we work out the third one (length) from the mean values of the others:

λ c = v m τ m k

where v represents the mean velocity (μm/min), τ the mean lifetime (s), and k is a conversion factor. In Table 6, such computed length values λ c are compared with those estimated by the algorithm, λ m . It is clear that the estimated length values do not significantly differ from those computed from velocity and lifetime mean values.

5.5. Performance Comparisons on True Test Data

As previously discussed, the algorithm assessment based on simulated data, even though significant, should be completed with comparisons of competing algorithms tested on the same real data. Obviously, in this case, a ground truth is not available, hence metrics such as α , β or Jaccard similarity coefficient cannot be worked out. Nevertheless, in literature other algorithms are described, that estimate average MT parameters in cell cultures comparable to those addressed in this paper, and/or whose software implementation is publicly available. Hence, we have compared the average velocity and length and the respective standard deviations yielded by AMicro with:

· Algorithm 2 (plusTipTracker) described in [28] . The software, publicly available, has been run on a subset of the same stacks of the present work.

· Algorithm 3 described in [16] . As the software is not available, we replicate the results reported by the authors, obtained on a HeLa Kyoto cell line stably expressing EB3-EGFP, and doped with 0, 80 nM nocodazole and 0, 20, 100

Table 6. Estimated vs. computed MT track length.

nM taxol concentrations.

· Algorithm 4 described in [18] . This algorithm has been directly tested on a subset of the same stacks of the present paper.

The results are reported in Table 7 for both nocodazole and taxol-doped cell cultures. If data are not available for a given stack, the table reports NA.

It can be noticed that Algorithm 2 generally yields higher mean velocity if compared to AMicro, whereas Algorithms 3 exhibit similar values in control stacks and at 100 nM drug concentrations. Moreover, Algorithm 3 yields very low standard deviation, due to an algorithmic choice that suppresses variability to a large extent. As for length, AMicro and Algorithm 3 yield nearly the same average values; instead, Algorithm 2 is able to detect much shorter tracks. This can be justified by the removal of shorter tracks carried out by both AMicro and Algorithm 2, but not by Algorithm 3.

If we focus on the comparison between AMicro and Algorithm 4, we can appreciate that they yield coherent trends as for mean velocity and length, with AMicro generally providing slightly lower values, due to the fact that [18] is tuned so as to yield a larger number of detected tracks. We point out the significant difference in average length at 100 nM concentration of both taxol and nocodazole. This is due to the screening process implemented in AMicro, which removes all short tracks not matching the uniform linear motion assumption.

Table 7. Comparison among different algorithms.

*Concentration 80 nM. **Concentration 20 nM.

Figure 7. Nocodazole-doped cells: comparison between AMicro and Algorithms 2, 3 and 4.

Table 8. Average growth velocity: comparison between AMicro and manually scored data.

The velocity and length yielded by the algorithms, in the case of nocodazole-doped cell cultures, are reported in Figure 7. We point out the fact that, whereas AMicro and Algorithm 4, at intermediate concentrations (i.e. 10 nM), shows an increase of MTs dynamicity in terms of velocity, this behaviour is not revealed by Algorithms 2 and 3, which yield a monotonic decreasing velocity curve.

Finally, the results of AMicro have been compared with those computed by hand by expert biologists of the Department of Molecular Biotechnology and Health Sciences of the University of Turin. For the sake of brevity, this has been done on a small subset of the same data stacks, referring to nocodazole-doped cell cultures, and only means velocity values have been taken into account. The available comparisons are listed in Table 8.

Paired t-test has been evaluated, which has turned out to be at the limits of statistical significance (p = 0.052). The Pearson correlation coefficient, adjusted in order to keep into account the small sample size [29] , turned out to be 0.89, so denoting a strong correlation between the two sample measures.

6. Conclusions

We have presented AMicro, an automatic tool for tracking and analyzing astral MTs in fluorescence confocal microscopy images. The algorithm has been validated using data and metrics provided by the International Competition [8] ; then, it has been run on real experimental data. Despite the lack of a ground truth, the validation process has provided encouraging results, which are also well-substantiated by the expected drug effects at high concentrations. An important aspect is related to computational time; indeed the time spent on analyzing samples is almost negligible if compared to the manual labour (several hours compared to few minutes). Moreover, the automatic software is not affected by human errors, due to tiredness or attention deficit, and can provide a valid support for biological experiment evaluation.

The main achievement of our method is that it is extremely easy to use, and all parameters are automatically set up without requiring the user intervention. At present, it is being routinely employed by the biologists of the Department of Molecular Biotechnology and Health Sciences of University of Turin; other experiments, different from those employed in this work, are being carried on. As for future developments, the preliminary results related to the impact on the velocity trend of intermediate concentrations of nocodazole will be refined using cell cultures doped with more nocodazole concentrations (as well as other MTAs). AMicro will be tested also on mitotic cell cultures. Finally, to ensure a better portability, it is planned to leave the MathWorks environment developing an ImageJ plug-in.

Cite this paper: Varrecchia, M. , Levine, J. , Olmo, G. , Grangetto, M. , Gai, M. and Cunto, F. (2019) A Method for Astral Microtubule Tracking in Fluorescence Images of Cells Doped with Taxol and Nocodazole. Advances in Molecular Imaging, 9, 60-86. doi: 10.4236/ami.2019.94009.

[1]   Prosser, S.L. and Pelletier, L. (2017) Mitotic Spindle Assembly in Animal Cells: A fine Balancing Act. Nature Reviews Molecular Cell Biology, 18, 187-201.

[2]   Cooper, G.M. (2000) The Cell: A Molecular Approach, Sinauer Associates, Inc., Sunderland, MA.

[3]   Ito, A. and Goshima, G. (2015) Microcephaly Protein ASP Focuses the Minus Ends of Spindle Microtubules at the Pole and within the Spindle. Journal of Cell Biology, 211, 999-1009.

[4]   White, D., Honoré, S. and Hubert, F. (2017) Exploring the Effect of End-Binding Proteins and Microtubule Targeting Chemotherapy Drugs on Microtubule Dynamic Instability. Journal of Theoretical Biology, 429, 18-34.

[5]   Vonesh, C., Aguet, F., Vonesh, J.L. and Unser, M. (2006) The Colored Revolution of Bioimaging. IEEE Signal Processing Magazine, 23, 20-31.

[6]   E. Meijering, I. Smal, and G. Danuser (2006) Tracking in Molecular Bioimaging. IEEE Signal Processing Magazine, 23, 46-53.

[7]   Matov, A., Applegate, K., Kumar, P., Thoma, C., Krek, W., Danuser, G. and Wittmann, T. (2010) Analysis of Microtubule Dynamic Instability Using a Plus End Growth Marker. Nature Methods, 7, 761-768.

[8]   Chenouard, N., Smal, I., De Chaumont, F., Maska, M., Sbalzarini, I.F., Gong, Y., Cardinale, J., Carthel, C., Coraluppi, S., Winter, M., Cohen, A.R., Godinez, W.J., Rohr, K., Kalaidzidis, Y., Liang, L., Duncan, J., Shen, H., Xu, Y., Magnusson, K.E.G., Jalden, J., Blau, H.M., Paul-Gilloteaux, P., Roudot, P., Kervrann, C., Waharte, F., Tinevez, J.Y., Shorte, S.L., Willemse, J., Celler, K., Van Wezel, G.P., Dan, H.W., Tsai, Y.S., De Solórzano, C.O., Olivo-Marin, J.C. and Meijering, E. (2014) Objective Comparison of Particle Tracking Methods. Nature Methods, 11, 281-289.

[9]   Dzyubachyk, O., Meijering, E. and Smal, I. (2012) Chapter Nine-Methods for Cell and Particle Tracking. In: Michael Conn, P., Ed., Imaging and Spectroscopic Analysis of Living Cells, Methods in Enzymology, Academic Press, London, 183-200.

[10]   Godinez, W.J., Lampe, M., Wörz, S., Müller, B., Eils, R. and Rohr, K. (2009) Deterministic and Probabilistic Approaches for Tracking Virus Particles in Time-Lapse Fluorescence Microscopy Image Sequences. Medical Image Analysis, 13, 325-342.

[11]   Roudot, P., Ding, L., Jaqaman, K., Kervrann, C. and Danuser, G. (2017) Piecewise-Stationary Motion Modeling and Iterative Smoothing to Track Heterogeneous Particle Motions in Dense Environments. IEEE Transactions on Image Processing, 26, 5395-5410.

[12]   Smal, I., Loog, M., Niessen, W. and Meijering, E. (2010) Quantitative Comparison of Spot Detection Methods in Fluorescence Microscopy. IEEE Transactions on Medical Imaging, 29, 282-301.

[13]   Boulanger, J., Kervrann, C., Bouthemy, P., Elbau, P., Sibarita, J. and Salamero, J. (2010) Patch-Based Nonlocal Functional for Denoising Fluorescence Microscopy Image Sequences. IEEE Transactions on Medical Imaging, 29, 442-454.

[14]   Sbalzarini, I.F. and Koumoutsakos, P. (2005) Feature Point Tracking and Trajectory Analysis for Video Imaging in Cell Biology. Journal of Structural Biology, 151, 182-195.

[15]   Smal, I. and Meijering, E. (2015) Quantitative Comparison of Multiframe Data Association Techniques for Particle Tracking in Time-Lapse Fluorescence Microscopy. Medical Image Analysis, 24, 163-189.

[16]   Sironi, L., Solon, J., Conrad, C., Mayer, T.U., Brunner, D. and Ellen-berg, J. (2011) Automatic Quantification of Microtubule Dynamics Enables Screening of New Mitotic Spindle Regulators. Cytoskeleton, 68, 266-278.

[17]   Mahemuti, B., Inoue, D., Kakugo, A. and Konagaya, A. (2016) Investigation of the Microtubule Dynamics with Probabilistic Data Association Filter. 2016 IEEE 11th Annual International Conference on Nano/Micro Engineered and Molecular Systems, Sendai, Japan, 17-20 April 2016, 101-106.

[18]   Levine, J., Grangetto, M., Varrecchia, M. and Olmo, G. (2018) Detection and Tracking of Astral Microtubules in Fluorescence Microscopy Images. 2018 25th IEEE International Conference on Image Processing, Athens, 7-10 October 2018, 361-365.

[19]   Smal, I., Basu, S., Sayas, L., Galjart, N. and Meijering, E. (2017) Gaussian Processes for Trajectory Analysis in Microtubule Tracking Applications. 2017 IEEE 14th International Symposium on Biomedical Imaging, Melbourne, 18-21 April 2017, 206-209.

[20]   Smal, I., Galjart, N. and Meijering, E. (2018) Accurate Estimation of Intracellular Dynamics and Underlying Spatial Structures Using Hierarchical Trajectory Smoothing. 2018 IEEE 15th International Symposium on Biomedical Imaging, Washington DC, 4-7 April 2018, 973-976.

[21]   Samuylov, D.K., Szekely, G. and Paul, G. (2017) Tracking Microtubule Ends Is More Than Point Tracking. 2017 IEEE 14th International Symposium on Biomedical Imaging, Melbourne, 18-21 April 2017, 808-812.

[22]   Newby, J.M., Schaefer, A.M., Lee, P.T., Forest, M.G. and Lai, S.K. (2018) Convolutional Neural Networks Automate Detection for Tracking of Submicron-Scale Particles in 2D and 3D. Proceedings of the National Academy of Sciences of the United States of America, 115, 9026-9031.

[23]   Sorokin, D.V., Peterlk, I., Ulman, V., Svoboda, D., Neasov, T., Morgaenko, K., Eiselleov, L., Tesaov, L. and Maka, M. (2018) Filogen: A Model-Based Generator of Synthetic 3-D Time-Lapse Sequences of Single Motile Cells with Growing and Branching Filopodia. IEEE Transactions on Medical Imaging, 37, 2630-2641.

[24]   Svoboda, D. and Ulman, V. (2017) Mitogen: A Framework for Generating 3D Synthetic Time-Lapse Sequences of Cell Populations in Fluorescence Microscopy. IEEE Transactions on Medical Imaging, 36, 310-321.

[25]   Yao, Y., Smal, I. and Meijering, E. (2018) Deep Neural Networks for Data Association in Particle Tracking. 2018 IEEE 15th International Symposium on Biomedical Imaging, Washington DC, 4-7 April 2018, 458-461.

[26]   Lucey, B.P., Nelson-Rees, W.A. and Hutchins, G.M. (2009) Henrietta Lacks, Hela Cells, and Cell Culture Contamination. Archives of Pathology & Laboratory Medicine, 133, 1463-1467.

[27]   Garg, R. and Kumar, A. (2012) Comparison of Various Noise Removals Using Bayesian Framework. International Journal of Modern Engineering Research, 2, 265-270.

[28]   Applegate, K.T., Besson, S., Matov, A., Bagonis, M.H., Jaqaman, K. and Danuser, G. (2011) Plustiptracker: Quantitative Image Analysis Software for the Measurement of Microtubule Dynamics. Journal of Structural Biology, 176, 168-184.

[29]   Soper, H.E., Young, A.W., Cave, B.M., Lee, A. and Pearson, K. (1917) On the Distribution of the Correlation Coefficient in Small Samples. Biometrika, 11, 328-413.