ExTrack characterizes transition kinetics and diffusion in noisy single-particle tracks

Studying the motion of proteins by single-particle tracking (SPT) allows to characterize how proteins dynamically interact with their environment (Manley et al., 2008; Kusumi et al., 2014). Notably, SPT can reveal if and where proteins are diffusive or immobile (Persson et al., 2013; Uphoff, 2016; Martens et al., 2019). This information has significantly improved our understanding of important biological processes such as transcription-factor binding dynamics, antibody recognition, cytoskeletal dynamics, or intracellular transport (Kusumi et al., 2014; Pierobon et al., 2009; Monnier et al., 2015; Stracy et al., 2016; Jézéquel et al., 2018; Callegari et al., 2019; Li and Xie, 2011; Izeddin et al., 2014; Özbaykal et al., 2020).

Molecules often transition between different motion states. If transitions happen rarely and if trajectories are long, different states such as immobile or diffusive states are reliably detected from time-averaged quantities such as the mean-squared displacement (Michalet, 2010; Elliott et al., 2011; Bosch et al., 2014; Hansen et al., 2018). However, molecules often undergo rapid transitions between different states (Martens et al., 2019; Callegari et al., 2019; Pierobon et al., 2009). Furthermore, tracks are often short as particles can bleach or diffuse out of the field of view or focal plane (Hansen et al., 2018). In such situations, probabilistic methods are better suited to determine global parameters such as diffusion coefficients and transition rates (Dempster et al., 1977; Persson et al., 2013; Meent et al., 2013; Calderon, 2014; Slator et al., 2015; Monnier et al., 2015; Smith et al., 2019; Falcao and Coombs, 2020; Karslake et al., 2021; Vink et al., 2020a; Rahm et al., 2021). Some of these methods can also predict the motion states of individual molecules at every time point (Persson et al., 2013; Monnier et al., 2015; Briane et al., 2020; Rahm et al., 2021), which can reveal the locations of binding sites, spatial correlations, and complex, potentially non-Markovian dynamics (Mahmutovic et al., 2012).

Previous probabilistic methods for diffusive models shown to correctly estimate diffusion and transition parameters (Persson et al., 2013; Vink et al., 2020a) are based on absolute distances between subsequent localizations. These methods have been developed for situations where physical displacements are large in comparison with the localization uncertainty for each molecule. However, when molecules transition rapidly between states, high time resolution is needed, which results in small physical displacements, which, in turn, make identifying different motion states hard or impossible (Fig. 1, a–c). On the contrary, the whole track still allows the distinction of states (Fig. 1 a), simply because subsequent positions of immobile or slowly diffusing molecules fall in the same small area determined by localization error, while subsequent positions of fast-diffusing molecules are nearly uncorrelated.

To account for those spatial correlations, the full sequence of track positions must be taken into account. This approach has been used to characterize a single population of diffusing molecules (Berglund, 2010; Relich et al., 2016). However, if molecules transition between states, this approach becomes computationally demanding because all possible sequences of single-molecule states need to be considered. To avoid this computational complexity, different mean-field approximations (Calderon, 2014; Slator et al., 2015; Bernstein and Fricks, 2016; Lindén et al., 2017; Lindén and Elf, 2018) and machine-learning approaches (Kowalek et al., 2019; Chen et al., 2021) have been proposed. However, their performance across model parameters remains to be investigated.

Here, we propose an alternative probabilistic method to extract diffusive motion states and transitions: We tackle the combinatorial problem of different motion states by introducing a sliding window approximation that maintains the most important spatio-temporal correlations. The method is fast and accurate for a large range of parameters, even if physical displacements are similar to the localization error. The method is also robust with respect to deviations between data and model assumptions. Additionally, the method annotates the state probabilities at the single-molecule level, refines localizations (Lindén and Elf, 2018), and extracts distributions of state durations. We demonstrated its versatility by analyzing two bacterial membrane proteins that diffuse slowly and transition rapidly between immobile and diffusive states.

We developed ExTrack, a maximum-likelihood estimation (MLE) method that contains two main modules: A fitting module fits a multi-state Markovian diffusion model to a data set of noisy single-molecule tracks. This module infers global model parameters including localization error, diffusion lengths, transition rates, and the initial fractions of molecules (at the first time point of all tracks). Part of these global parameters can also be provided by the user, and localization error can even be provided for each peak (Thompson et al., 2002; Quan et al., 2010) if desired. ExTrack is flexible in terms of the number of states and spatial dimensions. Additionally, it can explicitly consider molecules leaving the field of view, which otherwise introduces bias (Hansen et al., 2018). Based on global parameters, a single-molecule annotation module then estimates state probabilities for molecules to reside in each state at each time point. To characterize single-molecule tracks further, we developed two additional modules: A position-refinement module refines molecule positions by taking advantage of spatial correlations between subsequent localizations, conceptually similar to Lindén and Elf (2018). This feature allows to maintain high temporal resolution for state transitions, while attaining high spatial resolution for immobile molecules. A fourth module produces histograms of state durations to reveal potential non-Markovian transitions.

