Main

Single-particle tracking (SPT) is a valuable tool for characterizing protein activity (1–3). While Brownian diffusion is commonly observed inside of cells, deviations from this behavior are of major interest and can provide key biophysical insight (4–6). Classically, non-Brownian behaviors have been identified by fitting a power law to the mean-squared displacements (MSD) as a function of the time lag τ, (MSD = Γ·τ^α) (6), where subdiffusive motion has an α exponent < 1 and superdiffusion has an exponent > 1. While fitting the MSD curve can therefore be used for categorizing motion, the anomalous exponent and the generalized diffusion coefficient provide little insight into the combination of underlying biological forces, e.g., diffusion, directed motion, and confinement.

To accurately model the motion encountered in cellular contexts, physics-based models have been developed to explicitly describe interactions that give rise to nonlinear MSDs. For example, superdiffusion can arise from the linear movement powered by molecular motors (7), and subdiffusive motion can arise from confinement (8, 9). In such models, the observed track positions result from a stochastic process described by hidden (unobserved) variables that evolve with time.

In the context of random processes, Maximum Likelihood Estimation (MLE) is used to determine the underlying parameters of the model. MLE consists in computing the probability of observing a track given the model and its parameters. The main challenge in models with hidden variables is that they require computing the integral of a joint probability density over all possible hidden variables. This integration step is currently realized with coarse-grained approximations or sampling methods that are either slow or inaccurate. Thus, designing a hidden-variable method that benefits from accurate and efficient integration is an important challenge for improving the reliably and ease-of-use of the physics-based models. Further, when the underlying physical model is undecided, statistical tests can be applied to identify the most appropriate model (10–12).

Current hidden-variable models applied to non-Brownian diffusion have three disadvantages: 1) They perform the integration step using either coarse-grained approximations (7) or sampling method making them either inaccurate or slow; 2) they are specialized in either confinement or directed motion but do not treat both types of motion; 3) they typically do not allow variations of the potential well position or the speed of the directed motion, properties that are frequently encountered in biological systems.

Here, we present aTrack, an analysis tool that alleviates the limitations mentioned above. Our approach uses a versatile motion model that considers the relationships between the observed track, the real particle positions (localization error and Brownian motion), as well as the influence of a non-Brownian variable that can either be the potential well for confined diffusion or the velocity vector for directed motion. The main innovations of this model are: 1) it uses analytical recurrence formulas to perform the integration step for complex motion, improving speed and accuracy; 2) it handles both confined and directed motion; 3) anomalous parameters, such as the center of the potential well and the velocity vector are allowed to change through time to better represent tracks with changing directed motion or confinement area; and lastly 4) for a given track or set of tracks, aTrack can determine whether tracks can be statistically categorized as confined or directed, and the parameters that best describe their behavior, for example, diffusion coefficient, radius of confinement, and speed of directed motion.

We validate the approach on simulated data and demonstrate its versatility for analyzing a variety of experimental SPT data, including particle diffusion in an optical trap, detection of motile bacteria with gold nanoparticles, and motion characterization of spindle pole bodies in budding yeast.

Results

Accounting for hidden variables that characterize confined and persistent motion

Modeling non-Brownian motion

We model noisy tracks undergoing confined, Brownian, and directed motion by considering four relations at each time step: (1) there is a Brownian diffusion step followed by (2) an anomalous step (Fig. 1a-b); (3) the hidden anomalous variable, h, can evolve according to a Gaussian distribution; and (4) localization error is incorporated as a Gaussian-distributed noise term added to the underlying real position to produce the observed positions. This model encompasses a variety of motion types depending on the model parameters, as illustrated in Fig. 1c. For example, particles can be immobile, where the observed displacements are only due to localization error; tracks can undergo Brownian motion, as well as anomalous super- and subdiffusion; or multiple motion mechanisms can also occur simultaneously, (e.g. diffusion and drift, or changes in directed motion direction and speed). The model can also account for confinement in an area of variable radius (8, 9, 13) as well as a diffusing potential well. In our model, the velocity is the characteristic parameter of directed motion and the confinement factor represents the force within a potential well. More precisely, the confinement factor l is defined such that at each time step the particle position is updated by l times the distance particle/potential well center (see the Methods section for more details).

Figure 1.

Determining the type of motion

To categorize the type of motion from a measured trajectory, we first calculate the likelihood that a track belongs to each considered motion class (diffusion, directed, or confinement) and then perform a statistical comparison between the likelihoods. To do the former, we must integrate the joint probability density function of a track for all the modeled hidden variables, e.g. all the real positions and the potential well positions or the velocity vectors. Since our model is a multivariate Gaussian process expressed as the product of univariate Gaussian functions, we can perform the integration step using analytical recurrence formulas (see Supplementary information) (14, 15). The recurrence formula enables our model computation time to scale linearly with the number of time points.

