A maximum mean discrepancy approach reveals subtle changes in α-synuclein dynamics

Numerous models have been developed to account for the complex properties of the random walks of biomolecules. However, when analysing experimental data, conditions are rarely met to ensure model identification. The dynamics may simultaneously be influenced by spatial and temporal heterogeneities of the environment, out-of-equilibrium fluxes and conformal changes of the tracked molecules. Recorded trajectories are often too short to reliably discern such multi-scale dynamics, which precludes unambiguous assessment of the type of random walk and its parameters. Furthermore, the motion of biomolecules may not be well described by a single, canonical random walk model. Here, we develop a methodology for comparing biomolecule dynamics observed in different experimental conditions without beforehand identifying the model generating the recorded random walks. We introduce a two-step statistical testing scheme. We first use simulation-based inference to train a graph neural network to learn a fixed-length latent representation of recorded random walks. As a second step, we use a maximum mean discrepancy statistical test on the vectors of learnt features to compare biological conditions. This procedure allows us to characterise sets of random walks regardless of their generating models. We initially tested our approach on numerical trajectories. We then demonstrated its ability to detect changes in α-synuclein dynamics at synapses in cultured cortical neurons in response to membrane depolarisation. Using our methodology, we identify the domains in the latent space where the variations between conditions are the most significant, which provides a way of interpreting the detected differences in terms of single trajectory characteristics. Our data show that changes in α-synuclein dynamics between the chosen conditions are largely driven by increased protein mobility in the depolarised state. Author summary The continuous refinement of methods for single molecule tracking in live cells advance our understanding of how biomolecules move inside cells. Analysing the trajectories of single molecules is complicated by their highly erratic and noisy nature and thus requires the use of statistical models of their motion. However, it is often not possible to unambiguously determine a model from a set of short and noisy trajectories. Furthermore, the heterogeneous nature of the cellular environment means that the molecules’ motion is often not properly described by a single model. In this paper we develop a new statistical testing scheme to detect changes in biomolecule dynamics within organelles without needing to identify a model of their motion. We train a graph neural network on large-scale simulations of random walks to learn a latent representation that captures relevant physical properties of a trajectory. We use a kernel-based statistical test within that latent space to compare the properties of two sets of trajectories recorded under different biological conditions. We apply our approach to detect differences in the dynamics of α-synuclein, a presynaptic protein, in axons and boutons during synaptic stimulation. This represents an important step towards automated single-molecule-based read-out of pharmacological action.