ExTrack is based on a Hidden Markov Model (HMM) that approximates a continuous-time process by a discrete-time Markov model (Das et al., 2009). The method calculates the probability density of observing each track given a set of global parameters (Lindén et al., 2017). In principle, this calculation requires to integrate the joint probability density of states, real positions, and observed positions over all possible sequences of hidden diffusive states and over all sequences of hidden real (physical) molecule positions (see Materials and methods, ExTrack fitting module). This integration can be performed analytically by computing the sequences μ and s (Eq. 6, Materials and methods). However, this problem is computationally intractable for long tracks: Indeed, the computation time scales with the number of possible sequences of states which increase exponentially with the number of time points. To reduce the computation time, we took advantage of the fact that the real position at a given time step is little influenced by the actual state several time points away. This allows us to introduce a sliding-window approximation that greatly limits the computation time of ExTrack. More specifically, when recursively computing the probability density of observing a track, we need to consider all possible combinations of states up to the current time point of the integration. To save computational time, the window method averages over different sequences of states outside of the sliding window, while explicitly considering different combinations of states within the window (see Materials and methods, Using a window...). Choosing the window length is a compromise between accuracy and speed. We suggest using a window length of 3 to 7 depending on the expected diffusion lengths. We will show later that a window length of 3 can be sufficiently large if the ratio between diffusion length and localization precision is greater than about 2.

In many HMMs, it is assumed that state transitions can only occur at the time points of the measurement (Das et al., 2009). However, this approximation introduces a bias toward higher fractions of fast diffusing molecules. Instead, we assume state transitions to occur at the middle of steps (see Materials and methods, Parameter fitting...). Additionally, ExTrack can consider sub-steps to further reduce bias at high transition rates.

ExTrack is available both as a Python library (Simon, 2022) and as a TrackMate module (Ershov et al., 2022) on Fiji.

First, we tested the performance of the ExTrack fitting module by applying it to computationally simulated noisy tracks of molecules (10,000 tracks of 10 positions each, if not stated otherwise) that transition between an immobile state (state 0) and a slowly diffusive state (state 1). The latter is characterized by a small diffusion length d₁ = 2 σ, where σ is the localization error (Fig. 1 a). The diffusion length is the typical physical displacement: $d1=2D1Δt$ with D₁ the diffusion coefficient and Δt the time step. Here, we assume symmetric binding and unbinding rates k_u = k_b = 0.1 Δt⁻¹. Thus, on average, molecules reside in each state for 10 time steps.

The dimensionless parameter $d1*=d1/σ$ can be regarded as a signal-to-noise ratio. For a typical experiment, with σ = 20 nm and a time step of Δt = 20 ms, a rescaled diffusion length of d* = 2 corresponds to a diffusion coefficient of D₁ = 0.04 µm² s⁻¹, which is representative of typical membrane proteins in vivo (Kumar et al., 2010; Oswald et al., 2016).

ExTrack reliably estimates all global model parameters (Fig. 1 d) despite similar observed distances for immobile and diffusive molecules (Fig. 1 c) and despite a low number of 10 localizations per track. These estimates are robust with respect to the initial parameters (for example, if varying the initial d^* between 0.3 and 10). Since molecules are considered at steady state in our example, the initial fraction of immobile molecules is given by $F0=kb ku+ kb.$

Next, we simulated tracks for a three-state model, with $d0*$ = 0, $d1*$ = 2, and $d2*$ = 5, where transition rates are k_i,j = 0.1 Δt⁻¹ for all pairs of states. Fig. 1 e demonstrates that ExTrack estimates global model parameters reliably. If the data contain long enough tracks, ExTrack can also correctly predict two immobile states of different lifetimes and their associated transition rates (Table S1). We will revisit more complex data sets below.

Returning to the simpler two-state model, ExTrack is capable of predicting global model parameters reliably for a large range of model parameters (Fig. 1 f). Predictions are accurate for diffusion lengths as low as the localization error and transition rates as high as 0.5 Δt⁻¹ (for independent variations of k_b and k_u, see Fig. S1 a). To account for rapid transitions and low diffusion lengths, we employed ExTrack considering two sub-steps and a window length of 10. However, the method predicts parameters almost equally reliably without sub-steps (Fig. S1), or with a smaller window size of 3 (Fig. S1 c) while achieving improved computation time (Fig. S1 d).

Next, we compared ExTrack with the two MLE-based methods, vbSPT (Persson et al., 2013) and anaDDA (Vink et al., 2020b) that use absolute distances between localizations for parameter estimation. While vbSPT uses an HMM for the likelihood estimate (Persson et al., 2013), anaDDA is based on an analytical form of the distributions of apparent diffusion coefficients from short tracks (Vink et al., 2020b). Both methods are restricted to a smaller parameter range than ExTrack (Fig. 1, f and g; and Fig. S2 a) in the tested regime. The errors of parameter estimation by vbSPT are largely due to systematic bias, while the error of anaDDA is predominantly stochastic (Fig. S2 d; Vink et al., 2020b). We also tested a mean-field approximation based on track positions and considering hidden particle positions, the variational method UncertainSPT (Lindén et al., 2017). We found that UncertainSPT performs worse and takes more computation time than ExTrack, anaDDA, or vbSPT (Fig. S2 b).

Single-molecule tracks in real cells often deviate from our basic model assumptions. Here, we investigated three different types of such deviations: (i) variations of diffusion coefficients or localization precision, (ii) finite track lengths due to a finite field of view or focal depth, and (iii) physical confinement.