We can apply these formulas to determine the best sets of parameters and the maximum likelihoods, then use the ratio between the maximum likelihood assuming Brownian diffusion (null hypothesis) and the maximum likelihood assuming confinement or directed motion as the alternative hypothesis to build a likelihood ratio test (16, 17). Fig. S2 shows that these likelihood ratios, ρ = l_Brownian/l_confined or ρ = l_Brownian/l_directed are systematically skewed towards 1 when particles follow Brownian diffusion. Conversely, ρ is skewed toward zero when applying the directed-diffusion test to directed tracks or the confined-diffusion test to confined tracks. These properties enable the likelihood ratio to serve as a robust proxy for the p-value as it overestimates it. As a consequence, if ρ < X with X the type I error rate (e.g. 0.05), we can reject the null hypothesis with a confidence of 1 − X. Relying on the skewness of the likelihood ratio to obtain an upper bound of the p-value is a simple way to categorize the type of motion, but sometimes a more sensitive test is needed. In such cases, one can use simulations to better estimate the p-value (See Supplemental Information: Statistical test for more details).

To estimate the impact of the track length on the classification, we simulated tracks either in confined diffusion or in linear motion (without diffusion) of varying lengths (5-400 steps for confined, 4-100 steps for directed) and computed the likelihood ratios. As expected, the classification certainty increases with the track length; where the inverse of log ρ increases with the number of steps (Fig 2a). Of course, the statistical certainty depends on the track parameters (Fig 2b). In the range of evaluated anomalous parameters, the significance of the test increased with the confinement factor and the directed motion velocity. To determine the useful range more systematically, we varied the track length and either the confinement factor for confined motion or the velocity for directed motion, and computed the average likelihood ratios (Fig 2c). A low average ratio indicates significantly low p-values for most tracks as this ratio is an overestimate of the p-value. We see that the test is significant as long as the anomalous parameter is high enough or the track length is high enough. Note that increasing the confinement factor so much that the confinement radius becomes similar or smaller than the localization error will impair the capacity of the test.

Figure 2.

Characterizing confinement

To characterize confined trajectories, aTrack estimates several parameters, namely the diffusion coefficient D and the diffusion length where Δt is the single-frame time step; the confinement factor l, which is proportional to the spring constant of the potential well; the diffusion coefficient of the confinement area D_c; and the localization error σ; and The confinement radius, which is proportional to , can also be calculated (See the methods section for more details).

To measure the precision of parameter estimation, we simulated tracks with different confinement factors, l, (Fig. 3a). We then used aTrack to estimate the parameters for each track of 200 time steps (Fig. 3b). The average diffusion coefficient, confinement factor, and confinement radius were accurate over the range of confinement factors, 0 to 1 per time step. Panel Fig. 3c shows the working range for calculating the confinement factor as a function of track length and confinement factor. Longer tracks result in better parameter estimation. Similarly, our method correctly estimates the confinement radius (Fig. S3d). In a second set of simulations, we tested the reliability of our predictions when the confinement area is moving (Fig. 3d, Fig. S3e).

Figure 3.

Characterizing directed motion

Superdiffusive behavior, in particular directed motion, occurs in a variety of circumstances, such as molecular motormediated active transport (7, 18) and polymerase processivity (19). Localization error complicates velocity estimation, especially when the localization error is relatively large. We simulated noisy tracks undergoing linear motion at various speeds (Fig. 4a) with a localization error of 20 nm·Δt⁻¹ and estimated the velocity per tracks. Fig. 4b shows the velocity estimates for tracks with 30 time points. By varying both the track length and the velocity, we found our method to be reliable for a wide range of velocities as long as tracks were long enough (Fig. 4c).

Figure 4.

In real experiments, persistent motion is rarely perfectly linear, that is, changes in direction and speed are very common (7, 18–20). Naturally, characterizing directed motion with direction (orientation) changes is more difficult when localization error is non-negligible (21, 22). To verify our method’s capacity to accurately quantify tracks with such behaviors, we simulated tracks with constant speed and random changes of orientation (rotational diffusion). We previously observed similar behavior when analyzing the directed motion of the Rod complex in bacteria (23), motor-driven directed motion in mammalian cells (21, 24), and cell motility (25). We varied the rate of orientation changes and the track velocity to determine the working range of aTrack for this type of directed motion (Fig. 4d,e,f). In Fig. 4e,f, we find accurate estimates of the velocity and the rotational diffusion for a wide range of parameters. As the velocity increases, we can distinguish directed motion from Brownian motion for an increasing range of rotational diffusion (Fig. 4e). However, we also see that fast changes in orientation (high rotational diffusion) make estimating the velocity and rotational diffusion more difficult. This tradeoff is expected, as rapid changes in direction make the track appear more diffusive as shown in the likelihood ratio heatmap Fig. S5a. The diffusion length heatmap (Fig. S5a) explains why the rotational diffusion coefficient is poorly estimated when high: the model interprets the high rotational diffusion as simple diffusion since the two types of motion are very difficult to distinguish in this regime.