Numerical models make it possible to generate synthetic observations of biological 2 systems across a broad range of parameters. However, the computational cost of 3 directly using these simulations to perform statistical inference is often prohibitive [1]. 4 The reason is that both the likelihood and evidence are intractable in most systems, and 5 these inferences must therefore be addressed as likelihood-free simulation-based 6 inferences [2]. The primary approach to simulation-based inference is approximate 7 Bayesian computation (ABC) [2]. ABC relies on comparing user-defined summary 8 statistics from experimentally recorded and simulated data using a chosen distance 9 metric. The ever-growing amount of available data, along with recent advances in deep 10 learning [3] make it possible to capture more and more detailed properties of 11 experimental systems, and have thus boosted the development of simulation-based 12 inference. We refer the interested reader to [1] for a detailed taxonomy of 13 simulation-based inference methods. 14 Typical inference schemes developed to analyse biomolecule trajectories focus on 15 estimating physical parameters such as the diffusion coefficient, the anomalous diffusion 16 exponent, the type of random walk model, or other ad-hoc quantities measuring 17 particular aspects of the dynamics. Here, instead of describing trajectories using a set of 18 explicitly defined features, we rely on an encoder neural network, in order to 19 characterise each trajectory by a latent vector of features. The goal of this encoder is to 20 automatically learn optimised features that describe random walks beyond predefined 21 canonical models and features. We employ the encoder network to project recorded 22 trajectories into the fixed-dimensional latent space. We develop a statistical test on this 23 latent space to test for differences in dynamics between two sets of trajectories. Our 24 methodology can in particular be used to compare dynamics observed in different 25 biological conditions and different cell organelles, by comparing the sets of latent vectors 26 computed from trajectories observed in the respective microscopy recordings or regions 27 of the cell. The central advantage of the testing methodology we propose is that it is 28 not dependent on the specification and selection of a model of the recorded random 29 walks. This enables statistically robust testing of differing biological conditions, which 30 are likely to induce different levels of cellular heterogeneity and do not necessarily 31 generate canonical random walks. 32 We train the encoder network using a simulation-based inference framework 33 (detailed below), allowing it to provide a representation of trajectories without assuming 34 that they are realisations of a canonical random walk model. The subsequent statistical 35 test seeks to differentiate the distributions of latent vectors coming from different 36 conditions or organelles (or both). It is based on the maximum mean discrepancy 37 (MMD) test [4], which uses a kernel approach to compare two distributions. This test 38 allows us to compare sets composed of different numbers of trajectories, and provides a 39 means for interpreting the differences between biological conditions. Finally, we show 40 the robustness of the approach to the intrinsic variability of biological observations, and 41 demonstrate that the statistical differences do not stem from a single experiment, nor 42 from an outlier composed of a minority of trajectories. 43 We demonstrate our methodology by studying the dynamics of α-synuclein inside 44 and outside of synapses. α-synuclein is a small, soluble, and highly mobile protein (140 45 amino acid residues) that is strongly accumulated in presynaptic boutons (reviewed 46 in [5]). Experiments based on fluorescence recovery after photobleaching [6] have shown 47 the existence of at least two main modes of diffusion, one in which α-synuclein is 48 transiently bound to synaptic vesicles in the synaptic bouton, and another in which the 49 protein diffuses freely both in axons and in synaptic regions. The existence of an 50 immobile population of α-synuclein molecules, taking the form of protein aggregates at 51 synapses, has also been proposed [6]. In response to strong depolarising signals the 52 bound population of α-synuclein dissociates from its synaptic binding sites and 53 disperses in the neighbouring axon [7]. In agreement with these earlier studies, we found 54 that α-synuclein dynamics differ between synapses and axons. Furthermore, Primary murine cortical neuron cultures were prepared at embryonic day E17 as 66 described previously [8]. Cortices were dissected, the tissue was dissociated and the cells 67 were seeded at a concentration of 5 x 10 4 cm −2 on glass coverslips that had been coated 68 with poly-D,L-ornithine. Neurons were kept at 37 o C and 5% CO 2 in neurobasal 69 medium supplemented with Glutamax, antibiotics and B27 (all from Gibco, Thermo [mM]: 120 NaCl, 2.5 KCl, 2 CaCl 2 , 2 MgCl 2 , 25 glucose, 5 pyruvate, 25 HEPES,80 adjusted to pH 7.4) at room temperature, using an inverted Nikon Eclipse Ti microscope 81 equipped with a 100x oil objective (NA 1.49), an Andor iXon EMCCD camera (16 bit, 82 image pixel size 160 nm), and Nikon NIS acquisition software. First, an image of the 83 chosen field of view (average of 10 image frames taken with 100 ms exposure time) was 84 taken in the green channel (non-converted mEos4b fluorescence) using a mercury lamp 85 A reference image (left panel) was taken in the green channel, followed by a SMLM movie of 25000 frames in the red channel (panel two). (B) A second recording of the same field of view (image and movie) were then acquired in the presence of elevated KCl concentration. Note the dispersal of αSyn:Eos4 in response to depolarisation compared to the control (left panels). (C) A third image and movie were acquired after addition of 2% paraformaldehyde. The third column of images shows zoomed SMLM reconstructions of the synaptic bouton indicated in the first image. The fourth column depicts a subset of trajectories from the same synaptic terminal. and specific excitation (485/20 nm) and emission filters (525/30 nm). This was followed 86 by a streamed acquisition of 25 000 movie frames recorded with 15 ms exposure and 87 ∆t = 15.4 ms time lapse (total duration: 6 min 25 s) in the red channel using a 561 nm 88 laser at a nominal power of 150 mW for excitation (inclined illumination), together with 89 pulsed activation lasers applied during the off time of the camera (405 nm, approx.  Finally, the neurons were fixed with the addition of phosphate buffer at pH 7.4 100 containing 4% paraformaldehyde and 1% sucrose (final concentration 2% PFA), and a 101 third reference image (green) and SMLM movie (red channel) were acquired in the presence of the fixative (Fig. 1C).

103
Image processing and analysis 104 SMLM image stacks (tiff files) were pre-processed in order to remove background 105 fluorescence using a quantile filter computed on a sliding window. Then, localisations 106 were detected using the algorithm described in [9], based on a wavelet analysis. Subpixel 107 localisation was performed using the radial symmetry center algorithm introduced 108 in [10]. Sample drift was corrected by subtracting the displacement that yielded the 109 best correlation between densities of successive temporal slices grouping 10 000 110 localisations each. To isolate axons, we applied a Sato filter [11] with a width of 3 pixels 111 on the logarithm of the pixel-wise mean intensity. Then, we used a local thresholding 112 algorithm, provided by [12] to compute a mask over the image. All steps of the analysis 113 were implemented in Python, the code is available at http://gitlab.pasteur.fr.

114
Synapses were manually detoured using an ad-hoc graphic user interface. In total, our 115 analysis includes 321 synapses for which more than 150 trajectories were recorded, 116 coming from 10 different fields of view. A synapse is counted twice if it appears in two 117 or three recordings based on the same field of view but done in different conditions (e.g. 118 some synapses appear in control, KCl and fixed conditions).

119
In the analysis of experimental trajectories, we considered only trajectories located 120 in the axons, and we split them into two groups: those located outside the synaptic 121 region and those located inside. Synaptic regions were delimited by a density threshold 122 of one tenth of the maximum density of detections in the synapse. The density was 123 estimated using a Gaussian kernel method with a bandwidth of 150 nm. We estimated 124 the apparent effective diffusivity [13] of a trajectory from the sample variance of its  The first step of the analysis is to build a latent representation of random walks that 130 does not require strong assumptions about the underlying generative models. In this 131 section, we present the simulation-based inference scheme, the architecture of the neural 132 network used to compute latent vectors from trajectories, and the visualisation of these 133 vectors.

134
Simulation-based inference 135 In order to ensure that our characterisation of trajectories is accurate, robust and 136 length-independent, it should be trained on an as wide as possible variety of random 137 walks. Hence, we chose to rely on a simulation-based inference procedure [1]. It consists 138 in generating data on which the neural network is subsequently trained. In our case, 139 this amounts to simulating trajectories of a variety of models known to encapsulate 140 different properties of biomolecule dynamics in cells. The physical parameters chosen to 141 simulate these trajectories should at least cover the range of the experimentally 142 observed ones, in order to ensure that the network is able to encode relevant 143 information about the recorded trajectories on which the inference will eventually be 144 performed after its training.

145
To ensure the diversity of the training set, we simulated trajectories using five 146 different canonical random walk models covering a wide spectrum of possible random 147 walk characteristics:

158
The models' parameters were drawn from the same distributions throughout the 159 entire study and were chosen to cover the entire ranges observed experimentally. Var(log 10 (D/D 0 )) = 0.5 2 . For the OU model, the relaxation rate θ was drawn from ... 164 We added uncorrelated localisation noise to each point of the trajectories, drawn 165 from a centred Gaussian distribution with standard deviation drawn uniformly between 166 15 and 40 nm (to include the signal intensity dependence of the localisation precision). 167 The time lapse between recordings was set equal to that of the camera (15,4 ms).

168
The neural network (the architecture of which is detailed below) was then trained to 169 infer two characteristics of interest from the trajectories: their anomalous diffusion 170 exponent (if applicable), and the random walk model from which they were generated 171 among the five described above. Throughout the training, the network processed ∼ 10 6 172 independent simulated trajectories.

173
Graph neural network and random walks 174 Here, we explain how we construct an informative size-independent representation of 175 random walks. We feed trajectories into a neural network, which we train to compute a 176 vector of summary statistics from each trajectory. It contains information relevant to 177 the characterisation of the random walk. Although observed trajectories are not all of 178 the same length, the summary statistics is a vector of constant size. Thus, this vector 179 can be used to compare trajectories of different sizes. Details about the subsequent 180 analyses are found in the next subsections. Here, we focus on the neural network 181 architecture, which is based on our previous work [24]. We refer the interested reader to 182 the Supporting Information for implementation details 183 Graphical models are methods of choice to handle complex inferences [2,25], model 184 large scale causal relationships [26] and provide inductive biases in Bayesian 2 ∈ Y, of size n y , contains information about the trajectory's course between pairs of 201 nodes. The graph construction is illustrated in Fig. 2A.

202
A trajectory's graph and its associated feature vectors are then passed to a series of 203 graph convolution layers, as illustrated in Fig. 2B. Node features vectors are aggregated 204 in the pooling layer such that in the second part of the network, each trajectory is 205 represented by a fixed-size vector. The output of the multi-layer perceptron directly 206 downstream of the pooling layer provides a 16-dimensional vector representing each 207 trajectory, i.e. a summary statistics, which we designate in the following as the "latent 208 representation", or "latent vector". During the training phase, the latent vector is fed to 209 two separate MLPs predicting the anomalous exponent and the underlying model of the 210 random walk, respectively. The loss minimised during the network's training is the sum 211 of two task-specific losses: the mean squared error of the prediction of the anomalous Once processed by the encoder, each trajectory is reduced to a 16-dimensional vector.

218
For visualisation purposes, we projected this vector on a 2D plane using a 219 parametric-UMAP to perform the projection [32]. This variation of UMAP allows us to 220 learn the transformation projecting the data from 16 to 2 dimensions solely on 221 simulated trajectories, so that it is independent of the experimental trajectories and is 222 only trained once. The GNN was trained first, then the parametric-UMAP projection 223 was learnt and both their sets of weights were frozen. 224 We designate as "2D latent representation" the two-dimensional vector, output by 225 the parametric-UMAP, which represents a trajectory in the latent space. Figure 3A   226 shows that latent representations of simulated and experimental data largely overlap, 227 and that the experimental trajectories fall within the region covered by the simulated 228 ones. Figures 3B and 3C show that the random walk model and the diffusivity are 229 prominent determinants of the latent space structure. Figures 3D, 3E and 3F illustrate 230 the diversity of αSyn:Eos4 trajectories that can be found in a presynaptic bouton, and 231 how this diversity is captured by the latent representation. Figure 3F highlights the fact 232 that there is a high variability of trajectory dynamics even within a given synapse, 233 suggesting that α-synuclein molecules can transition between various dynamic modes.

234
Using the approach described above, we can associate to any set of trajectories, a set 235 of constant-sized vectors characterising their dynamics. Each microscope recording, or 236 organelle within it, can thus be characterised by a set of N 2-dimensional feature The graph is passed through a series of graph convolution layers (shown in blue), which propagate information along edges. The pooling operation (green) combines all node feature vectors from a graph into a vector of fixed size representing the graph. This vector is then passed to a multi-layer perceptron (MLP in orange), whose output we refer to as the "latent representation" of the trajectory. The latent representation is fed to two task-specific MLPs: one that predicts the trajectory's anomalous exponent α and one that assigns a vector of probabilities for the trajectory to have been generated by each of the models considered. organelle is characterised by the set of its latent vectors. Thus, we base our statistical 245 test on the comparison of the generating distributions of these vectors. In the absence 246 of a priori knowledge of these distributions, we employ a kernel-based approach: the 247 maximum mean discrepancy (MMD) test.
where E x and E y denote expectation w.r.t. p and q, respectively.

257
If the function class is the unit ball in a Reproducing Kernel Hilbert Space 258 (RKHS) [33] H, the square of the MMD can directly be estimated from data samples.

259
Denoting k the kernel operator such that ∀f ∈ F, f (x) = f, k(x, ·) , an unbiased 260 estimator of the square of the MMD between X and Y is given by: (2) In our case, X = R 2 , and we used the classical Gaussian kernel  The MMD is capable of detecting subtle differences such as the ones between data 266 generated by generative adversarial networks (GANs) and real data [34]. It has also 267 proved efficient in discovering which variables exhibit the greatest difference between 268 datasets [4,33].

269
Statistical test 270 We adapted the bootstrap test described in [4] to assess whether dynamics of greater than the 1 − α quantile of this distribution, we rejected H 0 . This test is said to 287 be of level α, because with probability α, we will reject the null hypothesis when it is 288 actually true.

290
Detecting differences between sets of simulated trajectories 291 To assess the performance of the full statistical testing framework, we applied it on 292 simulated data. We set the level of statistical significance to α = 0.05, and we simulated 293 trajectories as described in Material and Methods.

294
A first case of our test is to detect changes in the proportions of given types of 295 trajectories between two sets of observations. This is illustrated in Fig. 4, where we 296 show example comparisons in the 2D latent space between two sets with different 297 proportions of their trajectories generated by fBM and sBM. We compared fBM and 298 sBM, since they share numerous features and because for a large range of values of the 299 anomalous exponent they are challenging to distinguish [24]. Furthermore, these two 300 random walk models are highly representative of our experimental data, as can be seen 301 by their latent space occupations (compare Figs. 3A and B).

302
The difficulty of separating the two populations depends on their relative 303 proportions in the two datasets, and we see that both the amplitude of the witness 304 function (Fig. 4C) as well as the value of the test statistic (Fig. 4D) decrease as the 305 ratio is closer to 1:1. When the two sets are drawn from the same 50/50 distribution, 306 the test does not, and should not, find significant differences between them.

307
The other main factor determining the difficulty of detecting a difference, is the size 308 of the datasets. Experimentally, changes in biological conditions lead not only to 309 changes in the properties of the random walks. It also leads to changes in the total 310 number of trajectories of a given type. This causes challenges in performing proper 311 statistical testing. To quantitatively assess the effect of both the number of trajectories 312 and the relative proportions belong to different random walk classes, we conducted 313 numerical experiments where we varied these two parameters systematically (Fig. S2A). 314 Besides differences in the proportions of trajectories generated by different random 315 walk models, the sets may also differ in the models' parameter values. We thus 316 additionally evaluated the test's ability to distinguish two sets of fBMs, one with 317 anomalous diffusion exponent α = 1 − δ and the other with α = 1 + δ. Our results 318 indicate that in both cases 1 000 trajectories are sufficient to detect subtle changes 319 between distributions (Fig. S2). Fewer trajectories are needed to detect starker 320 differences. In cases where the compared sets are drawn from the same distribution 321 (ν = 0 and δ = 0), the null hypothesis is rejected in about 5% of cases, consistent with 322 our chosen α-level.

323
One way to further improve these results is to optimise the kernel used to compute 324 the MMD. We show in Fig. S1 how kernel bandwidth and shape affect the power of the 325 test. If the kernel bandwidth is too small, this weakens the test by making it too 326 sensitive to noise. Conversely, if its bandwidth is too large, this prevents the test from 327 detecting subtle changes. We tested the effect of the kernel characteristics in the same 328 setting as illustrated in Figs. 4 and S2A, with N = 200 trajectories in each set and 329 comparing sets with 70% fBM / 30% sBM and 30% fBM / 70% sBM. We observed that 330 a Gaussian kernel with radius σ equal to the median pairwise distance in the dataset 331 (i.e. σ between 1.5 and 2) yields a near-optimal test, in agreement with earlier 332 findings [4]. Finally, while we have here focused on optimizing type II error while 333 controlling type I error (i.e. fixed α-level), the parameters and the functional form of the 334    In all of the three cases illustrated in Fig. 5, the empirical MMD 2 u is significantly 352 higher than the top 1% quantile of the null distribution, meaning that our test detects a 353 significant difference in the properties of the trajectories of the two compared subsets at 354 the 1% α-level. Note the wide range of magnitudes spanned by these differences, which 355 is well judged by looking at the absolute intensity of the witness functions shown in 356 Fig. 5E. We see that the difference induced by fixation on α-synuclein mobility at 357 synapses (Fig. 5E, middle) is much more pronounced than the one induced by high KCl 358 April 8, 2022 13/23 treatment (Fig. 5E, top). This is also apparent when looking at how large MMD 2 u is 359 compared to its distribution under the null hypothesis: in the case of the fixed vs 360 control comparison (Fig. 5F, middle), the histogram of the MMD 2 u values obtained 361 under the null hypothesis is completely squeezed because the empirical MMD 2 u is more 362 than two orders of magnitude larger than the top 1% quantile of the distribution under 363 the null hypothesis. In comparison, KCl treatment produces a less drastic change in 364 α-synuclein trajectories located in synaptic terminals (Fig. 5F, top), although this effect 365 is also highly statistically significant (p 0.01). 366 We further observed that, while their magnitudes differ, the witness functions of the 367 control/KCl comparisons in axons and in synaptic boutons (Fig. 5E, top and bottom)   Furthermore, as illustrated in Fig. S3, we use this statistic to check that all acquired 389 fields of view contribute evenly to the difference between the control and KCl conditions. 390 We looked at the proportion of trajectories coming from each field of view and condition 391 in a region of the latent space which we define as "critical", based on the value of S (it 392 is the contiguous domain containing the maximum of S and where S is greater than half 393 of its maximum value). We could thus confirm that, on the one hand, trajectories 394 located in this region of the latent space originate from all considered fields of view in a 395 balanced manner, and on the other hand, that within each field of view, the difference 396 of representation of each condition is in the same direction (except for one field of view, 397 where there is almost no difference). This furthermore excludes the possibility that the 398 observed differences are solely due to a single abnormal recording.

399
Comparing synapses 400 The approach we propose is not restricted to inter-condition comparisons, but can be 401 used to compare any two subsets of trajectories. Hence, we can group trajectories by 402 synapse and external condition (control, KCl or fixed), and compute the value of 403 MMD 2 u of all the pairs of so-obtained synapses. This provides us with an inter-synapse 404 distance matrix, shown in Fig. 7A. Using these distances, we can embed the trajectory 405 subsets in an Euclidian space, i.e. summarise each subset by a vector of fixed dimension, 406 using for instance the multi-dimensional scaling (MDS) algorithm [2]. We adapted the 407 MDS algorithm in order to account for the uncertainty that we have in the estimation 408 of the squared distance, which notably depends on the number of observed trajectories 409 Our approach allows extracting physically and biologically relevant results without 425 having to assign the biomolecule motion to canonical models. Although these models 426 are instrumental for interpreting the properties of random walks, the complexity and 427 heterogeneity of biological environments at the nanometer/micrometer scale often 428 precludes an unambiguous model assignation. 429 An alternative approach that would also use a set of canonical random walks as a 430 possible basis for description is Bayesian Model Averaging (BMA). The BMA method 431 does not look for the best model to describe experimental data but rather evaluates the 432 parameters for each of the possible models and then averages the different results. It is 433 challenging to apply this approach to our current problem for two reasons: (i) the space 434 parameters of the different random walks are not identical and (ii) the evaluation of the 435 parameters associated with each random walk by Bayesian methods is not possible for 436 all models. It would be possible to develop Bayesian method variations to estimate the 437 parameters of models that do not have tractable likelihoods. However, the 438 marginalisation of the different models would be computationally intensive. Finally, 439 even though BMA could include information from multiple models, it would still rely on 440 the projection of the experimental date onto a set of canonical models. 441 We applied our approach to study the dynamics of α-synuclein molecules in axons 442 and presynaptic boutons. In agreement with earlier studies of the population dynamics 443 of α-synuclein [6,7], we found that the protein assumes differing dynamic states at   Beyond the automation of the analysis procedure between biological conditions, our 464 approach is well suited for exploratory data analysis. The capacity to project individual, 465 differently sized trajectories into finite sized vectors makes it possible to study precise 466 sub-cellular compartments or organelles in a standardised form, and thus allows to test 467 statistical differences between these regions. Hence, recorded single molecule data can 468 be searched in order to detect and characterise regions of the cell that have different 469 statistical properties. This exploration can be done even in regions with different 470 trajectory densities, as is the case for α-synuclein at synapses versus axonal domains.

471
One of the current limitations of the current approach is the difficulty in evaluating 472 the type II error [35] bounds on the statistical test. The MMD test is applied within the 473 latent space of the GNN. This manifold is built using a set of non-linear operations, 474 which depends both on the numerical trajectories seen during the training and on the 475 cost function being optimised. Hence, there may be domains within the latent space 476 that could lead to improper sensitivity of the statistical test. As can be seen in [24] and 477 in Fig. 3, different types of random walks occupy domains of different size and there is a 478 large overlap of the regions. Since our approach relies on a simulation based framework, 479 it is possible to use numerical simulations matching the experimental occupancy of the 480 latent space to evaluate the accuracy of the test. Furthermore, extensive simulations The funding sources had no role in study design, data collection and analysis, 497 decision to publish, or preparation of the manuscript.

498
Conflicts of interest. Hippolyte Verdier and Alhassan Cassé are Sanofi employees 499 and may hold shares and/or stock options in the company. The other authors declare to 500 have no financial or non-financial conflicts of interest. 501 the uncertainty, which can be evaluated by bootstrapping, is sometimes of the same 611 order of magnitude than the estimated value. Furthermore, the number of elements per 612 set (here, the number of trajectories per synapse), spans more than an order of 613 magnitude and uncertainty thus greatly varies from one measure to the other. This 614 should be taken into account when using MMD 2 u [F, X, Y ] to embed subsets of 615 trajectories.

616
Hence, starting from a matrix of squared distances D 2 between N sets of trajectories, and a matrix of uncertainties of these squared distances D 2 σ , we obtain a set of N Euclidian vectors {x 1 , . . . , x N } by maximizing the probability of the resulting squared distances, assuming that they follow Gaussian laws whose means are the coefficients of D 2 and standard deviations coefficients of D 2 σ . This amounts to solving the following optimisation problem : max x1,...,x N i<j We do so using a gradient ascent method. Prior to entering the encoder, each trajectory is turned into a graph as described in . To 620 each node is associated a vector containing the following features :    Origin of trajectories found in the critical region. These counts were obtained using a set composed of the same number n = 1, 000 of trajectories from each microscopy recording. We randomly subsampled those who had more intra-synaptic trajectories, and discarded those who had less than 1,000 (hence the column with missing number of KCl trajectories).