Diffusion coefficients can show intra- or inter-track variations (Stylianidou et al., 2014; Slator et al., 2015; El Beheiry et al., 2015), for example due to local variations of viscosity (Stylianidou et al., 2014), and localization error can vary, for example if molecules are out of focus. We thus simulated tracks of a two-state model with diffusion coefficients or localization precision drawn from a chi-squared distribution with fixed mean and variable coefficient of variation. First, we show that ExTrack gives very accurate predictions when localization error is specified for each peak instead of being treated as a single global fitting parameter (Fig. S3 a). However, even when no prior information on localization error is given, ExTrack reliably predicts the average model parameters for variations up to 30–50% (Fig. 2 a), in contrast to the distance-based methods anaDDA and vbSPT (Fig. S3 c). Track-to-track variations in diffusion coefficient of similar magnitude (up to about 50%) also do not affect predictions of average parameters (Fig. 2 b).

In situations, where the diffusion coefficient is even more broadly distributed, ExTrack can be used assuming a three-state model followed by aggregation of two diffusive states (Fig. S3 b). We tested this aggregation approach with simulations of one immobile and five diffusive states, mimicking a broad distribution of diffusion lengths and jump distances (Fig. S3 b and Table S2). The aggregated three-state approach reliably quantifies transitions between aggregated states and corresponding state fractions, thus providing a practical approach to the often-encountered difficulty of choosing the right number of diffusive states.

Second, molecules can leave the field of view depending on microscopy modality and substrate geometry. For example, cytoplasmic molecules studied by confocal or epi-fluorescence microscopy diffuse in and out of the focal plane, and proteins embedded or attached to a cylindrical membrane (for example, in bacteria) studied by Total Internal Reflection Fluorescence (TIRF) microscopy leave the illumination field (Fig. 2, c and d). Thus, immobile or slowly diffusing molecules are over-represented among long tracks, which have previously been described as “defocalization bias” (Hansen et al., 2018). We alleviate this bias by taking track termination into account explicitly (Materials and methods, Extension of ExTrack...) similarly to previous approaches (Kues and Kubitscheck, 2002; Hansen et al., 2018). In both free 3D diffusion and diffusion along a cylindrical membrane, ExTrack reliably estimates model parameters as long as the typical dimension (focal depth or width of the field of view) is at least twice the diffusion length (Fig. 2, c and d).

Finally, we tested the ability of ExTrack to analyze tracks of spatially confined molecules, as frequently found in membrane domains or small volumes such as bacteria or intracellular compartments. ExTrack performed robustly as long as confining dimensions are at least two to four times larger than the diffusion length (Fig. 2 e). Interestingly, while the diffusion length is underestimated for confining dimensions smaller than about 4 d, the transition rates are predicted reliably even if the confining box size is as small as 2 d.

Next, we tested the performance of the single-molecule probabilistic annotation module of ExTrack, which is based on global model parameters and annotates state probabilities for every time point (Bernstein and Fricks, 2016).

Fig. 3, a and b, shows tracks from the simulation of a two-state model with an immobile and a slowly diffusive state (d^* = 2). Despite the small value of d^*, motion states are reliably estimated in this example. Then, we estimated the predictive power of the probabilistic annotation module by applying it to large sets of tracks with different values of d^* (0.5, 1, or 2). As expected, higher d^* values result in higher confidence in predicting states (Fig. 3 c). To demonstrate the accuracy of the probabilistic annotation, we confirmed that among all molecules predicted to reside in the diffusive state with probability p₁, the fraction of molecules actually diffusive also equals p₁ (Fig. S4 a).

Previous methods often classify molecules categorically into the most likely state (Forney, 1973; Persson et al., 2013; Lindén et al., 2017). We used this approach to measure the performance of ExTrack and compare it to previous methods: First, we estimated the accuracy of categorical state annotations depending on the diffusion length and the transition rates (Fig. S4 b). Increasing transition rates led to worse estimates because states are easier to estimate at time points that are distant from transitions (Fig. S4 c). Second, we found ExTrack annotations to be robust with respect to wrongly chosen global parameters (Fig. S4 d). Finally, we found that the ExTrack state annotation module performs better than vbSPT (Fig. S4 e).

Next, we tested the capacity of ExTrack to refine positions by calculating the most likely physical position for each time point. Fig. 3, d and e, demonstrates that the position-refinement module effectively reduces the localization error of immobile molecules by $N,$ where N is the number of localizations in the immobile segment. This feature allows to obtain accurate positions of molecular binding sites inside cells, while still resolving state transitions dynamics.

ExTrack provides a histogram module that generates probability distributions of state durations. Instead of considering only the most likely set of states, ExTrack considers a large number of potential state vectors and weighs them with their corresponding probabilities. To test the histogram module, we first simulated a two-state model with Markovian transitions. The predicted diffusive and immobile state durations are distributed exponentially, as expected, and in agreement with the simulated data (Fig. 3 f). Therefore, any deviation from exponential decay can reveal more complex transition behavior: As an example, we simulated molecules that transition between two immobile states and one diffusive state (Fig. 3 g). The histogram of immobile state durations then accurately reveals two sub-populations, even though ExTrack considers a two-state HMM model. Our approach thus indicates the presence of a third state, as confirmed by the exponential distributions of state durations after fitting the data to a three-state model (Fig. 3 g).

The histogram approach is also relevant when the transitions are non-Markovian, for example if transition rates are spatially dependent (Mahmutovic et al., 2012; Laurent et al., 2019) or if states have minimum durations. In summary, the histogram module can help identify hidden states or states with non-Markovian transitions and thus guide model choice.

To test our approach on experimental data, we used TIRF microscopy to track single GFP (monomeric super-folder-GFP) fusions to two bacterial membrane proteins in Escherichia coli, each involved in one of the two major pathways of cell-wall synthesis.