Sometimes, particles undergo diffusion and directed motion simultaneously, for example, particles diffusing in a flowing medium (26). Sample drift can also introduce a combination of diffusion and directed motion in single-molecule tracking. Indeed, the thermal expansion of instrument components like microscope stages can induce steady motion that masks the biologically relevant diffusive motion (27). We first tested that our method could correctly differentiate directed motion from diffusion (Fig. 4g) for a range of diffusion coefficients and track lengths at a fixed velocity of 0.1 µm.Δt⁻¹. For mixed motion, our method accurately estimates the diffusion length and velocities even for short tracks, provided the diffusion length was low compared to the velocity. When the diffusion length is large compared to the velocity, parameters can still be predicted, but longer tracks are needed for reliable estimates.

Characterizing populations of tracks

The amount of information in an individual track is limited by its length making it difficult to extract the parameters precisely. One way to overcome this limitation is to consider a population of tracks that have the same state. The likelihood of the population can be easily computed by multiplying the probabilities of the individual tracks. To test our population approach, we simulated two populations: confined tracks and directed tracks. Fig. S6 shows that the expected linear increase in the log likelihood with the number of tracks and that using more tracks results in more precise parameter estimates.

Populations of tracks can often contain multiple states (e.g., free diffusion and directed motion). By taking into account the fraction of the particle in each state, we built a multi-state population model. We tested the capacity of our multistate-population model on groups of simulated tracks with 300 time points, where each track follows one of the 5 states shown in Fig. 5a.

Figure 5.

The first step is to determine the number of states. To this end, we compute the likelihood of the model depending on the number of states. As expected, the likelihood increases with the number of states, until the correct number of states is reached and the function plateaus, in this case 5 states (Fig. 5b). This plateau in likelihood is a good indicator that the appropriate number of states is reached; however, increasing the number of states further usually results in higher likelihood, albeit marginally.

Quantitative criteria such as the Akaike Information Criterion (AIC) (28) and the Bayesian Information Criterion are often used to determine the number of parameters of a model by placing a trade-off between increasing the likelihood and increasing the number of parameters. When increasing the number of tracks, from 50 to 12,800, we found these criteria to be unreliable (Fig. S7). In theory, these criteria only work if the true underlying model is included in the alternative models, which is never the case for real tracks. As an alternative, we found that adding a small penalization term proportional to the number of parameters and the log likelihood provides a reliable criterion for identifying the number of states for any dataset size, even with mismatches between the data and the model assumptions.

Once we have identified the number of states, the parameters of each state are estimated at the population level (Fig. 5b), even when individual tracks remain difficult to classify. For instance, classifying whether a track is in the diffusive state (orange) or diffusive plus directed state (green) is difficult.

Robustness to model mismatches

One of the most important features of a method is its robustness to deviations from its assumptions. Indeed, experimental tracking data will inevitably not match the model assumptions to some degree, and models need to be resilient to these small deviations. To test the generalizability of our approach to other types of motion, we simulated tracks using a different motion model, namely fractional Brownian motion (29). This was performed with three anomalous diffusion exponents, 0.5, 1.0, and 1.5, corresponding to sub-diffusion, Brownian diffusion, and super-diffusion, respectively Fig. 6a. The performance of aTrack was then measured by computing the difference between the likelihood of the directed-motion model and the confined-motion model. This metric is closely related to the likelihood ratio mentioned earlier, and has the advantage of showing all three motion behaviors. We expect the likelihood difference to be <0 for sub-diffusive motion, near 0 for Brownian diffusion, and >0 for super-diffusion. We verified this by plotting a histogram of the log-likelihood differences for each type of motion and estimated the classification accuracy in detecting anomalous diffusion from Brownian diffusion. aTrack was 88% accurate when differentiating subdiffusive fractional Brownian motion (anomalous diffusion exponent of α = 0.5) from Brownian motion (α = 1), and 95% accurate for differentiating superdiffusive fractional Brownian motion (α = 1.5) from Brownian motion.

Figure 6.

To compare our approach to one of the leading methods specifically designed to characterize fractional Brownian motion, we performed the same test with Randi (30), the best-performing machine learning method in a head-to-head comparison of available techniques (31). On the same fractional Brownian motion dataset, Randi and aTrack achieved similar accuracies.

We then used the same approach to compare aTrack and Randi on tracks generated with our motion model, which differs from Randi’s anomalous motion models (Fig. 6b). Here, we created a dataset of anomalous diffusion with parameters that result in tracks with qualitatively similar properties and MSD curves to those observed with fractional Brownian motion (see tracks in Fig. 6a&b). For these data, our model performs with high accuracy, at least 99%. In contrast, Randi shows a lower classification accuracy. In particular, Randi has difficulty differentiating diffusive tracks from tracks with both diffusion and directed motion, 61% accuracy (only 11% better than random labeling). Curiously, directed versus Brownian tracks were also surprisingly inaccurate (83%), considering the striking difference in track behaviors (see Fig. 6b, diffusive versus directed examples).

To test another type of mismatch between our hidden variable model assumptions, we simulated diffusive tracks confined within rigid boundaries. This differs from our model, which uses a potential well to model confinement. We varied the radius of the rigid boundary for simulated tracks and measured the effect on the estimated confinement radius (Fig. 6c). aTrack accurately determines the confinement radius for a wide range of confinement radii, where the calculated confinement radius is estimated as 3-times the standard deviation of the potential well. The lower radius bound of the operating range relates to the diffusion length per step, while the upper bound is limited by the extent of the particle’s exploration of the confinement area. Note that the exploration distance is determined by the track length and by the diffusion length). This calculation was relatively independent of track length (Fig. 6d).

