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

Summary 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.


In brief
The rapid development of superresolution single-particle trajectories microscopy requires high-throughput analysis techniques. Parutto et al. propose methods to reconstruct and quantify organelle networks and membrane reorganization from trajectories of any molecules, receptors, or channels of interest that can be tagged. The methods are available as an ImageJ plugin (SPTAnalysis).

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(Wu et al., , 2018Petkovic et al., 2014), the 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.
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-trajectoriesusing-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.

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 DX = Xðt + DtÞ À XðtÞ, where XðtÞ is the  (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.
Cell Reports Methods 2, 100277, August 22, 2022 3 Article ll OPEN ACCESS position of the particle at time t (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 ðDXÞ 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: 1. Automatic determination of bins with the highest density of points ( Figure S3Ba). 2. For each such bin, we iterate over square regions of increasing width w k around the bin center. For each iteration k, we apply a principal-component analysis to estimate the semi-axis a k ; b k , the center c k , and a score L (Equation 24) based on the points falling in the square of size w k (Figures S3Bb and S3Bc). We iterate until we reach the maximum size of the well set by the user ( Figure S3Bc). 3. 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 E = 0 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 D47 and CaV2.1 + 47, 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 a = 143 ± 51 and b = 104 ± 33 nm for Cav2.1 D47. These parameters are similar for Cav2.1 + 47, while they were smaller for endogenous Cav2.1 with a = 100 ± 47 and b = 73 ± 31 nm ( Figure 2D). The diffusion coefficient was D = 0:091 ± 0:052 mm 2 =s for the D47 variant and D = 0:087 ± 0:045 mm 2 =s for the + 47 variant and was smaller for endogenous CaV2.1, with D = 0:069 ± 0:051 mm 2 =s (Figure 2E). Interestingly although the associated energy was $ 3:9 kT for the two variants and 3:3 kT 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 A 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).
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  Article ll OPEN ACCESS a well across successive time windows, we consider that two wells identified at times t k and t k + 1 are the same if the distance between their centers is less than 200 nm. The ensemble of consecutive times ðt q ; .; t r Þ, where a well is first detected at time t q and disappears at time t r + 1 , is used to define the stability duration t = t r À t q . 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 D47 and + 47 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. 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 $ 0:3 mm 2 =s, and attraction coefficients $ 1 mm 2 =s, corresponding to energies $ 2:9 kT and residence times $ 90 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.

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 groundtruth 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 D app = 1:56 ± 0:83 mm 2 =s.
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 d nodes = 0:92 ± 0:32 mm ( Figure 4G) and the nodes area A nodes = 0:15 ± 0:12 mm 2 ( 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  Figure 5A). These trajectories are characterized by a distribution of instantaneous velocities in the range ½0 À 3:5 mm=s ( 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 fðvÞ of velocities (Figure 5B) can be fitted by a sum of two exponentials: where a best fit approximation gives v 0 = 0:043 mm=s (95% confidence interval ½0:041; 0:044) and v 1 = 0:388 mm=s (95% confidence interval ½0:370;0:407) (coefficient of determination R 2 = 0:999), with A = 0:780 and B = 0:075. 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 0:5 mm=s, and 12.7% are in the range of [0.5-1.5] mm=s.
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 ( 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 a = 516 ± 196 (large) and b = 278 ± 143 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 v = 1:03 ± 0:32 mm=s, 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).
To recover the size of the confinement areas, we fitted ellipses over these regions and obtained average semi-axes lengths (Figure 6D) of a = 232 ± 77 (large) and b = 94 ± 47 nm (small). Furthermore, this approach allowed us to estimate the confinement strength l 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 l = 0:123 ± 0:025 s À 1 ( Figure 6E), associated with an average local diffusion coefficient of D = 0:032 ± 0:020 mm 2 =s ( 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 t = 5:35 s ( Figure 6G). The average residence time of lysosomes in these regions is t = 30 ± 12 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 v e = 30:2 ± 10:2 mm=s 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. 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 5:4 ± 2:44 mm ( Figure S8J), much longer than the $ 1 mm 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 t ll = 79 ± 76 ms ( Figure S8K), the transition time between different positions t t = 27 ± 15 ms ( Figure S8L), the length of a transition step D ll = 0:53 ± 0:45 mm ( Figure S8M) and finally the SD around the stable positions SD = 0:14 ± 0:07 mm ( 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 . 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 Article ll OPEN ACCESS 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 ($ 0:388 mm=s) and slow ($ 0:043 mm=s) 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 $ 5 s, and (2) deterministic motion between HDRs ( Figures 6F and 6G), characterized by a distribution of velocity with an average of 1:03 mm=s. 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 to-pology 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
Detailed methods are provided in the online version of this paper and include the following:

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.
Data and code availability d All data reported in this paper will be shared by the lead contact upon request. d All original code has been deposited at github.com/holcman-lab/SPTAnalysis and is publicly available as of the date of publication. DOIs are listed in the key resources table. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.

METHOD DETAILS
Diffusion model, velocity, vector fields and empirical estimators In the Smoluchowski's limit of the Langevin equation (Schuss, 2009), the position XðtÞ of a molecule is described by where FðXÞ is a field of force, W is a white noise, g is the friction coefficient (Schuss, 2009) and D 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 aðXÞ is the drift field and BðXÞ the diffusion matrix. The effective diffusion tensor is given by DðXÞ = 1 2 BðXÞB T ðXÞ (: T denotes the transposition) (Schuss, 2009(Schuss, , 2010. The drift and diffusion fields from Equation 3 can be recovered from SPTs acquired at any infinitesimal time step Dt by estimating the conditional moments of the trajectory displacements DX = Xðt + DtÞ À XðtÞ The notation E½,jXðtÞ = x represents averaging over all trajectories that are passing at point x at time t. To estimate the local drift aðXÞ and diffusion coefficients DðXÞ at each point X of the membrane and at a fixed time resolution Dt, we use a procedure based on a square grid. REAGENT  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 Sðx k ; DxÞ centered at x k and of width Dx. For an ensemble of N two-dimensional trajectories fX i ðt j Þ = ðx ð1Þ i ðt j Þ; x ð2Þ i ðt j ÞÞ; i = 1.N; j = 1.M i g with M i the number of points in trajectory X i and successive points recorded with an acquisition time t j + 1 À t j = Dt.,the discretization of Equation 4 for the drift aðx k Þ = ða ð1Þ ðx k Þ; a ð2Þ ðx k ÞÞ in a bin centered at position x k is ( Equation 6) where u = 1::2 and N k is the number of points x i ðt j Þ falling in the square Sðx k ;rÞ. Similarly, the components of the effective diffusion tensor Dðx k Þ are approximated by the empirical sums The centers of the bin and their size Dx 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 now a ðuÞ ðx k Þz with N k the number of points of the trajectories falling in the disk Dðx k ; rÞ of radius r and centered at x k . The weight of a displacement starting at X i ðt j Þ with respect to the disk Dðx k ; DxÞ is given by w i;j ðx k ; rÞ = cos p 2 kX i ðt j Þ À x k k r ; (Equation 9) with k:k the Euclidean norm. In that case, we can choose a refined grid Sðx k ; ðDxÞ 0 Þ with bin size ðDxÞ 0 = Dx=l sc , where l sc is a scaling factor. The role of the cosine weights w is to decrease continuously the influence of the points falling near the boundary. Similarly, the generalized formula for the effective diffusion tensor Dðx k Þ are given by the weighted sums D ðu;vÞ ðx k Þz where the weights w are given by Equation (9).

Local point density estimation
The local density of points r can be determined using a procedure similar to the drift or diffusion estimation, dividing the image plane into a square bin Sðx k ; DxÞ. We then compute for each square of s centered at x k where N k is the number of trajectory points falling into the bin centered at x k . In practice, it usually helps to smooth this density estimation by applying a local average using a small d3d kernel with d $ 1; 3; 5.

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 E = A=D where A is the attraction coefficient and D 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 E > 10 kT. with r the radius of the well, A its attraction coefficient and D its diffusion coefficient. In the case of an elliptic well, we obtain an approximate circular boundary using the harmonic mean of the semi-axes r = ffiffiffiffiffi ffi ab p , where a and b are the large and the smallaxes lengths respectively. This approximation holds for azb.

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 function where X = ½x ð1Þ ;x ð2Þ , m = ½m ð1Þ ; m ð2Þ is the center of the well, a; b are the elliptic semi-axes lengths and the elliptic boundary is defined by :

(Equation 13)
Recovering the center m The center of the nanodomain B is estimated as the center of mass of the cloud of points falling inside the HDR. We use the empirical averaging formula where N is the total number of points such that X i˛B .

Covariance matrix S
We use the sample estimator of the covariance matrix S defined for a cloud of N two-dimensional points X i = ðx with j x is given by Pðv < j x Þ = x for a Chi-Squared distribution with two degrees of freedom (for example j 99 = 9:210, j 95 = 5:991 and j 90 = 4:605). Finally, the orientation 4 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 m, with an attraction coefficient l and diffusion coefficient s = ffiffiffiffiffiffi ffi 2D p . The probability density function pðXðt j + 1 Þ Xðt j ÞÞ for j = 1::ðM À 1Þ of observing two successive positions of the same trajectory Xðt j Þ and Xðt j + 1 Þ, separated by a time step t j + 1 À t j = Dt is given by

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 m, the semi-axes lengths a R b, the orientation 4, the attraction coefficient A, the diffusion coefficient D and the energy E.
Fixed spatial scale hybrid density-drift algorithm Initiation Search for high-density regions: the image is partitioned by a grid G Dx with square bins of size Dx. from which we compute the density map r Dx ðxÞ (Equation 10). We then select the bins from r Dx ðxÞ with the highest d% 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 k = 1::K, we keep only the trajectories contained inside the square G k;Dx of size ½ð2k + 1Þ 3ð2k + 1ÞðDxÞ 2 and centered at point m k À 1 . The point m k is the center of mass of points falling in the square G k (Equation 14) and m 0 is the center of the initial high-density bin. The elliptic semi-axes a k , b k are computed as the x% confidence ellipse (Equation 17, ellPerc parameter as defined in Methods S1) from the covariance matrix C k (Equation 15) and the angle 4 k is the orientation of the ellipse (Equation 18). These parameters define the elliptic boundary of the well at iteration k The attraction coefficient A k and diffusion coefficient D k are computed from Equations 26 and 27 respectively, for the trajectories contained inside ε k . Specifically, we obtain A k from We repeat this procedure K times, with K = AE Ms Dx Ç for the spatial parameter Dx and the maximum region size M s , defined by the user. Termination This step consists in selecting the best iteration among K: we evaluate for each iteration k > 1 the likelihood score L k = Lðm k ; l k ; s k X 1 ; ::X p˛εk À ε k À 1 Þ defined by Equation 24 but computed for sub-trajectories falling in the ring formed by the ellipses ε k À 1 and ε k . In practice, we filter out rings with < ringMinPts trajectories which is a value specified by the user (see Methods S1). The best iteration k Ã is selected as the first local maximum of the curve L k .

Multiscale hybrid algorithm
We generalize the hybrid density-drift algorithm defined above for a fixed spatial scale, by now varying the grid size Dx, in the range Dx 1 < Dx 2 < :: < Dx N selected by the user. The purpose of this algorithm is to select the optimal size Dx i Ã that maximizes the likelihood is the iteration that maximizes L across all the spatial scales Dx i , i = 1::N. In practice the range of grid sizes is specified by a minimal value Dx min , a step Dx step and a maximum value Dx max 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 G a with respect to a local maximum M Ã as the ensemble of all trajectory points falling in bins with a density greater than aM Ã : G a = fX i such that r e ðxÞ > aM Ã g; ( Equation 31) where r e is the empirical point density, estimated over the bins of a square grid (Equation 10) and a˛½0; 1 is a density threshold. For the points X i = ðx We now define the density-based algorithm: Initiation Search for high-density regions: the image is partitioned by a grid G Dx with square bins of size Dx from which we compute the density map r Dx ðxÞ (Equation 10). We then select the bins from r Dx ðxÞ with the highest d% density as possible regions containing a potential well.

Iterations
For each selected high-density bin, we initialize the well centerm 0 at the center of the bin. We then construct a refined grid centered at m 0 and with bin size Dx 0 < Dx (parameter locGridDx as defined in Methods S1). In this grid, we compute the centers m 0;ak , for different values of a: a 1 < . < a N , selected by the user. The refined center m 0 is obtained as the center of mass of the centers m 0;ak for k = 1::N.
We then apply an iterative procedure that considers increasingly larger concentric annulus of center m 0 , width Dr and radius r k for k = 1::K. Where the number of iterations K is determined based on a minimal r 1 = r min and maximal r K = r max distances defined by the user (see Methods S1). For each iteration k, we compute the confidence ellipse ε k = ðm 0 ; a k ; b k ; 4 k Þ (see sub-section Method details) obtained from the covariance matrix S k (Equation 33) computed only for the points that fall in the annulus of radius r k . We then search for the iteration r Ã that maximizes the ratio C v ðr k Þ = ffiffiffiffiffiffiffiffiffiffiffiffi a k =b k p (in practice, we limit the maximal distance for finding C v to r k < ratMaxDist, which value is specified by the user (see Methods S1)) and use it to define the refined distance to the center where k = C v ðr Ã Þ, that transforms an ellipse into a circle with the same center. Finally, we compute the density of points N e ðr k Þ falling in the annulus of radius r k based on the refined distance measure r e .

Termination
We select the first iteration k Ã such that N e ðr Ã k Þ > N e ðr k Ã À 1 Þ: 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 m 0 , has semi-axis lengths a Ã ; b Ã given by a Ã = r k Ã and b Ã = kr k Ã and orientation 4 k Ã . We then use the ML estimator (Equation 27) to estimate the constant diffusion coefficient D inside ε Ã . Finally, to compute A Ã 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: Initiation Search for high-density regions: the image is partitioned by a grid G Dx with square bins of size Dx from which we compute the density map r Dx ðxÞ (Equation 10). We then select the bins from r Dx ðxÞ with the highest d% density as possible regions containing a potential well.

Iterations
For each region selected in the initiation, we apply the following iterative procedure for k = 1::K: (a) We select only the sub-trajectories falling into a square G k ðm k ; DxÞ centered at m k À 1 and of size ð2k + 1ÞDx 3 ð2k + 1ÞDx. We take m 0 to be the center of the high-density bin. (b) We estimate the elliptic well boundary ε k = ½m k ; a k ; b k ; 4 k as the x% confidence ellipse (parameter ellPerc as discussed in Methods S1) from the cloud of points in G k ðc k ; DxÞ. Where x is a parameter selected by the user (usually 90; 95 or 99). (c) We then compute a new grid G Dx ðm k Þ centered at m k , that we use to estimate the local drift map (Equation 6) a k ðXÞ = ½a (d) Finally, we estimate the quality of the well (parabolic index) based on the residual least square error: The index S k˛½ 0; 1 is defined such that S k /0 for a drift field generated by a parabolic potential well and S k /1 for a random drift vector field as observed for diffusive motion.
The number of iterations is given by N = bw max =Dxc where w max is the maximum size of the region to consider and is given by the user.

Termination
We select the iteration k Ã that minimizes the parabolic index S: k Ã = argmin k = 1.K S k ða k ; A k Þ. We estimate the diffusion coefficient inside the well using the local estimator (Equation 4) for all the displacements inside the ellipse ε k Ã .
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 v th 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: 1. An ensemble of distances R = fR u ; u = 1::Ug (in mm ) defining the neighborhood radius around points.
2. An ensemble of counts N = fN v ; v = 1::Vg defining the numbers of points required in the neighborhood to form a cluster (Ester et al., 1996).
A pair of these two parameters ðR u ; N v Þ for any u; v defines a local density Nv Ru (points/ mm 2 ) inside each cluster. The values of R and N 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: 1. We form the ensemble of points belonging to low-velocity trajectory fragments S = 2. We apply the dbscan procedure with parameters ðR U ; N 1 Þ to obtain a first ensemble of K clusters c 1 ; .; c K from the points in S. 3. We then perform a ''recursive'' step where we refine the initial clusters possessing more than maxClustNpts points. For any such cluster c k : (a) We iteratively re-apply the dbscan algorithm with the more stringent parameter pair ðR u ; N v Þ for u = ðU À 1Þ::1 and v = 2::V 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. 4. 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. 5. 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 K nodes of the network. In the second step, we define tubules by constructing a connectivity matrix C of size K3K where c i;j ð1 % i; j % KÞ is the number of first or second order trajectory segments connecting nodes i and j. Specifically, we increment c i;j for each data point X k ðt l Þ (1 % k % N t , 0 % l < M k À 1) in the following cases: 1. X k ðt l Þ is located in node i and X k ðt l + 1 Þ in node j 2. X k ðt l Þ is located in node i, X k ðt l + 1 Þ does not belong to any node and X k ðt l + 2 Þ is located in node j (in this case 0 % l < M k À 2).

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 Xðt + DtÞ À XðtÞ by v = Xðt + DtÞ À XðtÞ Dt ; (Equation 39) where Dt = 1:5 s. We approximate the distribution of instantaneous velocities using a two exponential model obtained by fitting fðvÞ = Ae v 1 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 Dx = 480 nm. From this density map, we select only the 5% 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 c, we computed a refined density map of size 535 squares, centered at c and with bin size Dx 0 = 200 nm. From this local map, we collected trajectory points falling into the bins that have a density > 80% of the maximal bin value and use them to estimate the elliptic boundary of the region from the 95% 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 Xðt j Þ of a trajectory, we considered the ensemble of its successors e tj;n = fXðt j Þ;.; Xðt j + n Þg, where initially n = N nh is set by the user. We then computed the center of mass m tj;n and checked that all the points Xðt k Þ˛e tj ;n have a distance to the center of mass kXðt j Þ À m tk ;n k < R nh , for a chosen distance threshold R nh (k:k is the Euclidean norm). We then iterate the procedure, considering increasingly larger ensembles of successors n = fN nh ; N nh + 1; .; N nh + Kg until either reaching the end of the recorded trajectory or when the next point Xðt j + n + 1 Þ do not fall into the neighborhood of the center of mass m tj ;n . The confinement duration is then computed by considering the difference t j + n À t j in time between the two endpoints of the ensemble. Finally, the spring constant l and diffusion coefficient D 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.   Figure S1: Evaluation of diffusion field estimation, related to STAR Methods.
A. Perlin noise image used to generate a non-uniform diffusion field. B. The non-uniform diffusion field generated from A by setting a pixel size of 64.5 nm, D min = 1e −4 µm 2 /s and D max = 0.1µm 2 /s. C. Individual Brownian trajectories simulated from the diffusion field presented in B. D. Diffusion maps estimated from the trajectories presented in C without filtering (top, parameters ∆x = 100 nm, minP ts = 5) and with cosine-filtering (bottom, parameters ∆x = 50 nm, minP ts = 10, r f ilt = 500 nm). E. Heat maps evaluating the quality of the reconstructed diffusion field for a uniform field with D = 0.05µm 2 /s obtained over 100simulations. The top line (blue shade) reports the average error (see eq.1 of Data 1) between the reconstructed ground truth field while the bottom line (brown shade) shows the average percentage of the field area recovered. F. Same as E but for the nonuniform field presented in B. Black squares correspond to conditions where not enough data could be recovered.  Figure S2: Comparing the diffusion maps for CaV2.1 SPTs from three estimators, related to STAR Methods. A. Individual CaV2.1+47 trajectories overlaid on top of the synaptotagmin pre-synaptic terminal marker's signal. B. Diffusion map computed from the classical estimator (Eq. 7 main text). C. Diffusion map computed with the cosine-estimator (Eq. 9 Main text). The cosine estimator is computed over a sliding window with a moving disk of radius r f ilt = 0.12 µm and a grid size ∆x = 0.01 µm. D. Map computed using a local average filter (using a 3 × 3 window around each bin). E. Comparison of the tree approaches from a representative regions. F. Violin plots of the diffusion coefficients computed in each bin for the three methods (n represents the number of bins).  Figure S3: Reconstruction algorithm of nanodomains from superresolution SPTs, related to Figures 2 and 3. A. a-Simulated trajectories attracted in a well of known boundary. b-density map c-drift field. B. Algorithm to reconstruct the center, boundary and energy of the nanodomain characterized as a potential well. Initiation of the center c 0 with an ellipse γ 0 and the center of mass µ 0 . The error between the vector field and the reconstructed field is measured by the score S 0 . The spatial resolution is fixed by ∆x. k-iteration step: the domain is enlarged to Γ k , the center is revaluated at point c k and the error S k is recomputed.
In the termination step, the optimal value k * is selected that minimizes the error S k * , leading to the optimal center and ellipse γ k * , for which the energy is computed as the ratio A/D of the field to the diffusion coefficient. C. Three examples of potential wells reconstructed from trajectories. D as a function of the energy of the well and depending on the acquisition time. The simulation, well detection algorithm and its parameters are the same as in Table S3.  Table S5).  Network segmentation algorithm Step 1: Clustering low-velocity displacements Step 2: Constructing high-velocity regions Step 3: Constructing the ER graph A. Individually color-coded ER-luminal trajectories. B. the same trajectories as in A but where each displacement is color-coded by its instantaneous velocity. C. Distribution of instantaneous velocities color-coded according to their amplitude as in B. D-G. Extraction of the network nodes from low-velocity regions: a threshold v th is selected on the distribution of instantaneous velocities (D) and the corresponding displacements are shown (E), these points are clustered using the dbscan algorithm (or a recursive version of it, see Method) and the convex hull of each cluster is computed to obtain its boundary (F). Finally, the corresponding boundaries are overlaid on top of the original trajectories (G). H-K. Extraction of the graph links from high-velocity jumps: we now select only the displacements with velocities above the threshold used in step 1 (H) shown in (I). We then filter them to keep only those that either: i. have their endpoints that fall into two separate ellipses (single-jump) extracted from step 1 or ii. successive displacements where the first point falls in an ellipse, the middle point is outside and the third point fall in a different ellipse (double-jump) (J) from which we obtained the reconstructed graph (K). L. A directionality analysis applied to the graph recovered from step 2 allowing to recover an oriented graph of the ER network.  Table S3: Estimation of the potential well characteristics as a function of the energy, related to STAR Methods. Based on a circular well (a = b = 150 nm) and using hybrid MLE algorithm and an acquisition time ∆t = 50 ms.