First, we studied the cell-wall-inserting penicillin-binding protein PBP1b (Video 1), which was previously described to reside in immobile or diffusive states (Cho et al., 2016; Vigouroux et al., 2020). However, transition rates and potentially hidden states remain unknown. When assuming a two-state model, ExTrack indeed reveals an immobile and a diffusive fraction (Fig. 4 a), with the immobile fraction increasing with decreasing expression level (Fig. 4 b) as expected (Vigouroux et al., 2020). However, distributions of state durations obtained through the histogram module suggest the presence of at least two immobile populations with distinct unbinding rates (Fig. 4 c). Since applying ExTrack assuming a three-state model revealed one immobile and two diffusive states (rather than two immobile states, Fig. S5 a), we also applied ExTrack with a four-state model (Fig. 4, d and e; and Fig. S5 a). The four-state model confirmed two diffusive states and two immobile states: Among the immobile states we found a long-lived state (lifetime of around 0.5 s) that is highly dependent on expression level (Fig. 4 e), likely reflecting enzymatically active PBP1b, and a short-lived state with a lifetime of about 50 ms, likely reflecting non-specific associations with the cell wall. Thus, PBP1b displays rapid transitions between at least four different states.

Next, we investigated the motion of RodZ, a trans-membrane protein that physically links cytoplasmic MreB-actin filaments to a multi-enzyme complex that inserts new peptidoglycan while continuously moving around the cell circumference over minutes (Van Teeffelen et al., 2011; Morgenstein et al., 2015; Van Den Ent et al., 2010). Here, we studied the motion of GFP-RodZ on short time scales of seconds, where continuously moving complexes appear as immobile (Video 2). Assuming a two-state model, the fitting module reveals that 75% of RodZ molecules are immobile (Fig. S5 b based on rates), with a lifetime of about 0.7 s. This time scale is much shorter than the minute-long lifetime of the rod complex (Özbaykal et al., 2020; Cho et al., 2016; Van Teeffelen et al., 2011) demonstrating that a majority of immobile RodZ molecules is not stably associated. Instead, these molecules might transiently bind the MreB-actin cytoskeleton. Interestingly, RodZ molecules seem to often unbind and rebind in very close vicinity from the initial binding site (Fig. 4 f). Such behavior would be expected if RodZ could bind anywhere along extended MreB filaments, since filaments constrain diffusion in two dimensions (Fig. 4 g). To test whether proximal rebinding occurs more often than randomly, we compared tracks that were initially bound, then diffusive for four steps, and then either rebound or remained diffusive (Fig. 4 g). Short distances were indeed over-represented among rebinding molecules compared to molecules that remained diffusive (Fig. 4 h). This behavior contrasts with PBP1b, which appears to bind to random sites (Fig. S5 e). The annotation module of ExTrack thus allows us to identify spatial patterns of molecule binding that can be responsible for non-Markovian binding (Mahmutovic et al., 2012).

In summary, ExTrack provides a suite of robust tools to characterize single-particle tracks, extracting global model parameters, state probabilities at every time point, refined positions, and histograms of state durations, even if tracks are noisy, transitions are rapid, and tracks deviate from idealized model assumptions.

In tracking experiments, a major challenge is to identify the relevant number of immobile and diffusive states. Multiple previous methods obtain this number automatically (Persson et al., 2013; Meent et al., 2013; Smith et al., 2019; Lindén and Elf, 2018; Karslake et al., 2021). However, at least some of these approaches have been reported to overestimate the number of states (Lindén and Elf, 2018; Vink et al., 2020a; Smith et al., 2019). In more recent approaches, a high number of states are fixed, followed by aggregation into one aggregated immobile state and one aggregated diffusive state, based on a user-defined diffusion-coefficient threshold (Metelev et al., 2022). Here, we propose an alternative and iterative approach to complex tracking data: Data are initially fit to a coarse-grained two- or three-state model that can subsequently be expanded depending on desired variables and fitting results. For example, if one is predominantly interested in the exchange between immobile and diffusive molecules but not in the presence of multiple diffusive states, we propose a coarse-grained two-state or an aggregated three-state model that reliably predicts immobile-diffusive transitions, even if diffusion coefficient is variable or if molecules transition between different diffusive states (Fig. 2 b and Fig. S3 b). At the same time, ExTrack can also distinguish multiple diffusive states explicitly (Fig. 1 e). Additionally, the distribution of immobile state durations can reveal the presence of multiple immobile fractions, which can then motivate the increase of the number of states.

The capacity of ExTrack to work with noisy single-molecule tracks is based on the explicit consideration of all sequences of states within a sliding window when computing the probability of every track, while states outside the sliding window are taken into account through averaging to limit computation time. In the future, this versatile principle can be extended to capture different and more complex dynamics, for example by considering persistent motion (Pierobon et al., 2009; Monnier et al., 2015), anomalous diffusion (Chen et al., 2021), and spatial maps of diffusion coefficients or states (El Beheiry et al., 2015).

In the following sections, we will describe the ExTrack method with its four different modules: the fitting module, the annotation module, the position-refinement module, and the histogram module. Subsequently, we will describe implementation and computational time, the generation of computationally simulated data sets, the interpretation of results by vbSPT, and the experimental methods.