Motion blur is another experimental effect in single-particle tracking that can bias parameter estimation. While our model does not explicitly account for it, the estimated diffusion coefficient can be easily corrected to adjust for that effect (Fig.S 8). As described by Berglund, static and dynamic localization errors have antagonistic effects on the offset term of the MSD. Our model, which explicitly models static localization error but not dynamic error, yields good estimates of the diffusion length if (MSD curve with a positive offset), but it underestimates d by a factor if . To explicitly adapt our tool to motion blur, one can include motion blur using our new framework for model design (33).

Implementation on experimental data

To test the usefulness of aTrack on experimental data, we performed classification, population characterization, and parameter-estimation experiments.

First, we applied aTrack to analyze the movement of the spindle pole body (SPB) in Saccharomyces cerevisiae (Fig. 7a). The SPB is the microtubule-organizing center in yeast and is embedded in the nuclear envelope. Directed motion of the SPB occurs during S phase, where actin is required to establish spindle orientation by directing the spindle pole body towards the bud neck as well as during mitosis when the spindle elongates to separate the chromosome masses of mother and daughter cell (34). To visualize SPB dynamics, we imaged unsynchronized cells expressing Spc42-mCherry at a time resolution of 100 ms and used treatment with Latrunculin A (Lat A) to disrupt actin polymerization. Computing the MSD of the population of tracks showed that in average tracks appear diffusive (Fig S9a). To go beyond this ensemble metric, we used aTrack to compute the likelihood ratio (l_b/l_d) of each track of 99 time points and to infer the associated p-value from the ratio and the distribution of the likelihood ratio under the null hypothesis (Fig. 7b). Then, we computed the fractions of significantly directed tracks (type I error rate of 5%) (Fig. 7c). In untreated cells, 17 % of the tracks exhibited significant directed motion. In contrast, cells treated with LatA showed significantly lower fractions of directed tracks (10%, p-value = 0.0108). This drop in the fraction of directed tracks was not observed in a LatA-resistant strain. Thus, aTrack can reliably detect directed actin-dependent motion at timescales of 10 s.

Figure 7.

To test the applicability of aTrack to the field of biosensing, we performed a tracking experiment with highly visible gold nanoparticles (AuNP) and E. coli. Specifically engineered gold nanoparticles (AuNPs) can attach to cells, and many species are motile. Zapata-Farfan et al.. Thus, detecting a directed fraction of AuNPs could be used for sensitive and fast biodetection in complex environments. AuNPs were added to a diluted E. coli culture and imaged. While free AuNPs diffuse rapidly, our method identified directed tracks in a 5.9% population with very high certainty (likelihood ratio < 0.001, Fig. 7d). This directed fraction is readily apparent by manual data inspection (see supplemental video). Notably, while the directed fraction also exhibited lateral oscillations, consistent with cell rotation (36), our model was robust to this deviation from the model assumptions. Next, we analyzed the tracks at the population level and computed the number of states as a function of likelihood (Fig. 7e). The likelihood function shows a noticeable difference between 4 and 5 states. In the 5-state model, there are 2 directed states and 3 diffusive states. The diffusive states likely represent free AuNPs as well as NPs conjugated to diffusing debris of different sizes. The 2 directed populations are likely caused by the abrupt tumbling motion present in some tracks. These directed states comprised 6.2% of tracks, similar to the analysis of individual tracks.

Notably, in this dataset, the likelihood plot shows an ambiguous number of states, where the likelihood increases marginally after 5 states. Such ambiguities are expected when the dataset and the model assumptions do not perfectly match. To ensure that the directed fraction is well estimated independently of the number of states, we varied the number of states to see how the parameter estimates changed, as in (37). We found that using 3-7 states, resulted in the same fraction of clearly directed tracks with similar parameters.

Finally, we tested our method’s capability to detect confinement. To do so, we confined 1-micron beads in an optical trap (38), and varied the laser power to control the trap stiffness. Based on populations of 300 tracks of 20 time points, we measured the confinement factor l and confirmed that it increases with the laser power of the optical trap, while the calculated confinement diameter u decreases, where (Fig. 7g).

Discussion

aTrack is a new tool for classifying and characterizing noisy tracks, which performs well in a diverse range of conditions. Specifically, the framework’s flexibility is relevant for a wide variety of diffusive, confined, and directed motion types. Our tool classifies the motion type as Brownian or not and quantifies biologically relevant motion parameters, such as the diffusion coefficient, confinement diameter, and velocities. Importantly, this approach calculates how statistically robust a classification is regarding the likelihood difference from a Brownian diffusion model. The flexibility of our approach in capturing a variety of motion-type behaviors contains the usefulness of the catch-all approach of using an anomalous exponent; however, unlike the anomalous exponent that fits multiple underlying motion models (6), our framework has the advantage of outputting interpretable parameters that describe the motion, e.g. velocity or confinement radii.

