High-throughput super-resolution single-particle trajectory analysis reconstructs organelle dynamics and membrane reorganization.

Super-resolution imaging can generate thousands of single-particle trajectories. These data can potentially reconstruct subcellular organization and dynamics, as well as measure disease-linked changes. However, computational methods that can derive quantitative information from such massive datasets are currently lacking. We present data analysis and algorithms that are broadly applicable to reveal local binding and trafficking interactions and organization of dynamic subcellular sites. We applied this analysis to the endoplasmic reticulum and neuronal membrane. The method is based on spatiotemporal segmentation that explores data at multiple levels and detects the architecture and boundaries of high-density regions in areas measuring hundreds of nanometers. By connecting dense regions, we reconstructed the network topology of the endoplasmic reticulum (ER), as well as molecular flow redistribution and the local space explored by trajectories. The presented methods are available as an ImageJ plugin that can be applied to large datasets of overlapping trajectories offering a standard of single-particle trajectory (SPT) metrics.

Introduction

Subcellular compartments are focused sites where large numbers of molecules dynamically interact to support cellular function (Cole et al., 1996). The trajectories of ions and proteins as they move between the cytoplasm, plasma membrane, and organelles are critical to cellular function (Betzig et al., 2006; Cohen et al., 2018). These dynamics occur at the endoplasmic reticulum (ER) (Wu et al., 2017, 2018; Petkovic et al., 2014), the mitochondrial network, endosomes and lysosomes, and microtubules and are impacted by the local properties of these different environments (Lu et al., 2020). Several experimental paradigms measure these constitutive molecular motions at subcellular sites, including fluorescence recovery after photobleaching (FRAP) (Axelrod et al., 1976; Saxton and Jacobson, 1997), which locally depletes fluorescence and measures the timescale and fraction of recovery. Analysis of these data can reveal trafficking at a population level. In contrast, photoactivation (Wang et al., 2014) consists of activating molecules in a local region of the cell and reveals their spread over a transient time frame. Combined with diffusion modeling and stochastic simulations, various biophysical properties can be measured, including diffusion coefficients and the fraction and timescale of recovery (Lippincott-Schwartz et al., 2003). These methods provide information on the dynamic function of organelles but are insufficient to identify and reconstruct high-density regions. These approaches also cannot examine phase separation stability or the local spaces explored by molecules at a nanoscale resolution. Statistical analysis (Manley et al., 2008; Hoze et al., 2012) of a large ensemble of super-resolution single-particle trajectories (SPTs) (Figures 1A and 1B) has the potential to reveal local molecular interactions. Molecules are not uniformly distributed inside a cell but instead form heterogeneous aggregates, possibly in phase-separated nanodomains, characterized by high-density regions (HDRs). Such regions are characterized by reduced velocity movement of molecules and confinement. These local areas can be enriched with calcium channels at neuronal synapses and store-operated calcium entry receptors such as STIM1 on spine apparatus and can also be found at ER nodes (Heck et al., 2019; Wu et al., 2014; Heine and Holcman, 2020; Holcman et al., 2018). Interestingly, these ubiquitous structures are transient yet persist, with a timescale that is longer than that associated with molecular trafficking. In short, many interactions that are critical to the cell are regulated and controlled by subcellular mechanisms that currently cannot be easily captured by quantitative analysis.

Figure 1

High-throughput SPT analysis pipeline

For a Figure360 author presentation of this figure, see https://doi.org/10.1016/j.crmeth.2022.100277.

High-throughput automated trajectory analysis workflow.

(A) Acquisition device and raw data of a single-particle experiment.

(B) Raw data from (A) are transformed into trajectories using classical software such as Trackmate, available as an ImageJ plugin.

(C) Schematic description of the high-throughput analysis implemented here: trajectories are first discretized both spatially, using a square grid, and temporally, using temporal binning (time-windows analysis). We then interpret the trajectories based on the Langevin equation, allowing us to generate high-resolution maps of the local trajectories motion. High-density/low-velocity regions of the maps are extracted by automated algorithms to detect potential wells/reconstruct the network. Finally, the outputs consist in statitstics associated with well locations and with reconstructed network. These characteristics allow to analyze how trajectories locally explore their environment at the nanometer scale.

To determine the underlying physical properties of molecular trafficking, various computational modeling methods have been developed to analyze SPTs. These models include those based on classical free (Crank, 1979) and confined (Kusumi et al., 1993; Saxton, 1995; Saxton and Jacobson, 1997; Holcman and Schuss, 2004) diffusion, active deterministic motion, or a mixture of deterministic and stochastic models (Hozé and Holcman, 2017). Based on this theoretical framework, analysis of SPTs has revealed the dynamics of local chromatin organization in the nucleus (Gasser, 2016; Amitai and Holcman, 2017; Hauer et al., 2017), synaptic receptor trafficking at neuronal synapses (Dupuis and Groc, 2020), and vesicular stomatitis virus G protein (VSVG) virus assembly (Manley et al., 2008). A significant recent advance is the analysis of massive numbers of overlapping SPTs. This statistical analyses can reveal the properties of molecular trajectories. However, analyzing these data can potentially also allow quantification of membrane dynamics (Holcman et al., 2015) and may give insight into organelle organization.