ExTrack fits a multi-state diffusion model to noisy single-particle tracking data. We assume that tracks come about according to a continuous-time Markov model, where molecules transition randomly between N diffusive states at rates k_i. As long as molecules reside in state i, they undergo Brownian diffusion with diffusion coefficient D_i. Additionally, observed positions c_i are displaced from real positions r_i according to a Gaussian distribution f_σ(k_i – r_i), where the SD equals the localization error σ. Here and in the following, f_x(y) generally denotes a Gaussian distribution of SD x. The N-state diffusion model is thus characterized by the parameters θ = (σ, D_i, F_i, k_i,j) for all states i, j ∈ (1, …, N). Here, F_i are the fractions of molecules residing in state i at the first position of the track. Later, we will introduce additional parameters for additional spatial dimensions and for the treatment of non-constant track lengths.

Parameters are estimated based on an MLE approach, which, in turn, is based on accurately computing f_C(C | θ), the probability density of observing a track of positions C = (c₁, c₂, ..., c_n). The likelihood of the parameters given the data (θ | all C) then equals the product of f_C(C | θ) for all tracks C. By maximizing this function, we can find θ^∗, the optimal estimator of the underlying parameters. Optimal parameters θ^∗ are found by MLE using the Powell method. ExTrack also allows to fix individual or multiple parameters. This generally speeds up the fitting process and reduces variations in the remaining parameters. In this realm, we also found that fixing the localization error to a slightly wrong value has little impact on the fitting of the other parameters as long as it does not deviate by more than about 20–30%. Here and in the following, we treat localization error as a model parameter, but the user can also provide spot-specific localization errors based on photon counts (Thompson et al., 2002).

In the following sections, we will first compute f_C. This calculation is presented in one spatial dimension (1D). However, the model is easily extendable to 2D or 3D due to the independence of the displacements and localization error in each axis, as we will see below.

This assumption effectively decreases the bias inherent to the discrete approximation of continuous tracks (Fig. S2 c). Later, we will also introduce sub-steps between time points that improve the approximation (see Approximating continuous transitions...).

The recursion process can then be summarized by the sequences s₁ ∶ s_n-1 and μ₁ ∶ μ_n-1 which depend on C, B, and θ. At the last step (integration over r_n), we integrate the product of the two remaining Gaussian distributions: the previous term $fsn−1(rn−μn−1)$ and the localization error term $fσ(cn−rn)$ as described for the first step to compute the density function $fC,BC,B | θ$⁠.

Finally, we compute the value of the probability density function of the observed track f_C(C | θ) as the sum of f_C,B(C, B | θ) over all possible B.

Since diffusive motion is independent in each spatial dimension, the principle described above for one dimension can simply be extended to two or three dimensions by multiplication of independent distribution functions. For example, in 2D, the function f_C,B(C, B | θ) is simply replaced by the product f_C,B(C_x, B | θ) f_C,B(C_y, B | θ). In principle, each axis can have a different localization error and different diffusion lengths for each state. This is especially true for localization in the direction of the optical axis compared to the lateral axes. ExTrack, therefore, allows independent localization errors for each axis.

Alternatively, the user can also provide localization error for each peak, for example using the Cramer–Rao lower bound estimate (Ober et al., 2004; Lindén et al., 2017). Since these and other estimators might underestimate the true localization error, peak-wise localization estimates can also be implemented as scaling factors that are then assumed to be linearly related to the true localization error estimated by ExTrack.

This method has a number of operations which initially scales with Nⁿ, where N is the number of states and n is the number of time points. This means the calculation time can become unrealistically long when analyzing long tracks. To alleviate this issue, we developed a window method to allow it to work with longer tracks in a reasonable time scaling with nN^m+1 (m the window length of minimal value 1). For computational reasons, we advise to use a window length of 7 for two-state models, 5 for three-state models, and 3 or 4 for more states.

Here, we briefly motivate and describe the implementation of the window method: During the recurrence process described above, Eqs. 5 and 6, f_s can be regarded as a density probability function of the position r_i+1 knowing the previous observed positions c₁, … , c_i and states from positions 1 to i + 1. We realized that the current localization of a particle is very little affected by its state m steps ago when m >> 1. Thus, for a given track, two sequences of states varying only in their first state should give very similar $fsi$⁠. The values μ_m and s_m of these normal distributions should also be similar.

For example, if the track has been diffusive during at least one of the positions from steps i−m to i, the current observed position c_i is much more informative for the real position r_i+1 than the first observed position c₁. If the molecule has been immobile from position i−m to i, all observed positions c_i−m to c_i are equally informative. However, even the past 5–7 positions are likely sufficient to predict the distribution $fsi$⁠.

As we saw in Parameter fitting..., for a given sequence of states B, computing f_C,B(C, B | θ) is nothing but computing three sequences s₁:s_n, μ₁:μ_n, and K₁:K_n until the last step where we simply have to compute $fsn$ (c_n).

In this way, two sequences of states are merged (for example, the sequences starting with [0,0,0,1,1,1] and [1,0,0,1,1,1]). We thus reduce the number of μ_m, $sm2,$ and K_m from 2^m+1 to 2^m. By recursion of this principle over all steps from m to n−1 we limit the computation time to 2^m+1 for a two-state model, or, more generally, to N^m+1 for a N-state model. In the following section, we will introduce sub-steps between the discrete observation time points. Our approach is easily generalized to sub-steps by considering state vectors. Applying our approach with a window length m = 5–7, we observed similar functional dependencies of the likelihood on the model parameters θ, allowing us to drastically speed up our method without losing accuracy.