aTrack has several features that make it advantageous for analyzing directed and confined motion. For example, allowing the anomalous variable to change over time enables more flexibility. For directed tracks in cells, motion is rarely exactly straight over long distances. Allowing changes in velocity makes it possible to classify and characterize these curved trajectories robustly. Analogously for confined motion, a confined particle may leave its local environment or hop between environments (39–41). Our model accounts for this by allowing the center of the confinement well to diffuse as well. Interestingly, we found that directed motion can be readily identified from very short tracks due to the deterministic nature of this type of motion. Longer tracks are needed to classify confined motion with the same statistical certainty. This is to be expected, as a confined particle should reach the boundary of the confinement area several times to be distinguishable from a freely diffusive particle.

As with all classification tools, a source of ambiguity arises when motion types resemble one another. For example, this confusion occurs between immobile particles and very tightly confined particles, which can sometimes be indistin-guishable and not necessarily insightful. Indeed, immobilized fluorophores can be modeled as confined to an area around a substrate to which they are bound. To better distinguish these ambiguous behaviors, experimental modifications are needed, such as adapting the experimental framerate or improving the localization precision using with brighter fluo-rophores (42, 43).

While machine-learning tools have proven to be effective for subdiffusion characterization (31), our approach shows it is possible to achieve similarly high performance with many fewer parameters. For our tool, these parameters map onto the stochastic physical behaviors of particle motion, making interpretation more straightforward. Compared to the machine-learning tool we tested, we found that aTrack is more robust to small model mismatches. This finding is consistent with the well-documented issues of machine-learning models generalizing poorly to new data. Of course, in the context of tracking, the generalizability of a model to new data is a key factor, as experimentally-obtained data never perfectly match the model assumptions or the training data set. One solution to make a machine learning model that generalizes better is to use a physics-informed neural network (44). Such a network would use probabilistic relationships to efficiently learn the physical properties of unlabeled tracks and would contain far fewer parameters than classical networks.

It is often useful to assume a finite number of states with fixed parameters to model the various molecular states in a sample (43). Selecting the right number of states is difficult, but can be automated by different methods. The criterion, e.g., AIC and BIC, used to determine the number of states is reliable when the model and the data are in perfect agreement. However, experimental data never match the model assumptions perfectly, and even discrepancies between models and simulated data can affect reliability, such as continuous-time simulations and discrete models. This usually results in overestimating the number of states for large data sets (45). To avoid this flaw of classical criteria, we showed that a penalization term proportional to the number of states and to the log likelihood of the data can prevent the spurious increase of the number of states when increasing the number of tracks. The drawback of this approach is the addition of a tunable parameter that influences the number of estimated parameters of the model, and it may still be necessary to limit the number of states to the biologically relevant system, or consider groups of states, e.g. all significantly directed states.

An important limitation of our approach is that it presumes that a given track follows a unique underlying model with fixed parameters. In biological systems, particles often transition from one motion type to another; for example, a diffusive particle can bind to a static substrate or molecular motor (46). In such cases, or in cases of significant mislinkings, our model is not suitable. However, this limitation can be alleviated by implicitly allowing state transitions with a hidden Markov Model (15) or alternatives such as change-point approaches (30, 47, 48), and spatial approaches (49).

In conclusion, we have shown that aTrack can identify anomalous diffusion and parameterize the motion over a broad range of motion types using a robust probabilistic framework. As the motion-model parameters estimated by the method represent physical phenomena, these variables are readily interpretable, for example the diffusion coefficient, confinement radius, and velocities. Finally, the employed motion models were selected to permit analytical integration, which makes calculating the model parameters fast and accurate; of course, this integration strategy can be implemented for a variety of motion models with hidden states, further expanding the applicability of this approach to other motion types.

Methods

Modeling particle motion with observed and hidden variables

Probabilistic model for confined motion

At each time step i, our confined-motion model consists of (1) a Brownian motion step, (2) a confinement step, (3) an update to parameters, and (4) a localization error step. The Brownian motion step updates the particle’s position r_i to an intermediate position z_i. The variable z_i −r_i follows a Gaussian distribution centered at 0 with standard deviation d, where d is the diffusion length, is the diffusion coefficient, and Δt is the time step. Next, the confinement step is modeled by an attractive force between the particle with a potential well centered at h_i. More precisely, the particle moves toward the center of the potential well proportionally to a confinement factor l and to the distance r_i −h_i. The next real position r_i+1 is thus determined by the following relationship r_i+1 = (1 −l) · z_i +l ·h_i. To allow the potential well to move, h is updated at each step such that h_i+1 − h_i follows a Gaussian distribution of mean 0 and standard deviation q. Finally, to model localization error, the observed positions c_i and the real positions r_i are related by the localization precision, σ, where c_i − r_i follows a Gaussian distribution of mean 0 and standard deviation σ.