Software developed to analyze SPTs generally fall into two categories: (1) those where parameters are extracted along individual trajectories, and (2) those based on spatially combining these trajectories to recover intrinsic properties. The first category includes SpotOn (Hansen et al., 2018), used to fit multi-states diffusion models from the distribution of displacements, the 4P-parameter algorithm, used to reconstruct chromatin dynamic (Amitai et al., 2017; Shukron et al., 2019), and, more recently, the 4P-Gaussian-Mixture (Basu et al., 2021), used to subsegment trajectories into confined and unconfined regions. For the second category, common methods include diffusion/drift maps and potential well extractions pioneered in Hoze et al. (2012) (https://bionewmetrics.org/super-resolution-single-particle-trajectories-using-stochastic-analysis/) or the SR-Tesseler algorithm (Levet et al., 2015), used to reconstruct and extract biophysical parameters from the local density of points. Another package is TRamWAy (Beheiry et al., 2015), used to reconstruct diffusion and energy maps, but this approach does not reconstruct the potential well boundaries.

Here, we have developed a method based on hybrid algorithms and automated an analysis pipeline for SPTs (Figure 1C) based on the Langevin (Schuss, 2009) equation of motion to provide statistical analysis of these data. This method estimates biophysical characteristics and is capable of reconstructing nanodomain sizes and boundaries using the classical physical model of a potential well (Kramers, 1940; Schuss, 2009), well known since Kramer’s work in 1940 (Kramers, 1940). This method allowed us to characterize calcium channel organization in the membranes of hippocampal neurons. We also present an algorithm that reconstructs a network from SPTs based on the clustering of low-velocity trajectory fragments and use it to reconstruct ER network topology as well as estimate the timescale of lysosome trafficking and ER network interactions. Finally, these methods allowed us to extract the motion of trajectories relative to the reconstructed network, thereby revealing the redistribution of trajectories inside the ER of normal and atlastin-null mutant cells (Niu et al., 2019). The diversity of datasets used here demonstrate the broad applicability of these methods to cell biological processes. The algorithms we developed here are available in an ImageJ plugin with a graphical interface for processing individual experiments and a programming interface for batch processing.

Results

Cosine filtering improves diffusion map accuracy

The diffusion coefficient is extracted from SPTs locally and is used to construct the diffusion map. Such maps have been previously constructed from a square grid overlaid on top of the trajectories (Hoze et al., 2012; Heck et al., 2019), but this estimation is strongly affected by the number of displacements falling in each bin. To improve the accuracy of the diffusion map, we propose to use a cosine filter that computes for each bin of a square grid the diffusion coefficient based on all the displacements starting within some radius of the bin center (larger than the bin size). Each displacement is weighted by a cosine function depending on its distance to the center of the bin (see Method details).

We tested the capacity of the method to reconstruct a synthetic diffusion field, either with a constant or with local variations using numerical simulations (see Methods S1). We find that the cosine-filter method has a smaller overall error and much greater coverage than the classical diffusion map (Figure S1). We then compared the results obtained by raw diffusion map, cosine filtering, and local averaging on empirical data (Figure S2; Data S1) and found that they all give similar results, but cosine filtering recovers more local details compared with the other methods.

Algorithm to reconstruct nanodomains in HDRs

In order to extract subcellular regions where a high density of molecular interactions are occurring, we developed a novel computational approach and automated pipeline (Figure 1) based on stochastic equations, multi-scale analysis, optimal estimators, and maximum likelihood statistics. In contrast to methods used to extract density of points (Sibarita, 2014; Levet et al., 2015) or to compute the maximum likelihood estimator (MLE) (Parutto et al., 2019; Briane et al., 2020), the present approach combines density of points with local dynamics associated with the elementary displacement , where is the position of the particle at time (Hoze et al., 2012; Heck et al., 2019) to extract the field of forces. Indeed, trajectories following the stochastic Equation (2), where the forces defining the nanodomain given by Equation (12) are characterized by a local accumulation approximated by a Gaussian density of point, and converging arrows for the local drift field (Figures S3Aa–S3Ac; Data S1). We could thereby assess the nature, organization, and stability of a large amount of HDRs by collecting statistics that can reveal hidden cellular organization. HDRs could previously only be characterized by manual curation based on extracting parameters of potential wells, a concept in classical physics (Chandrasekhar, 1943; Schuss, 2009) that describes the stability of dynamic systems, such as the motion of a bead attached to a spring. In contrast, using combined optimization procedures, the present method allows users to automatically extract local diffusion coefficients, energetics, local field potential, and, most importantly, local boundaries.

The method relies on an interesting observation that nanodomains of high density tend to have an elliptic shape. Thus, the first step to detect them is to recover their center and boundaries (see Method details). Non-automated classification algorithms have used the density of points (Parutto et al., 2019) or displacements separately, but these procedures often lead to parameter estimations that are not completely satisfactory due to a shallow minimum of the associated error function that leads to a large variability and possible mistakes in the estimation of most of the nanodomain parameters such as the boundary and the energy of nanodomains.

To overcome these difficulties, we developed a hybrid algorithm (Figure 1), described in the Method details, which combines two independent procedures starting with a principal-component analysis to recover the elliptic boundary and followed by a maximum likelihood estimation of the effective diffusion coefficient and drift properties. More precisely, the new algorithm comprises three steps:

Automatic determination of bins with the highest density of points (Figure S3Ba).

For each such bin, we iterate over square regions of increasing width around the bin center. For each iteration , we apply a principal-component analysis to estimate the semi-axis , the center , and a score (Equation 24) based on the points falling in the square of size (Figures S3Bb and S3Bc). We iterate until we reach the maximum size of the well set by the user (Figure S3Bc).

In the termination step, we select the optimal value of the iteration providing the optimal parameter (see Method details).

The reconstruction is illustrated in Figure S3C for three wells. To evaluate the performances of this hybrid algorithm, we constructed ground-truth datasets consisting of trajectories following Equation (2), with the energy field given by 12 for different types of potential wells (see Methods S1 for the detail of the procedures). In order to be able to compare different algorithms in a fair manner, we developed a parameter optimization procedure where a grid search is used to find the parameters of each algorithm that gives most accurate reconstruction of a known potential well. Based on this procedure, we showed that the hybrid algorithm performs better than two previous algorithms, the drift-based algorithm from Heck et al. (2019) and the density-based algorithm from Parutto et al. (2019), in estimating all parameters, especially at higher well energies (see Tables S1 and S2). The algorithm is effective at reconstructing wells over a large range of energies, from 1 to 10 kT (Table S3; Figure S3D), and can effectively discriminate the presence or absence of a parabolic field (see also Methods S1 for detection on wells with kT). Moreover, for sufficiently large regions, the algorithm can distinguish between a parabolic field and a Brownian motion confined by impenetrable walls (see Methods S1 and Data S1).

At this stage, we have thus validated the hybrid algorithm by using ground-truth datasets to identify and automatically reconstruct nanodomains. This hybrid algorithm outperformed the others when reconstructing the exact shape of the nanodomain and its energy.

Analysis of SPTs for the endogeneous voltage-gated calcium channels reveals their organization and weaker stability in nanodomains

We applied the new hybrid algorithm to reanalyze trajectories of calcium voltage channels (CaV2.1) on the surface of neuronal cells for two overexpressed splice variants, CaV2.1 and CaV2.1 , previously shown to shape synaptic short-term plasticity (Heck et al., 2019), as well a new set of endogenously tagged CaV2.1 channels. Using a large number of redundant SPTs, we were able to automatically detect the nanodomains defined as HDRs. Representative examples of such trajectories and their associated density and diffusion maps are presented in Figures 2A–2C. The algorithm allowed us to identify the characteristics of the nanodomains (reported in Table S4) approximated as ellipses with semi-axis lengths and nm for Cav2.1 . These parameters are similar for Cav2.1 , while they were smaller for endogenous Cav2.1 with and nm (Figure 2D). The diffusion coefficient was for the variant and for the variant and was smaller for endogenous CaV2.1, with (Figure 2E). Interestingly although the associated energy was for the two variants and for endogenous CaV2.1 (Figures 2F and 2G), we found that the associated residence times of a receptor in a nanodomain was around 174–181 and 94 ms, respectively (Figure 2H). The estimations of these residence times are larger than the ones reported in Heck et al. (2019), which could be due in part to differences in the boundary estimated by the present algorithm and in part due to the different estimators used for computing the attraction and diffusion coefficients of potential wells. Indeed, the estimators constructed from the statistical moments used in Heck et al. (2019) have a tendency to underestimate these parameters, whereas the MLEs used here give values closer to the ground truth (see Tables S1–S3 and Methods S1). Interestingly, we found a positive correlation between the diffusion coefficients and the sizes as well as the attraction coefficients but no correlation with the energy of the potential wells (Figures S6A–S6C; Methods S1). Finally, we conclude that the present approach captures the differences between the variant forms of the calcium channels, revealing a significant reduction in interactions for the endogenous dataset (see Methods S1).

Figure 2

Automated large statistics of CaV nanodomains located on dendrites

Super-resolution SPTs automated analysis reveals CaV nanodomain organization.

(A) Examples of endogenous CaV2.1 trajectories together with the detected potential well regions (red ellipses).(B) Associated density map, presenting the local point density (in (points/)) and computed over a grid with bin sizes of 50 nm and locally averaged over a 3 × 3 Gaussian kernel.

(C) Associated cosine-filtered diffusion maps (see cosine filtering improves diffusion map accuracy), displaying the local diffusion coefficients (in ) computed over a grid with bin sizes nm, for bins possessing at least 15 displacements and using a disk radius nm.

(D–H) Violin plots representing the distribution, median, and interquartile range of the well characteristics recovered using the hybrid algorithm from CaV2.1 47, CaV2.1 +47, or endogenous CaV2.1 trajectories. (D) Distribution of semi-axes lengths of elliptic well boundaries, (E) diffusion coefficients inside wells, (F) attraction coefficients , (G) energies of the well , and (H) residence time distribution of a trajectory inside a well.

See also Tables S4 and S5.

Time-lapse analysis reveals the stability of CaV nanodomains over time

To investigate the stability of nanodomains across time, we used a time-lapse analysis (Figure S4; Data S1) with sliding windows of 20 s and no overlap. This analysis allows determination of the lifetime of a nanodomain, which is given by the number of successive windows where it is detected (Figure S5A; Data S1). For example, the trajectories obtained during a 250 s experiment are split into 13 20 s windows (0–20, 20–40, …, 240–260 s). We searched for the presence of potential wells in each window (Figures S5A and S5B). To follow a well across successive time windows, we consider that two wells identified at times and are the same if the distance between their centers is less than 200 nm. The ensemble of consecutive times , where a well is first detected at time and disappears at time , is used to define the stability duration . This analysis allows us to follow the size of the small and large axes of the wells and the associated energy over time (Figures S5C and S5D). Finally, we found that ∼55% of the wells were present for more than 20 s (one time window) and that their average duration is ∼56 s for both and variants (Figure S5E; Table S5). However, nanodomain stability is reduced to 46 s for endogenous CaV2.1. All of these durations are longer than the ∼30 s that we previously reported (Heck et al., 2019), indicating that the present algorithm advances our ability to capture the dynamics of these high-density enigmatic subcellular domains.

The hybrid algorithm reveals that some ER nodes are nanodomains defined by an attracting field of force

To further explore the range of applicability of the present hybrid algorithm, we analyzed SPTs recorded from an ER luminal probe in COS-7 wild-type (WT) and HEK-293T cells, both presented in Holcman et al. (2018) and in COS-7 atlastin knockout (dATL) cells obtained as in Holcman et al. (2018) (see Methods S1). The ER networks of these cells generally form a flat tubular monolayer at their periphery (Schroeder et al., 2019), allowing the two-dimensional tracking of individual molecules as shown in Holcman et al. (2018). In the dATL mutant, the morphology of peripheral ER tubules is altered, but it is unclear how the ER flow is affected. Since nodes have been previously characterized as HDRs (Holcman et al., 2018), we asked here whether these nanodomains could further be defined as potential wells. Applying the hybrid algorithm reveals several potential wells (Figures 3A and 3B, red ellipses) precisely located at nodes forming HDRs. We further estimated the density, diffusion and drift maps, observing converging arrows patterns in these regions (Figures 3D and 3E, and S6D), a classical feature of potential wells (Hozé and Holcman, 2017). Interestingly, these HDRs were characterized by ellipses with large semi-axis lengths ∼220 nm, , small semi-axis lengths ∼160 nm, diffusion coefficients , and attraction coefficients , corresponding to energies and residence times ms (Figures 3F–3J; Table S6). Interestingly, although the elliptic parameters are not much different in the case of COS-7 WT and COS-7 dATL, a difference can be observed in the dynamical parameters characterizing the transport of the material across the ER network. To conclude, the present hybrid algorithm reveals that some ER nodes concentrate trafficking of luminal molecules by a spring-force type mechanism, the origin of which should be further explored.

Figure 3

HDR present in organelle such as the ER network characterized as potential wells

(A) ER network geometry (top), associated luminal trajectories (middle), and overlay (bottom).

(B) Individually color-coded trajectories recorded inside the ER for three different cells (HEK293t, COS-7, COS7-dATL). Red ellipses represent detected potential wells.

(C–E) Density, diffusion, and drift maps for the regions shown in (A). Arrows in the drift maps are colored according to the direction: West (purple), East (green), South (blue), and North (red).

(F–J) Violin plots representing the distribution, median, and interquartile range of the characteristics of the detected potential wells. (F) Semi-axis lengths (large and small ), (G) diffusion coefficient , (H) attraction coefficient , (I) energy of the well , and (J) estimated residence time.

See also Table S6.

Organelle network reconstruction from a large number of SPTs

We next wanted to check whether our approach can be used to define the structure of the ER and be generalized to other organelles. For these goals, we developed a novel method and algorithm to reconstruct the network from SPTs.

Graph reconstruction algorithm (GRA) to unravel the ER network

Although SPTs can be used to explore the ER network architecture (Holcman et al., 2018), we still lacked a method for the automated reconstruction of the ER, especially in cases where the local density is variable. We therefore developed an algorithm to see if the dynamic architecture of this complex organelle can be recovered from SPTs. We illustrate the principle of the automatic ER-reconstruction procedure in Figure S7 and (Methods S1) for ER-luminal trajectories. The starting point is to color trajectory displacements depending on their instantaneous velocity, which reveals a dynamical segregation of the ER into nodes and tubules (Figures S7A–S7C). Based on this segregation, we developed the GRA to recover the ER structure from SPTs (see Methods S1).

This procedure allows us to reconstruct a two-dimensional graph for the organelle network that can be used to study further statistical properties. Finally, we tested the GRA on a ground-truth dataset exacted from a live-cell image (Figure S8; Methods S1). We simulated trajectories on the extracted graph and applied our algorithm on these trajectories to reconstruct the topology of the underlying ER network. This procedure gave a satisfactory reconstruction result (8.1% error), thus validating the present algorithm (Table S7).

GRA of the ER network of COS-7 dATL cells reveals aberrant organization and trafficking

Next, we examined whether disruptions to organelle structure and dynamics can be captured using our algorithm. We analyzed SPTs recorded in the ER of COS-7 dATL cells (acquired as in Holcman et al. [2018]; see Methods S1), which lack the ER membrane-shaping protein atlastin (double knockout of the ATL2/3 genes [Niu et al., 2019]) and exhibit a disrupted peripheral ER morphology with elongated tubules (Niu et al., 2019). We present color-coded trajectories based on their instantaneous velocity (Figures 4A and 4B) and the density and diffusion maps (Figures 4C and 4D), as well as the histogram of the apparent diffusion coefficients (Figure 4E), revealing an average .

Figure 4

ER network reconstructed in a COS-7 dATL cell

(A and B) Individual trajectories color coded according to their instantaneous velocity shown in (B).

(C and D) Corresponding density (C) and drift (D) maps.

(E) Distribution of diffusion coefficients obtained from the individual bins of the diffusion map presented in (D), with the average SD (n is the number of bins).

(F) Reconstructed network, showing the nodes in red and links in blue overlaid on top of the individual trajectories (black), with the average SD (n is the number of links).

(G) Distribution of distances (i.e., tubule lengths) between connected nodes, with the average ± SD (n is the number of links) .

(H) Distribution of the areas covered by nodes, with the average SD (n is the number of nodes).

We applied the GRA to obtain a network from these trajectories (Figure 4F) where the HDRs (brown) are connected by blue segments. This approach allows the quantification of the distribution of distances between nodes with a mean (Figure 4G) and the nodes area (Figure 4H). Note that the GRA could miss non-explored ER regions or regions that are underrepresented in the trajectories. To conclude, our algorithm allows us to automatically reconstruct organelle networks based on SPT exploration once there are enough datapoints.

GRA and SPT segmentation reveal the duration of ER-lysosome interactions

To demonstrate the broad applicability of the algorithms, we analyzed an ensemble of trajectories of individual lysosomes obtained in COS-7 cells from Lu et al. (2020) (Figure 5A). These trajectories are characterized by a distribution of instantaneous velocities in the range (Figure 5B). However, low and high velocities are not segregated, as was the case for ER luminal trajectories, but are instead found in similar regions (Figure 5A, left and right). The distribution of velocities (Figure 5B) can be fitted by a sum of two exponentials:where a best fit approximation gives (95% confidence interval ) and (95% confidence interval ) (coefficient of determination ), with and . This fit suggests that the distribution of lysosomes is largely driven by low-velocity components. The rare appearance of high-velocity components suggests a possible switch between slow and fast motions. Finally, note that 86.1% of displacements are associated with a velocity of less than , and 12.7% are in the range of [0.5–1.5] .

Figure 5

Lysosome trajectories analysis

(A and B) Individual lysosome trajectories displacements color coded according to their instantaneous velocity shown in (B). The black line in (B) corresponds to a fit to a two-exponentials function (n is the number of displacements).

(C and D) Corresponding density (C) and diffusion (D) maps.

(E) Histogram of diffusion coefficients obtained from the individual bins presented in (D), with the average SD (n is the number of bins).

(F) Magnification of the density map of two regions of interest presented in (C), showing high-density regions, approximated by ellipses.

(G) Boxplot showing the median, interquartile range, and extreme data points of the distributions of semi-axis lengths of the high-density regions approximated as ellipses.

(H) Reconstruction of a lysosome graph, where nodes (red ellipses) correspond to high-density regions and a link (in yellow) is added when at least one trajectory starts in one node and enter to the other one in one or two frames.

(I) Average instantaneous velocities between pairs of connected nodes presented in (H), with the average SD and a fit to a Rayleigh distribution (n is the number of velocities between nodes).

(J) Percentage of displacements with a specific instantaneous velocity. Inset, percentage of displacements for the velocity regimes defined in (A) and (B).

To further study how lysosomes move in the cytoplasm, we computed relevant density and diffusion maps (Figures 5C and 5D) and found that the motion had a diffusion component (Figure 5E), with an average apparent diffusion coefficient of . Interestingly, regions of low-diffusion coefficients colocalized with regions of high density in the density map (Hoze and Holcman, 2014; Holcman et al., 2015) (Figures 5A–5C).

We then isolated regions of high density using a method based on the density of points (Method details), revealing an ensemble of n = 95 HDR subdomains, approximated by ellipses (black in Figure 5F) of semi-axis lengths (large) and nm (small) (Figure 5G). By considering the displacements connecting different regions, we reconstructed (Method details) a network explored by the lysosomes (Figure 5H), where HDRs (red circles) are connected by direct lines (yellow). Interestingly, the histogram of average velocities between these regions is not symmetric (Figure 5I), with a mean velocity , which clearly deviates from diffusion, as computed from the Rayleigh distribution. This deviation suggests that the transitions between these regions are driven by an active motion. Moreover, the overlay between ER (white) and the lysosome reconstructed network (Figure 5H) suggests that the lysosome trajectories follow the topology of the ER network (Lu et al., 2020). To conclude, our analysis reveals that lysosomes travel along a network that strongly colocalizes with the ER. However, high and low velocities occur in similar regions. Since lysosomes move along microtubules, this present statistical analysis suggests that the ER-microtubule network shapes lysosome trafficking.

Trajectory resynchronization approach reveals single local molecular dynamic exploration inside an ensemble

In the previous result sections, we reconstructed the networks hidden inside SPT data. We shall now introduce a last step in our method, which resynchronizes trajectories that fall inside the same subcellular area but that were acquired at random times. This approach allows us to study the dynamics of trajectories with respect to the ensemble of trajectories that visit the same spatial region. The approach also enables determination of the local spatiotemporal properties that trajectories explore at the single-unit level.

Trajectory segmentation reveals ER-lysosomes interaction timescale

To study the possible interactions between lysosomes and the ER (Figure 6A), we focused on the confined portion found along individual lysosome trajectories (see Method details and Figure 6A). We hypothesize that the lysosome motion can switch between a directed and confined motion (Figure 6B). We first show that the lysosome can indeed switch between directed and confined motion (Figure 6C).

Figure 6

ER network-lysosome interaction revealed by confined trajectories

(A) Schematic representation of a lysosome switching between a directed motion along microtubule and confined motion at ER nodes.

(B) Switching dynamics representation: the confined state is characterized by a spring constant , a diffusion coefficient , and a confined time constant .

(C) Three examples of confinement regions for different lysosome trajectories. Confined trajectories are colored, and unconfined trajectories are in gray.

(D–F) Boxplot showing the median, interquartile range, and extreme data points of the semi-axis lengths (small and large) of the ellipses fitted to the confinement regions (D), spring constants of the Ornstein-Ulenbeck process associated with the confinement regions (E), and diffusion coefficients estimated inside a confinement region (F). n = 818 confinement regions.

(G) Residence times inside a confinement region with a fit to a single exponential (n is the number of confinement regions).

(H) Fraction of time trajectories spend confined (relative to the trajectory length); n = 577 trajectories with at least one confinement event.

(I) Number of confinement events along individual trajectories; n = 521 trajectories with between 1 and 5 confinement events.

To recover the size of the confinement areas, we fitted ellipses over these regions and obtained average semi-axes lengths (Figure 6D) of (large) and nm (small). Furthermore, this approach allowed us to estimate the confinement strength by considering that the confined motion could be generated by a spring force, modeled by an Ornstein-Uhlenbeck process (Schuss, 2009). We found that the spring force is (Figure 6E), associated with an average local diffusion coefficient of (Figure 6F) for a total of n = 818 confinement regions. Finally, the distribution of times in confined regions could be well approximated by a single exponential with a time constant s (Figure 6G). The average residence time of lysosomes in these regions is s, which can be interpreted as a time where lysosomes could interact with the ER. These confinement events are quite common, with trajectories spending ∼30% of their lifetime confined (Figure 6H) in one or multiple confinement events (Figure 6I). To conclude, the present algorithm reveals that as lysosomes travel along a network that strongly colocalizes with the ER, the velocity can switch from large to small displacements, and the trajectories can become restricted into regions of size ∼200 nm, on a timescale of 5 s. This could correspond to a change in directionality of movement or a direct interaction with the ER. Our approach can therefore extract changes in lysosome dynamics that may reflect functional interactions from complex data.

Trajectory resynchronization approach shows how a single trajectory explores single nodes

After an ER network is reconstructed from SPTs by the GRA (Figure 7A), the node-tubule topology emerges. Thus, it becomes possible to study how trajectories locally explore the network by synchronizing them upon exit from a chosen node (Figure 7B). Interestingly, we found that the mean instantaneous velocity at exit is and keeps decreasing during the next 200 to 300 ms (Figure 7B). Escape occurs in equal directions (Figures 7C and 7D), as shown in four examples where we followed their dispersal. To characterize this dispersion, we plotted the dispersion index (as defined in Methods S1), revealing two phases (Figure 7B): one below 100 ms, showing a rapid dispersion, followed by a second phase with less expansion. These two phases can be interpreted as follows: for the first one, trajectories escape from a well, then in the second phase, trajectories tend to be recaptured for a certain time in nodes, thus preventing a fast exploration of the network.

Figure 7

Local space exploration after trajectory resynchronization

Local network exploration revealed by spatial resynchronization of super-resolution SPTs.

(A) Reconstructed ER network using the GRA with nodes (black polygons) connected by segments (blue). Four nodes (I, II, III, and IV) are selected (red arrows).

(B) Average distance between the points of trajectories exiting a given node versus time after escape. The four nodes highlighted in (A) are also highlighted in red here, and the average line is presented in blue with its SD in blue shade (n = 258 nodes). This curve exhibits a fast (<50 ms) and a slow phase.

(C) Local exploration of neighboring nodes from trajectories (various colors) exiting from the chosen node for example nodes presented in (A).

(D) Distribution of velocities for trajectories exiting a node for the four nodes presented in (C). The rapid desynchronization of trajectories causes the mean velocity to rapidly decrease after exit (n is the number of trajectories exiting each node).

Trajectory resynchronization reveals novel local dynamics within dATL ER-tubules

We next analyzed how the local space exploration of trajectories is modified in COS-7 dATL cells (Figure S8H). Under normal conditions, trajectories are mostly located in nodes (Holcman et al., 2018), while here trajectories predominantly explore long tubules (Figure S8I) with average lengths of (Figure S8J), much longer than the 1 found for the tubules of WT cells. Inside these tubules, we found that trajectories exhibit a “stuttering” behavior around different positions that lasts for seconds. To characterize this behavior, we estimated several parameters such as the duration that a trajectory spent around a given position ms (Figure S8K), the transition time between different positions ms (Figure S8L), the length of a transition step (Figure S8M) and finally the SD around the stable positions (Figure S8N).

To conclude, following the reconstruction of networks using our algorithm, we were able to reposition and resynchronize SPTs. Using analyses of the ER lysosome, ER in normal COS-7 cells, and ER in COS-7 dATL cells, the algorithm revealed trajectories explored by the local space and the associated time scales.

Discussion and limitations of the methods

We present here a general method and the associated algorithms that can automatically characterize nanodomains where trajectories accumulate. Our approach generates graph representations of organelle networks from SPTs. Automatically finding nanodomains is useful to extract large statistics (size, energy of potential wells, mean residence time of particles) and compare their properties. Further, by quantifying the trajectories inside and outside these nanodomains, we could recover membrane organization, as well as determine the local redistribution dynamics of organelles and proteins. By reconstructing a graph of ER or lysosome networks from SPTs, we can recover molecular flow at the nanoscale level. We found here that nanodomains could be characterized as an attractive region (potential wells), and this generic representation suggests a universal mechanism of molecular stabilization that probably requires further investigation. Interestingly, these structures can be transiently remodeled in time, as revealed by the present time-lapse analysis.

Universality of high-density nanoregions characterized as potential wells

High-density nanoregions have now been associated with potential wells for several receptors and channels such as CaV (Heck et al., 2019), AMPAR (Hoze et al., 2012), glycine receptors (Masson et al., 2014), or GAGs (Hoze and Holcman, 2017; Floderer et al., 2018). Interestingly, some nodes of the ER can also be characterized as potential wells, which may reflect retention of luminal content or could be interpreted as a nanoregion allowing protein maturation. This representation suggests a generic membrane organization to retain particles (receptors, channels, proteins, etc.) in a field of force with long-range interactions.

Interestingly, the geometry of these regions and their energy profiles are independent of the experimental conditions, further confirming again their stability. Note that the physical nature of the potential well remains unclear (Holcman, 2013). The present method could be applied to analyze molecular crowding and the dynamics of nanodomains, thus clarifying processes relevant for phase separation at synapses (Feng et al., 2019). With the development of new labeling methods, improved fluorophores, and the ability to tag endogenous populations of molecules via CRISPR-Cas9, it will soon be possible to investigate phase separation at a population level and with SPT, to track endogenous dynamics, offering novel opportunities for the present approach.

Trafficking in networks

We show here that we can reconstruct a network from lysosome SPTs that resembles the ER network (Figure 5). This reconstructed network is further segregated into nodes and links, but low and high velocities are mixed (Figure 5A) compared with the reconstruction obtained from luminal proteins (Figure S7B). It is possible that lysosomes follow the cytoskeleton network, which is correlated with the ER (Lu et al., 2020). In addition, we find that the distribution of lysosome velocity follows a double exponential (Figure 5B) with fast () and slow () components. However, a more detailed analysis revealed that these velocities can be further subdivided into (1) confined motion (Figure 6) characterized by a residence time of s, and (2) deterministic motion between HDRs (Figures 6F and 6G), characterized by a distribution of velocity with an average of . It would be interesting to better characterize the switch between confined and rectilinear motion. Regions of deterministic velocities and those where diffusion can be found are often not well separated, suggesting that lysosomes can use various modes of transport, independently of the subregions where there are located. We found, however, some regions characterized by a high density of trajectories, with a low velocity, suggesting that there are possible trapping mechanisms to retain lysosomes in specific subregions of the ER, possibly at exit sites (Manley et al., 2008). This mode of motion is quite different from the internal motion inside the ER lumen or on its membrane: in the first case, the node-tubule topology is associated with diffusion-drift while lysosomes may be trapped in interaction with the ER.

Future work should reveal interaction times between lysosomes and the ER. By applying our algorithms to different cells and organelles, we have shown that the boundaries and dynamics of subcellular interactions can be revealed from large SPT datasets. The automated algorithms presented here can be applied to analyze hundreds of thousands of trajectories and to study nanodomains with almost no human interventions and are available as an ImageJ plugin.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources should be directed to and will be fulfilled by the lead contact, David Holcman (david.holcman@ens.fr).

Materials availability

This study did not generate new unique reagents.

Method details

Diffusion model, velocity, vector fields and empirical estimators

In the Smoluchowski’s limit of the Langevin equation (Schuss, 2009), the position of a molecule is described bywhere is a field of force, is a white noise, is the friction coefficient (Schuss, 2009) and is the diffusion coefficient. At a coarser spatio-temporal scale, the motion can be coarse-grained as a stochastic process (Hoze et al., 2012; Hoze and Holcman, 2014)where is the drift field and the diffusion matrix. The effective diffusion tensor is given by ( denotes the transposition) (Schuss, 2009, 2010). The drift and diffusion fields from Equation 3 can be recovered from SPTs acquired at any infinitesimal time step by estimating the conditional moments of the trajectory displacements (Schuss, 2009; Friedrich and Peinke, 1997; Siegert et al., 1998; Hoze and Holcman, 2014; Hozé and Holcman, 2017)

The notation represents averaging over all trajectories that are passing at point at time . To estimate the local drift and diffusion coefficients at each point of the membrane and at a fixed time resolution , we use a procedure based on a square grid.

The local estimators to recover the vector field and diffusion tensor (Parutto et al., 2019) consist in grouping points of trajectories within a lattice of square bins centered at and of width . For an ensemble of two-dimensional trajectories with the number of points in trajectory and successive points recorded with an acquisition time .,the discretization of Equation 4 for the drift in a bin centered at position iswhere and is the number of points falling in the square . Similarly, the components of the effective diffusion tensor are approximated by the empirical sums

The centers of the bin and their size are free parameters that are optimized during the estimation procedure.

Cosine-filtered estimation of the diffusion coefficient and drift

To increase the accuracy of the diffusion and drift maps, we weighted the points in the moving windows with a cosine function (would also be possible to use wavelets). In that case, the new estimator for the drift field is nowwith the number of points of the trajectories falling in the disk of radius and centered at . The weight of a displacement starting at with respect to the disk is given bywith the Euclidean norm. In that case, we can choose a refined grid with bin size , where is a scaling factor. The role of the cosine weights is to decrease continuously the influence of the points falling near the boundary.

Similarly, the generalized formula for the effective diffusion tensor are given by the weighted sumswhere the weights are given by Equation (9).

Local point density estimation

The local density of points can be determined using a procedure similar to the drift or diffusion estimation, dividing the image plane into a square bin . We then compute for each square of centered at where is the number of trajectory points falling into the bin centered at . In practice, it usually helps to smooth this density estimation by applying a local average using a small kernel with .

Estimating potential well parameters

In this subsection, we present the estimators for the two main parameters of potential wells: the extent of their boundary and their associated energy (Hoze et al., 2012; Masson et al., 2014; Parutto et al., 2019). We recall that the diffusion coefficient inside a well is considered to be constant and the energy of the well given by where is the attraction coefficient and the diffusion coefficient. The correlations between potential wells characteristics for CaV2.1 and ER luminal data are presented in Figures S6A–S6C, S6E–S6G. In practice, for the potential wells reported in Figures 2 and 3, we filtered out the wells with energies .

Residence time of particles inside a well

The potential well model allows to estimate the residence time using the classical escape time formula (Schuss, 2009; Hoze et al., 2012) for a circular wellwith the radius of the well, its attraction coefficient and its diffusion coefficient. In the case of an elliptic well, we obtain an approximate circular boundary using the harmonic mean of the semi-axes , where and are the large and the small-axes lengths respectively. This approximation holds for .

Parabolic potential well representation

To extract the energy of a potential well, we consider the basin of attraction of a truncated elliptic parabola with the associated energy functionwhere , is the center of the well, are the elliptic semi-axes lengths and the elliptic boundary is defined by

Recovering the center

The center of the nanodomain is estimated as the center of mass of the cloud of points falling inside the HDR. We use the empirical averaging formulawhere is the total number of points such that .

Covariance matrix

We use the sample estimator of the covariance matrix Σ defined for a cloud of two-dimensional points as

Confidence ellipse estimation

We define the boundary of the well as the % confidence ellipse of the associated Gaussian density distribution of center and covariance . Using the singular value decomposition method, we decompose the covariance matrix aswhere are unitary matrices and is diagonal. The values in represent the variance in each dimension along the principal components. For a Gaussian density, the values of follow a Chi-Squared distribution with degrees of freedom. Therefore the semi-axes lengths can be obtained at the confidence level from aswith is given by for a Chi-Squared distribution with two degrees of freedom (for example , and ). Finally, the orientation of the ellipse is defined by the angle

Maximum likelihood estimators (MLE) based on an Ornstein-Uhlenbeck model

Using the potential well representation from Equation 12 in the stochastic model presented in Equation 2 leads to a truncated Ornstein-Uhlenbeck process, centered at , with an attraction coefficient and diffusion coefficient . The probability density function for of observing two successive positions of the same trajectory and , separated by a time step is given bywithandwhich we rewrite asand

The log-likelihood function of observing the successive pairs , , possibly from various trajectories, is given by

The corresponding maximum likelihood estimators for and , and are obtained by solving the system of equationsfrom which we obtain the empirical estimators for the drift coefficient (Tang and Chen, 2009)and the diffusion coefficient:

Note that the parameter in Equation 27 can be computed from the estimator .

Hybrid density-drift algorithm

In this sub-section, we present two variants of an algorithm to detect the main characteristics of a potential well from some observed trajectories: the center , the semi-axes lengths , the orientation , the attraction coefficient , the diffusion coefficient and the energy .

Fixed spatial scale hybrid density-drift algorithm

Initiation

Search for high-density regions: the image is partitioned by a grid with square bins of size . from which we compute the density map (Equation 10). We then select the bins from with the highest density as possible “seed” regions containing a potential well.

Iterations

For each seed region obtained in the initiation step, we apply an iterative procedure that is going to consider increasingly larger square neighborhoods around this region. For each iteration , we keep only the trajectories contained inside the square of size and centered at point . The point is the center of mass of points falling in the square (Equation 14) and is the center of the initial high-density bin. The elliptic semi-axes , are computed as the confidence ellipse (Equation 17, ellPerc parameter as defined in Methods S1) from the covariance matrix (Equation 15) and the angle is the orientation of the ellipse (Equation 18). These parameters define the elliptic boundary of the well at iteration

The attraction coefficient and diffusion coefficient are computed from Equations 26 and 27 respectively, for the trajectories contained inside . Specifically, we obtain from

We repeat this procedure times, with for the spatial parameter and the maximum region size , defined by the user.

Termination

This step consists in selecting the best iteration among : we evaluate for each iteration the likelihood score defined by Equation 24 but computed for sub-trajectories falling in the ring formed by the ellipses and . In practice, we filter out rings with ringMinPts trajectories which is a value specified by the user (see Methods S1). The best iteration is selected as the first local maximum of the curve .

Multiscale hybrid algorithm

We generalize the hybrid density-drift algorithm defined above for a fixed spatial scale, by now varying the grid size , in the range  selected by the user. The purpose of this algorithm is to select the optimal size that maximizes the likelihoodwhere is the iteration that maximizes across all the spatial scales , . In practice the range of grid sizes is specified by a minimal value , a step and a maximum value as presented in Methods S1.

Density-based algorithm

The Density-Based Algorithm uses the density of points estimated around the local density maximum of a high-density region and was first presented in (Parutto et al., 2019). The algorithm uses the level set of a Gaussian distribution of points inside a potential well. We define the level set with respect to a local maximum as the ensemble of all trajectory points falling in bins with a density greater than :where is the empirical point density, estimated over the bins of a square grid (Equation 10) and is a density threshold. For the points located in , the center of the distribution is approximated by the empirical estimators based on Equation 14 but restricted to the points in :with the number of points in the ensemble and . Similarly, we extend the estimator for the covariance matrix (Equation 15) to

We now define the density-based algorithm:

Search for high-density regions: the image is partitioned by a grid with square bins of size from which we compute the density map (Equation 10). We then select the bins from with the highest density as possible regions containing a potential well.

For each selected high-density bin, we initialize the well center at the center of the bin. We then construct a refined grid centered at and with bin size (parameter locGridDx as defined in Methods S1). In this grid, we compute the centers , for different values of : , selected by the user. The refined center is obtained as the center of mass of the centers for .

We then apply an iterative procedure that considers increasingly larger concentric annulus of center , width and radius for . Where the number of iterations is determined based on a minimal and maximal distances defined by the user (see Methods S1). For each iteration , we compute the confidence ellipse (see sub-section Method details) obtained from the covariance matrix (Equation 33) computed only for the points that fall in the annulus of radius . We then search for the iteration that maximizes the ratio (in practice, we limit the maximal distance for finding to ratMaxDist, which value is specified by the user (see Methods S1)) and use it to define the refined distance to the centerwhere , that transforms an ellipse into a circle with the same center. Finally, we compute the density of points falling in the annulus of radius based on the refined distance measure .

We select the first iteration such that : it is the first iteration where the derivative of the density with respect to the distance to the center, stops decreasing. This criteria is more stable on empirical data than searching for the minimum of the density (see Figures 3 and 4 panel B3 of (Parutto et al., 2019)). The elliptic boundary of the well is centered at , has semi-axis lengths given by and and orientation . We then use the ML estimator (Equation 27) to estimate the constant diffusion coefficient inside . Finally, to compute we use the diagonal form of the covariance matrix estimated (Equation 33) for all the points falling in :and estimate

Drift-based algorithm

The drift based algorithm uses an error function in the space of the vector field to estimate the characteristics of a well and was introduced in (Heck et al., 2019). Its principle is as follows:

Search for high-density regions: the image is partitioned by a grid with square bins of size from which we compute the density map (Equation 10). We then select the bins from with the highest density as possible regions containing a potential well.

For each region selected in the initiation, we apply the following iterative procedure for :where are the centers of the bins from that are contained inside the ellipse .

We select only the sub-trajectories falling into a square centered at and of size . We take to be the center of the high-density bin.

We estimate the elliptic well boundary as the confidence ellipse (parameter ellPerc as discussed in Methods S1) from the cloud of points in . Where is a parameter selected by the user (usually or 99).

We then compute a new grid centered at , that we use to estimate the local drift map (Equation 6) and estimate the attraction coefficient of the well using the least-square regression formula

Finally, we estimate the quality of the well (parabolic index) based on the residual least square error:

The index is defined such that for a drift field generated by a parabolic potential well and for a random drift vector field as observed for diffusive motion.

The number of iterations is given by where is the maximum size of the region to consider and is given by the user.

We select the iteration that minimizes the parabolic index : . We estimate the diffusion coefficient inside the well using the local estimator (Equation 4) for all the displacements inside the ellipse .

Sliding window analysis to study the stability of the wells over time

To determine the stability of the potential wells, we use a non-overlapping sliding window of 20 s (Heck et al., 2019), to recover the ellipse characteristics, as shown on different examples in Figure S6A and S6C. When a well disappears in a given time window, but reappears latter, we kept the well for the entire period.

Reconstructing the graph of the network explored by SPTs

We describe here a new variant of the algorithm to reconstruct a graph from SPTs exploring a network. This algorithm is based on the graph reconstruction algorithm introduced in (Holcman et al., 2018) and exploits the heterogeneous distribution of points caused by trajectories spending more time inside nodes than in tubules.

Recursive velocity based graph reconstruction algorithm of a network explored by SPTs

The recursive graph reconstruction algorithm first generates an ensemble of points from slow trajectory displacements based on a maximum displacement length threshold and then recursively uses the dbscan algorithm (Ester et al., 1996) to cluster these points based on their local density. The algorithm requires specification of two ensembles of parameters:

An ensemble of distances (in ) defining the neighborhood radius around points.

An ensemble of counts defining the numbers of points required in the neighborhood to form a cluster (Ester et al., 1996).

A pair of these two parameters for any defines a local density (points/ m2) inside each cluster. The values of and depend on the local number of recorded trajectories and can vary inside the image. For each dataset, these values can be determined such that the computed clusters overlap with the structure of the organelle formed by the trajectories.

We now present the steps of the algorithm:

We form the ensemble of points belonging to low-velocity trajectory fragments .

We apply the dbscan procedure with parameters to obtain a first ensemble of clusters from the points in .

We then perform a “recursive” step where we refine the initial clusters possessing more than maxClustNpts points. For any such cluster :

We iteratively re-apply the dbscan algorithm with the more stringent parameter pair for and and replace the initial cluster with the resulting sub-cluster(s). We continue iterating over the generated sub-clusters until they all possess less than maxClustNpts points.

We then approximate the boundary of each cluster either by its minimum volume ellipsoid (Todd, 2016) or its convex hull polygon and assign each point discarded in step 1 to the cluster in which they fall if possible.

Finally, we merge any overlapping pair of clusters by computing the boundary of the combined ensemble of points (either elliptic or the convex hull) and we iterate this procedure until no more clusters overlap.

This first step allows to find the nodes of the network. In the second step, we define tubules by constructing a connectivity matrix of size where is the number of first or second order trajectory segments connecting nodes and . Specifically, we increment for each data point (, ) in the following cases:

is located in node and in node

is located in node , does not belong to any node and is located in node (in this case ).

Lysosome analysis

Trajectories analysis for lysosome SPTs

To study the dynamics of lysosomes, we plotted the distribution of instantaneous velocities, computed from each trajectory displacement by

where s. We approximate the distribution of instantaneous velocities using a two exponential model obtained by fitting to the distribution using MATLAB’s fit toolbox. The density and diffusion maps were computed using the estimators from Equation 4 described above.

Local high-density region analysis: Ellipse approximation of the boundary

High-density regions are extracted from trajectories as follows: we construct a density map (Equation 10) based on a square grid with bin size nm. From this density map, we select only the highest density bins and in case multiple such bins appear within a distance of two squares of each other, only the one with the highest value was kept. For each selected bin of center , we computed a refined density map of size squares, centered at and with bin size nm. From this local map, we collected trajectory points falling into the bins that have a density of the maximal bin value and use them to estimate the elliptic boundary of the region from the confidence ellipse (see sub-section Method details). Finally, when a pair of ellipses overlap, we replaced them by the ellipse computed over their combined ensemble of points and iterated this procedure until no more overlaps could be found.

Transient confinement detection

To detect transient confinement periods along individual trajectories, we used the following procedure: for each point of a trajectory, we considered the ensemble of its successors , where initially is set by the user. We then computed the center of mass and checked that all the points have a distance to the center of mass , for a chosen distance threshold ( is the Euclidean norm). We then iterate the procedure, considering increasingly larger ensembles of successors until either reaching the end of the recorded trajectory or when the next point do not fall into the neighborhood of the center of mass . The confinement duration is then computed by considering the difference in time between the two endpoints of the ensemble. Finally, the spring constant and diffusion coefficient of the confinement is obtained by applying the Ornstein-Ulenbeck maximum likelihood estimators (Tang and Chen, 2009; Parutto et al., 2019), where the OU-process is given by

ImageJ plugin

The present method and algorithms are implemented into an ImageJ plugin called "SPTAnalysis". The plugin allows to reconstruct the various maps (trajectories, density, drift, diffusion), detect potential wells and reconstruct the graph associated to trajectories. It allows to extract various statistics such as the distribution of diffusion coefficients or the energy and the size of potential wells. The algorithms performances are summarized in Figure S9 and Method S1.

Quantification and statistical analysis

No statistical tests have been used in this paper. Data quantification is reported as average standard deviation, where the meanings of are indicated in the corresponding legends.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Software and algorithms
MATLAB2019b	mathworks	https://uk.mathworks.com/products/matlab.html
Python3.8.3	Python Software Foundation	https://www.python.org/
Fiji2.1.0	Open Source	https://fiji.sc/
SPTAnalysis1.1	This Paper	https://doi.org/10.5281/zenodo.6862643