ExTrack fits data of a continuous-time process to a discrete-time Markov model. Without the introduction of sub-steps, ExTrack assumes that transitions can only happen once per time step. It then estimates transition probabilities per time step, which must be translated into transition rates k_i, that describe the continuous-time Markov model. For continuous-time Markov processes, transition probabilities can be converted into rates according to a simple relationship P = e^GΔt, where P is the transition probability matrix, which contains the transition probabilities p_i,j from state i to j, and where G is the generator matrix with elements G_i,j = k_i,j for i ≠ j and $Gi,i=−∑i≠jki,j$ (Pishro-Nik, 2016). Here, the transition probabilities $pi,j$ allow the molecule to transition from state i to j through any number of intermediate states.

However, the implementation of this relation into ExTrack leads to a systematic overestimation of transition rates. The reason for this overestimation is found in our approximate representation of the distribution of physical displacements between time points (Eq. 4), which is based on the false assumption that transitions can only occur at the middle of steps, contrary to the continuous-time nature of the underlying physical process. We found that this error could be compensated for by using a slightly different approximation for state transitions $pi,j=1−exp−ki,j Δt$⁠. In the limit of small k_i, Δt, this approximation asymptotically equals the exact expression P = e^GΔt, which also asymptotically equals p_i,j = k_i,j Δt. We found that ExTrack using the approximate relationship performs better in the case of two-state and three-state models for a large range of transition rates. However, ExTrack (Python version) also allows using the generator-matrix-based relationship if the user desires.

When transition rates are high (when k Δt > 0.4 for the two-state model), our method allows to sub-divide time steps into a number of u sub-steps (where u = 2 corresponds to dividing each step into two). This allows ExTrack to account for multiple transitions and transition times that are different from the mid-points of time steps Δt. To take states at sub-steps into account, we introduce a new state vector B = (b_1,1, b_1,2, … b_1,u, … , b_n-1,u, b_n,1) and new physical positions, that require integration according to Eq. 1. This integration is straightforward: The probability β (Eq. 2) is simply replaced by the product of all state transitions between subsequent sub-positions. The probability distribution of real positions (Eq. 4) is replaced by the corresponding distribution of sub-positions. Since the physical displacement during Δt is the sum of Gaussian random variables (the sub-displacements), the functional form of Eq. 4 (a product of n−1 Gaussians) can be maintained while replacing the SDs $δi=di2+di+12$ by $δi=∑k=1uδi,k2/u,$ where $δi,k$ are the corresponding diffusion lengths for the sub-steps.

The use of a window length m will allow the user to do accurate computations for m sub-steps. Thus, at a given m, the number of observed positions c considered within the window equals floor(m,u). To consider the same number of observed positions per window, one thus needs to increase the window length. A trade-off between number of states, window length, and number of sub-steps has to be found (see Implementation and computation time).

Tracks can terminate due to different reasons: photobleaching, diffusive molecules leaving the field of view, or molecules transiently not being detected. The process of leaving the field of view requires diffusive motion. Observation of long-lived molecules within a finite field of view can thus show a bias toward non-moving or slowly moving molecules. An extension of ExTrack can take this bias into account by explicitly modeling the probability of track termination. We consider two contributions to track termination: first, a constant termination probability p_K, which is independently of the motion state. This probability summarizes photobleaching and the probability to not detect a molecule, for example because of low signal-to-noise ratio; second, a probability of leaving the field of view (or observation volume) p_L that depends on the diffusion length and the dimensions of the field of view. In the case of a cytoplasmic particle tracked through epi-fluorescence or confocal microscopy, the monitored length is the depth of field (or focal depth). In case of a membrane protein moving around a cylindrical cell imaged in TIRF microscopy, the monitored length is a fraction of the cell diameter (Fig. 2, c and d).

The annotation module can also take advantage of the window approximation described in Using a window... to reduce computation time and make the computation tractable in case of long tracks. Since the annotation module does not require parameter fitting and thus many iterations, the window length can be chosen larger than for the fitting module. A large window length is also more important for precise state prediction than for accurate global parameter fitting.

ExTrack can improve the estimation of molecule positions based on a track, in particular if molecules move slowly. Positions can be estimated by computing the probability density function of each real positions r_i (Lindén et al., 2017). To do so, we compute f_C(C | θ) without integrating at position r_i. This results in a probability density function f (r_i|C, θ) (Fig. 3, d and e), which is a sum of Gaussian functions for each sequence of states. While this probability density can be obtained explicitly, it is much faster to obtain the expected value and SD of the density function. Those values are computed by averaging the parameters of the Gaussian distributions associated with each sequence of states weighted by their respective probability. Like for the fitting method, the window method is applied (see Using a window...).

Computing state-duration histograms allows to assess non-Markovian transition behaviors or to reveal multiple hidden states with different transition rates. For a given state, the resulting rate is then the sum of track-termination rate (bleaching, track termination due to low signal-to-noise ratio, leaving the focal plan) and the transition rates to other states. If the track-termination rate is low, the histogram allows to identify one or multiple transition rates (see, for example, Fig. 3, f and g). Picking only long tracks can help removing the contribution from bleaching.

ExTrack is available as a Python (Python 3) package (Simon, 2022). A version with the core functionality is also available as a TrackMate module (Ershov et al., 2022) on Fiji. The TrackMate implementation can fit data to a two-state model and annotate states according to the results from the fit or manually chosen parameters. It then allows interactive visualization of tracks colored with state probabilities for each displacement. We allow parallelization with GPU (cupy library; Python version only) or multiple CPUs (both Python and TrackMate versions).