The distribution of positions for a diffusive particle in a fixed potential well is a Gaussian distribution with standard Deviation . As we have no prior information about the initial center of the potential well, we assume it is positioned according to a Gaussian probability density function centered on the initial observed position. While it would be even better to consider a Gaussian centered around r₀, we simplify it by approximating it to be centered around c₀ and of standard deviation . The joint probability density function corresponding to this model is the product of Gaussian functions shown in equation 1, which is integrated over all hidden parameters to calculate the likelihood, l_{conf ined}.

Where the real positions of the particle r_i, potential well h_i, and intermediate position z_i are hidden variables, and l is the confinement factor. By integrating this joint probability over the hidden variables, we can retrieve the probability of the track (the observed positions) given the model and its parameters.

While this integration step can be computed with a Monte Carlo approach (9), this is computationally expensive. Instead, we integrate using an analytical recurrence formula. This formula is allowed by the fact that the joint probability density function is a product of Gaussians and by the property that for two Gaussians, f and g, with means μ_f, μ_g and standard deviations σ_f, σ_g, the product is also Gaussian, f (x) · g(x) = ϕ · η(x) where ϕ and η are Gaussian distributions described by Equation 2.

See Supplementary information for more details.

Probabilistic model for directed motion

To consider directed motion, we use the same general framework as confinement. At each time step i, the real particle position r_i is first updated to an intermediate position z_i, where the variable z_i −r_i follows a Gaussian distribution centered at 0 with standard deviation d. Next, we add the directed-motion component with a vector w_i, such that the next real position is the sum of the two steps, r_i+1 = z_i + w_i. Combining these two substeps, we get z_i − r_i = r_i+1 − r_i − w_i, simplifying the integration process of the probability density function expressed in equation 3. Analogous to our confinement model, the velocity vector, w_i, is allowed to change over time, where the w_i+1 −w_i follows a Gaussian distribution with mean 0 and standard deviation q. Finally, we include the effect of localization error, where the observed position c_i is related to the real position r_i following a Gaussian distribution, where c_i −r_i is Gaussian distributed with mean 0 and standard deviation σ.

During the first time step, the orientation and length of the directed motion vector (speed) need to be initialized. To do so, we assume w₀ follows a Gaussian distribution function with mean 0 and standard deviation v. The resulting joint probability density function is shown in Equation 3.

As with confinement, the probability of a directed track given the model parameters can be calculated using an analytical recurrence formula to integrate over the hidden positions, r_0→n, and velocities w_0→n−1. The parameter v can, in principle, be used to estimate the velocity of a particle with a constant speed but changing orientation in the imaging plane; however, to estimate the average speed of a particle, we use another metric, k, that appears in our integration process (see Supplementary information).

Modeling time-varying velocities and changes in direction

In our model, each axis, x and y, is treated separately, and the velocities can evolve according to a Gaussian distribution. This allows a particle’s direction to change over time. To quantify how direction changes affect the analysis for particles with a constant speed, we simulated tracks with a fixed speed and time-dependent direction for a range of angular diffusion coefficients, D_θ. The model parameter, q, represents the standard deviation of the change in speed, which can be converted to an angular diffusion length using the following trigonometric relation and the estimated speed of the motion.

In the case of pure directed motion with direction changes, let us consider a single time step, i. We have A, the previous particle position r_i−1; B, the particle position assuming no changes of orientation or diffusion, r_i−1 + v_i−1; and C, the actual particle position after a change of orientation, r_i. ABC forms an isosceles triangle, which can be split into 2 right triangles. We find the following relationship between the scalar distances BC, AB, and , where θ is the orientation angle change.

Fitting method

For the pure Brownian model, the parameters are the diffusion coefficient and the localization error. For the confinement model, the parameters are the diffusion coefficient, the localization error, confinement factor, and the diffusion coefficient of the potential well. For the directed model, the parameters are the diffusion coefficient, the localization error, the initial velocity and the acceleration variance.

These parameters are estimated using the maximum likelihood approach which consists in finding the parameters that maximize the likelihood. We realize this fitting step using gradient descent via a TensorFlow model. All the estimates presented in this article are obtained from a single set of initial parameters to demonstrate that the convergence capacity of aTrack is robust to the initial parameter values.

Multi-state population model

We designed a specific algorithm to retrieve the number of states in a data set and estimate the parameters of each state as described in the code of the script atrack.py available on our Github page https://github.com/FrancoisSimon/aTrack. This multi-state population algorithm starts by performing individual fittings of tracks to get a type of motion and a set of parameters for each track. Then, we use a Gaussian mixture model on the parameters of the individual tracks to provide an overestimate of the number of states (e.g. 20 if we expect 5 states) and fit the model so that every actual state underlying the data is well represented by at least one of the the model state. Next, we iteratively remove the least useful model state and refit the model until we obtain a single-state model. At each iteration, the least useful state is determined as the state with the smallest negative impact on the likelihood.

Track simulations