As mentioned in Approximating continuous transitions..., a trade-off between number of states N, window length m, and the number of sub-steps u has to be found for reasonable computation time. When running the ExTrack fitting module on a computer with Intel Core i7-9700 processor (10,000 tracks of 10 positions) for 200 iterations using a two-state model, a window length of m = 2, and no sub-step (u = 1) the analysis can be as fast as 20 s. For the dependency of computation time on numbers of states, sub-steps, and window length, see Fig. S1 c.

To save computation time, we recommend to initially run ExTrack with low values of u and m and then to increase u if model predictions suggest high transition rates or m for low predicted diffusion lengths. Specifically, we suggest making the following adjustments: If localization error is negligible, for instance if there is no immobile state and all d_i > 2σ, window length m can be set to its minimal value of 1. Similarly, m = 1 should perform well if there is one immobile state and all diffusive states have large d_i > 5. In such cases, multiple sub-steps can be used at little computation cost. More generally, if predicted transition rates are larger than 0.4 Δt⁻¹ but localization error is not negligible, we suggest increasing u to 2 for most accurate estimates (Fig. S1 b vs. Fig. 1 f). In the hardest cases of small d^∗ ≲ 2 and high transition rates (≥0.4), we recommend using u = 2 and m ≥ 8.

To test the predictive power of the different methods, we conducted overdamped Brownian Dynamics simulations of tracks in two or three spatial dimensions with molecules transitioning randomly between the different states at discrete time points. To mimic a continuous-time Markov model for state transitions, we used a small time step τ = Δt/50 ≪ 1/k_i,j, where k_i,j are the transition rates. Brownian Dynamics simulations were carried out by randomly drawing physical displacements in each spatial dimension from Gaussian distributions of SD $2Dτ,$ where D is the diffusion coefficient corresponding to the diffusive state. An additional Gaussian distributed noise of SD σ was added to simulate localization uncertainty.

Except if specified otherwise, we simulated 10,000 tracks of 10 time points with localization error σ = 0.02 µm, Δt = 0.06 s, d₀ = 0 µm, d₁ = 0.1 µm, bound fraction F₀ = 0.5, transition fractions per steps k₀₁ = k₁₀ = 0.1 Δt⁻¹, with infinite field of view and perfectly stroboscopic tracks. We also assumed that molecules reached steady state, that is, $F0⋅k01=(1−F0)⋅k10.$

To test the robustness of ExTrack and the other methods to more complex behaviors, we also simulated tracks with variations of localization error or diffusion coefficients, a finite field of view or physical confinement as follows (see Fig. 2 for illustrations):

Track-to-track variations of localization error (or diffusion coefficients) were simulated with localization error σ (or diffusion coefficient D) following χ² distributions of given coefficients of variation and mean 0.02 µm (or 0.25 µm² s⁻¹ for D), Δt = 0.02 s. Models with multiple diffusion states were simulated as continuous-time transitions with model parameters detailed in Table S2.

To simulate a finite field of view in two dimensions, we simulated tracks in a box that is infinite in one spatial dimension (y) and finite in the other dimension (x) with size 3l, where l is the size of the field of view. All tracks or part of tracks that fall into the field of view are considered for further analysis. A single particle can thus result in several tracks if leaving the field of view and coming back. A finite field of view in three spatial dimensions was simulated analogously: The simulation box is infinite in x- and y-directions, while the box has periodic boundary conditions in the z-direction.

To simulate physical confinement, we considered tracks to move within a square area of indicated side length, using reflecting boundary conditions.

We also simulated tracks with one bound state and several intractable diffusive states to determine whether ExTrack could accurately estimate global binding and unbind with a coarse-grained three-state model. In Fig. S3, we show the results obtained from simulations of tracks with one bound state and five diffusive states of d^∗ = [2,3,4,5,6]. transition rates (in Δt⁻¹): k₀₁ = 0.03, k₀₂ = 0.02, k₀₃ = 0.06, k₀₄ = 0.04, k₀₅ = 0.05, k₁₀ = 0.04, k₁₂ = 0.05, k₁₃ = 0.1, k₁₄ = 0.07, k₁₅ = 0.04, k₂₀ = 0.19, k₂₁ = 0.01, k₂₃ = 0.04, k₂₄ = 0.06, k₂₅ = 0.1, k₃₀ = 0.02, k₃₁ = 0.03, k₃₂ = 0.05, k₃₄ = 0.05, k₃₅ = 0.1, k₄₀ = 0.1, k₄₁ = 0.04, k₄₂ = 0.04, k₄₃ = 0.01, k₄₅ = 0.05, k₅₀ = 0.15, k₅₁ = 0.05, k₅₂ = 0.01, k₅₃ = 0.02, k₅₄ = 0.05. The 10,000 tracks of 10 time points were then fitted to a three-state model, window length = 6.

To compare our results with vbSPT, we fixed the number of states to two so both algorithms performed exactly the same task. vbSPT does not consider localization error but a metric that we will call u. In case of pure diffusion, u = D ⋅ Δt but in case of immobile particle with localization error u = σ²/2 in principle. We can thus infer σ and D according to $σ=2⋅u0$ and D = (u₁ − u₀)/Δ.