We subdivided each time step into 20 substeps to simulate approximately continuous tracks. We applied the Brownian diffusion and anomalous movements to each of these substeps. For confined motion, the latter step moves the particle toward the center of the potential well center proportionally to the distance multiplied by a scaled confinement factor, l/20. In simulations with a moving potential well, the center moves according to the well’s diffusion coefficient, q. For directed motion, we applied a shift of constant velocity at each time step after the diffusion step. The orientation of the directed motion could vary according to a rotational diffusion coefficient. The particle’s position after the 20 substeps was set to the particle’s real position, r_i, and the localization error was added to create the observed position c_i. This process is equivalent for undivided frames for directed motion, as the diffusion and the directed motion steps are independent. However, for the confined model, where the diffusion of the particle influences the anomalous step, thus there is a dependence between the diffusion and confinement steps.

Fractional Brownian motion was simulated using the Python package ‘fbm’ (https://pypi.org/project/fbm/) with the Davies and Harte method (50).

Experimental methods

Spindle-pole body

Yeast cells expressing Spc42-mCherry with or without allelic mutation in actin act1-113 (51), (strains KWY10722 and KWY10328 described in (52)) were grown to exponential growth phase in synthetic complete medium with glucose (SCD). Cells were treated with 0.2 mM Latrunculin A (Enzo Life Sciences, BML-T119-0500) and DMSO as solvent control for 10 min prior to imaging. Matrical 384-well glass bottom plates coated with Concanavalin A were used to image treated cells on a temperature-controlled inverted Nipkow spinning disk microscope equipped with the Yokogawa Confocal Scanner Unit CSU-W1-T2 controlled by the VisiVIEW Software (Visitron). It was used in spinning disk mode with a pinhole diameter of 50 µm combined with a 1.45 NA, 100x objective. Images were acquired on an EMCCD Andor iXon Ultra camera (1024×1024 pixel, 13×13um pixel size); for one of three biological replicates, the data was acquired with dual camera settings. Imaging was performed at 30 °C with 80 % laser intensity of a Diode 561 nm, 200 mW laser. Timelapse data were acquired with 100 ms exposure time in stream mode for 300 frames. For all movies, tracks were obtained using the ImageJ (53) plugin TrackMate (54). Peak detection: LoG with radius = 0.45 µm and threshold = 20, linkage: simple lap tracker with a linking maxiumum distance of 0.5 µm and allowing gaps of one frame. The fractions of directed tracks were inferred by computing the likelihood ratio for each track of 99 time points. Then, we simulated Brownian tracks of diffusion length 0.1 µm which corresponds to the diffusion length of the most diffusive tracks (99 percentile) to estimate the distribution of the likelihood ratio of tracks following the null hypothesis. This is a conservative estimate as we found that higher diffusion lengths (compared to the localization error) produce likelihood ratio distributions that are less skewed toward 1. Based on these simulations, we found that tracks with a likelihood ratio lower than 0.295 are significantly directed with an error rate of 5%.

Bacteria detection and tracking

Bacteria detection with nanoparticles was performed following the protocol described in Zapata-Farfan et al. (2023). In brief, E. coli (strain 25922, ATCC) were cultured for 24 h in Trypticase soy broth (TSB) at 37°C and then incubated for 30 min with 100 nm spherical gold nanoparticles at 50 µg/mL (A11-100-CIT, Nanopartz) in 1x PBS. Samples were placed in a small chamber consisting of double-sided tape (5 µm thickness, Nitto) between a coverslip and glass slide. Trackings was performed using Trackmate 7.0 (54) with the LoG detection method with a diameter of 0.70 µm (6 pixels) and quality threshold of 24.

Optical tweezers

The tracking of microspheres captured in an optical trap was achieved using a custom instrument constructed on a modular inverted microscope (MIM/RAMM, ASI Imaging), incorporating optical trapping and imaging paths. The trapping diode laser (SNP-06E-100, Teem) emitting at 1064nm (TEM₀₀) had an average output power of 60 milliwatts. The laser beam is expanded using a 1:3 telescope (Achromatic doublets, Thorlabs) to overfill the objective’s back focal plane. This expanded beam was focused through a 100X, 1.45 NA objective (UPlanXApo 100X, Olympus) to trap dielectric particles. The same objective was used to collect bright-field illumination, which was imaged using a CMOS camera (ORCA - Flash4.0 LT3, Hamamatsu). A 3-axis piezo-driven stage (MicroScan SCXYZ100, Thorlabs) with a precision of 25 nm facilitated sample movement to capture particles in the trap. One µm polystyrene microbeads (Monodisperse fluorescent microspheres, Cromtech Research Center) were suspended in water. The solution was squeezed between two coverslips (No. 1.5, Thermo), and a single bead was brought into the laser trap by translating the sample. Data acquisition was performed at 315 frames per second for 2 minutes, these movies were then analyzed in shorter, 100-frame segments. Beads positions were tracked using TrackMate (54). aTrack analysis was performed on 300 tracks of 20 time points, inputting a fixed localization error and diffusion coefficient, which were determined from the lowest laser intensity (15mW) σ = 0.0116µm, d = 0.0259µm, respectively. Confinement .

Supplementary information

Statistical test

Statistical tests are important as they allow one model to be chosen over another. In the case of independent and identically distributed variables, the log likelihood ratio of a dataset drawn from the null hypothesis follows a chi-squared distribution that can be used to determine the p-value. The p-value can then be used to reject or not the null hypothesis.

In this article, we use the likelihood ratio as our testing metrics. However, the observed variables of our sequential models are not independent. Therefore, the log likelihood ratio is not expected to have a chi-squared form with a fixed degree of freedom. Nevertheless, several methods can still be used to estimate the p-value. Intuitively, the likelihood ratio tells how much more likely one model is compared to another, where the model with the highest likelihood should be preferred. However, this metric alone can be subject to over-fitting. In our framework, the Brownian diffusion model is a specific subset of both the directed and confined motion models. As a consequence, the likelihood of the directed or confined motion model is equal or higher than the likelihood ratio of our Brownian motion model, confining the ratio Brownian/Anomalous between 0 and 1. In Fig. S2, we tested our metrics on several data sets and we found that the likelihood ratio is indeed skewed towards 0 when the directed (resp. confined) test is applied to significantly directed (resp. confined) tracks. On the contrary, when the directed (resp. confined) test is applied to Brownian or confined (resp. directed) tracks, the distribution is consistently skewed toward 1. In comparison, a p-value is by definition uniformly distributed when the data are drawn from the null hypothesis. Therefore, our metrics can be effectively used as an overestimate of the underlying p-value. This means that for a given user-defined threshold representing the rate of false positives, the user can reject the null hypothesis for a given type I error rate. In heatmaps of the likelihood ratio like in Fig. 2c, we show the average likelihood ratio for a range of parameters. In these graphs, a low average ratio indicates a high test sensitivity.

While the likelihood ratio provides a reliable upper bound for the p-value, it is possible to create a more sensitive test by computing a better estimate of the p-value. Considering a track of n time points, running the likelihood ratio test results in a Brownian null hypothesis model of estimated parameters , an alternative hypothesis model of estimated parameters θ^* and a likelihood ratio ρ. By simulating Brownian tracks of length n and parameters , we can compute the empirical distribution of the likelihood ratio under the null hypothesis ρ₀. The p-value associated with a track of likelihood ratio ρ can then be estimated as the quantile ρ of the distribution of ρ_H0.

This statistical framework can be extended to other types of test: We can for example perform a conformity test by computing the ratio of likelihoods assuming a model with some fixed parameters over the likelihood of the model with fitted parameters. If the ratio is close from 0, we can reject the null hypothesis with fixed parameters.

Derivations

Likelihood calculation for confined motion

To simplify the integration of the product of Gaussians described in the Method section, Eq. (1), f_σ((1−l)z_i +lh_i −c_i+1) f_d ((1− l)z_i + lh_i − zi₊₁) can be refactored to of variances k² = σ² + d² and g^′2 = σ²d²/(σ² + d²). This can be further simplified into with and G represent Gaussian distributions of mean 0 with variances k² = σ² + d² and ,respectively. Then, we use the annotation N(x, V) to refer to other Gaussian probability density functions of mean 0 and of variance V.

Eq. (1) is equivalent to:

During the integration process, we obtain the following set of recurrence formulas:

- Initialization:

- Recurrence:

-Final step:

Likelihood calculation for directed motion

Similarly to the confined motion formula, the track probability can be expressed using a recurrence formula that depends on the parameters of the directed motion model, the localization error σ, the diffusion length d, the initial velocity v and the speed of the velocity change q.

- Initialization:

- Recurrence:

-Final steps (steps i = n − 2 and i = n − 1):

Supplementary Figure 1.

Supplementary Figure 2.

Supplementary Figure 3.

Supplementary Figure 4.

Supplementary Figure 5.

Supplementary Figure 6.

Supplementary Figure 7.

Supplementary Figure 8.

Supplementary Figure 9.

Data and code availability

Tracking data and the aTrack software is available online. aTrack is available as a stand-alone software for Windows and as a python package. Installation instructions are provided on the Github page https://github.com/FrancoisSimon/aTrack.

Acknowledgements

The authors thank Sven van Teeffelen for helpful discussions. This work was supported by the Natural Sciences and Engineering Research Council of Canada [NSERC Discovery grant to MM, (RGPIN-06404-2016) to CB, (RGPIN-2022-05142) to LEW], the Canada First Research Excellence Fund (TransMedTech Institute), the ETH research grant ETH-33 19-1 and the Swiss National Science Foundation (project number 320030-236124) to ED.

Additional information

Author contributions

Conceptualization: FS, LEW; Methodology: FS, LEW; Software: FS, IF; Validation: FS; Formal analysis: FS; Investigation: FS, GR, JZ, SP; Resources: MM, CB, ED, LEW; Data curation: FS, LEW; Writing - original draft: FS, LEW; Writing - review and editing: FS, GR, SP, ED, LEW; Visualization: FS, JZ, GR, ED; Supervision: LEW; Project administration: LEW; Funding acquisition: MM, CB, ED, LEW.

Funding

Natural Sciences and Engineering Research Council (RGPIN-06404-2016)

Natural Sciences and Engineering Research Council (RGPIN-2022-05142)

Swiss National Science Foundation (320030-236124)