To assess the propensity of RodZ molecules to rebind in close vicinity of their initial binding site, we first annotated tracks using parameters obtained from the ExTrack fitting module (three biological replicates). We considered the 16 first time points of tracks of at least 16 time points (pooling the tracks of the three replicates). Among tracks labeled as initially immobile for at least three time points (p_immobile > 0.5) then diffusive for four time points (with at least three time points of probability p_diffusive > 0.7), we grouped tracks into two sub-groups, the ones rebinding right after and the ones which continue to diffuse for at least one more time point. The histograms represent the distributions of distances between initially bound position and the position at the fourth time point after unbinding. If molecules were to rebind at random locations, the distributions of the distribution of distances for rebinding particles should be the as for particles, which continue to diffuse. Neither PBP1b (Fig. S5 e) nor tracks obtained from immobile fluorescent beads (of similar signal-to-noise ratio) showed any significant rebinding, which excludes wrong conclusions on RodZ data due to mis-annotations.

We used the IPTG-inducible GFP-RodZ strain FB60(iFB273; ΔrodZ, Plac:: gfp-rodZ) by Bendezú et al. (2009) and the GFP-PBP1b-containing strain AV51 (msfgfp-mrcB, ΔmrcA; Vigouroux et al., 2020). Cells were grown overnight at 37°C (shaking) in LB medium and then washed and diluted at least 1:1,000 in M63 minimal medium (Miller, 1,972) supplemented with 0.1% casamino acids, thiamine (5 × 10⁻⁵%), glucose (0.2%), and MgSO4 (1 mM) and grown for 6 h to early exponential phase (maximum OD600 of 0.1) at 30°C (shaking). Cells were then spread on an agar pad made from the same M63 media as described above. RodZ production was induced with 100 µM IPTG. In the strain AV51, CRISPR repression of msfGFP-PBP1b is induced with 100 ng/ml of anhydro-tetracycline (Acros Organics). When necessary, strains were supplemented with kanamycine (50 µg/ml) or carbenicillin (100 µg/ml) during overnight cultures. Biological replicates result from independent cultures starting from separate colonies.

Cells were all positioned in the same focal plan in between an agar pad (1%) and a coverslip to be imaged in TIRF microscopy. Coverslips were cleaned by 60 min sonication in saturated KOH solution followed by two washing steps (15 min sonication in milli-Q water). Single-particle tracking of GFP-PBP1b was performed with a custom-designed fluorescence microscope based on an ASI Rapid Automated Modular Microscope System, equipped with a 100× TIRF objective (Apo TIRF, 100×, NA 1.49, Nikon), Coherent Sapphire 488–200 laser, and a dichroic beamsplitter (Di03-R488/561-t3-25 × 36, Semrock). Excitation was controlled with an acousto-optic tunable filter (AA Optoelectronics) through an Arduino (15 ms light exposure per frame). Images were acquired using an Andor iXon Ultra EMCCD camera with an effective pixel size of 130 nm. Image acquisition was supervised with MicroManager.

Tracks were built from movies using TrackMate with LoG peak detection (estimated blob diameter of 0.5 µm, quality threshold of 15, with subpixel localization). Then, peaks were linked into tracks using the “Simple LAP Tracker” (max distances of 0.4 µm, and “gap” of 2).

Data analysis with ExTrack was restricted to tracks with at least three positions. For the longest tracks, only the first 50 positions were analyzed.

Fig. S1 shows ExTrack parameter fitting. Fig. S2 shows error on two-state model parameters for different methods. Fig. S3 shows robustness of ExTrack to biases due to wrong model assumptions. Fig. S4 shows capacity of the annotation module. Fig. S5 shows complementary results for GFP-PBP1b and GFP-RodZ tracks. Video 1 shows msfGFP-PBP1b single-particle tracking. Video 2 shows GFP-RodZ single-particle tracking. Table S1 shows using ExTrack to fit parameters of a three-state model to different simulations of three-state data with qualitatively different types of transitions. Table S2 shows two-state and three-state fits of tracks from simulated particles either in immobile state (state 0) or in one of five diffusive states (states 1–5).

Single-particle tracks, ExTrack results, and selected movies of the experiments corresponding to Fig. 4 are available on the Zenondo repository (https://doi.org/10.5281/zenodo.7548793).

We thank the former lab members Andrey Aristov for setting up the TIRF microscope, Antoine Vigouroux for the GFP-PBP1b strains, and Gizem Özbaykal for her guidance for experiments. We also thank Felipe Bendezú (Umeå University, Sweden) and Piet De Boer (Case Western Reserve University, Cleveland, OH, USA) for the GFP-RodZ fusion.

This work was supported by the European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant Agreement No. 679980) to S. van Teeffelen, the French Government’s Investissement d’Avenir program Laboratoire d’Excellence “Integrative Biology of Emerging Infectious Diseases” (ANR-10-LABX-62-IBEID, S. van Teeffelen), France BioImaging (ANR-10-INBS-04, J.-Y. Tinevez), the Mairie de Paris “Emergence(s)” program to S. van Teeffelen, a Natural Sciences and Engineering Research Council of Canada Discovery Grant to S. van Teeffelen, a Fonds de Recherche du Québec - Santé Salary Fellowship to S. van Teeffelen, as well as support from the Volkswagen Foundation to S. van Teeffelen. Open Access funding provided by Université de Montréal.

Author contributions: F. Simon: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Writing. He conceived and developed ExTrack, carried out all experiments, conducted the analysis, and wrote the paper. J.-Y. Tinevez: Software. He ported ExTrack to the Fiji plugin TrackMate. S. van Teeffelen: Conceptualization, Methodology, Writing, Supervision, Project administration, Funding acquisition. All authors reviewed the final manuscript.