Neuron cultures and αSyn:Eos4 expression.

Primary murine cortical neuron cultures were prepared at embryonic day E17 as described previously [25]. Cortices were dissected, the tissue was dissociated and the cells were seeded at a concentration of 5 × 10⁴ cm⁻² on glass coverslips that had been coated with poly-D,L-ornithine. Neurons were kept at 37°C and 5% CO₂ in neurobasal medium supplemented with Glutamax, antibiotics and B27 (all from Gibco, Thermo Fisher Scientific), infected at day in vitro (DIV) 11–24 with lentivirus driving the expression of α-synuclein tagged at its C-terminus with the photoconvertible fluorescent protein mEos4b (αSyn:Eos4) under the control of a ubiquitin promotor, and used for experiments 7 days later. All cell culture and imaging experiments were conducted at the Laboratory for cellular synapse biology at IBENS (Paris). Procedures involving animals were performed according to the guidelines set out by the local veterinary and administrative authorities.

Single molecule localisation microscopy (SMLM).

Living neurons expressing αSyn:Eos4 were imaged in modified Tyrode’s solution (in [mM]: 120 NaCl, 2.5 KCl, 2 CaCl₂, 2 MgCl₂, 25 glucose, 5 pyruvate, 25 HEPES, adjusted to pH 7.4) at room temperature, using an inverted Nikon Eclipse Ti microscope equipped with a 100x oil objective (NA 1.49), an Andor iXon EMCCD camera (16 bit, image pixel size 160 nm), and Nikon NIS acquisition software. First, an image of the chosen field of view (average of 10 image frames taken with 100 ms exposure time) was taken in the green channel (non-converted mEos4b fluorescence) using a mercury lamp and specific excitation (485/20 nm) and emission filters (525/30 nm). This was followed by a streamed acquisition of 25 000 movie frames recorded with 15 ms exposure and Δt = 15.4 ms time lapse (total duration: 6 min 25 s) in the red channel using a 561 nm laser at a nominal power of 150 mW for excitation (inclined illumination), together with pulsed activation lasers applied during the off time of the camera (405 nm, approx. 1–5 mW; 488 nm, 10 mW; 0.45 ms pulse). The red emission of the photo-converted mEos4b fluorophores was detected with a 607/36 nm filter (Fig 1A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Single molecule localisation microscopy (SMLM) of α-synuclein in cortical neurons.

(A) Neurons expressing αSyn:Eos4 were first imaged in control condition. 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.

https://doi.org/10.1371/journal.pcbi.1010088.g001

After recording of the baseline dynamics of αSyn:Eos4, the buffer composition was changed with the addition of Tyrode’s solution containing elevated KCl at the expense of NaCl (final concentrations in [mM]: 78 NaCl, 44.5 KCl, 2 CaCl₂, 2 MgCl₂, 25 glucose, 5 pyruvate, 25 HEPES, pH 7.4). This treatment causes the depolarisation of the neurons leading to the dissociation of α-synuclein from its binding sites in the synaptic bouton [24]. A reference image was taken in the green channel, followed by a second SMLM recording (Fig 1B), starting approximately 7.5 min after the first acquisition. Finally, the neurons were fixed with the addition of phosphate buffer at pH 7.4 containing 4% paraformaldehyde and 1% sucrose (final concentration 2% PFA), and a third reference image (green) and SMLM movie (red channel) were acquired in the presence of the fixative (Fig 1C).

Image processing and analysis.

SMLM image stacks (tiff files) were pre-processed in order to remove background fluorescence using a quantile filter computed on a sliding window. Then, localisations were detected using the algorithm described in [26], based on a wavelet analysis. Subpixel localisation was performed using the radial symmetry center algorithm introduced in [27]. Sample drift was corrected by subtracting the displacement that yielded the best correlation between densities of successive temporal slices grouping 10 000 localisations each. To isolate axons, we applied a Sato filter [28] with a width of 3 pixels on the logarithm of the pixel-wise mean intensity. Then, we used a local thresholding algorithm, provided by [29] to compute a mask over the image. All steps of the processing were implemented in Python using the palmari package. Synapses were manually detoured using an ad-hoc graphic user interface. In total, our analysis includes 321 synapses for which more than 150 trajectories were recorded, coming from 10 different fields of view. A synapse appearing in several recordings based on the same field of view but done in different conditions (e.g. some synapses appear in control, KCl and fixed conditions) is counted several times.

In the analysis of experimental trajectories, we considered only trajectories located in the axons, and we split them into two groups: those located outside the synaptic region and those located inside. Synaptic regions were delimited by a density threshold of one tenth of the maximum density of detections in the synapse. The density was estimated using a Gaussian kernel method with a bandwidth of 150 nm. We estimated the apparent effective diffusivity [30] of a trajectory from the sample variance of its single-time-lapse displacements, i.e. , where Δr_i = r_i − r_i−1 is the displacement between the (i−1)th and ith recorded positions and is the average displacement.

Simulation-based inference.

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

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

• the Levy walk (LW) [31–33], which has non-Gaussian increments and exhibits weak ergodicity breaking;
• scaled Brownian motion (sBM) [34–36], which is Gaussian, non-stationary and weakly non-ergodic;
• the Ornstein Uhlenbeck process (OU) [37], a Gaussian, stationary process with exponentially decaying autocorrelations;
• fractional Brownian motion (fBM) [38], which is Gaussian, stationary and exhibits slowly decaying temporal correlations;
• and the continuous time random walk (CTRW) [10, 39], which is non-Gaussian, shows weak ergodicity breaking, ageing, and has discontinuous paths.

The models’ parameters were drawn from the same distributions throughout the entire study and were chosen to cover the entire ranges observed experimentally. Trajectory lengths were drawn from a log-uniform distribution between 7 and 25 points, which corresponds to a mean length of 14 points. The effective diffusivity was drawn from a log-normal distribution with D₀ = 1μm²/s, 〈log₁₀(D/D₀)〉 = −0.5 and Var(log₁₀(D/D₀)) = 0.5². For the OU model, the relaxation rate θ was log-uniformly drawn between 0.01 and 1.

We added uncorrelated localisation noise to each point of the trajectories, drawn from a centred Gaussian distribution with standard deviation drawn uniformly between 15 and 40 nm (to include the signal intensity dependence of the localisation precision). The time lapse between recordings was set equal to that of the camera (15,4 ms).

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

Graph neural network and random walks.

Here, we detail the architecture of the neural network which we use to compute summary statistics vectors describing trajectories. Named GRATIN (for “Graphs on trajectories for inference”), it is inspired from the architecture presented in our previous works [40, 41]. An important feature of this encoder network is that the size of its output is independent of the trajectory length, facilitating the comparison of trajectories of different sizes. Details about the subsequent analyses are found in the next subsections, and we refer the interested reader to the Supporting Information for implementation details. The code for the encoder network (architecture and training) is available at https://github.com/hippover/gratin.

Graphical models are methods of choice to handle complex inferences [15, 42], model large scale causal relationships [43] and provide inductive biases in Bayesian inferences [44]. Over the last five years, graph-based analysis methods have been complemented by graph neural networks (GNNs), which extend classical neural network approaches to learn and process relations between data points in a flexible manner, using graph-convolutions operations [45]. GNNs are efficient at representation learning [46–48] and naturally apply to data of varying dimensions. Furthermore, GNNs are ideally suited to encode natural symmetries of physical systems [49]. For these reasons, we chose to use a such network to process trajectories.

To do so, we represent a trajectory R = (r₁, r₂, …, r_N) by a directed graph G = (V, E, X, Y). Here V = {1, 2, …, N} are the nodes, each associated with a position in the trajectory, E ⊆ {(i, j)|(i, j) ∈ V²} is the set of edges connecting pairs of nodes, are node feature vectors, and are edge features. Each node feature vector , of size n_x = 6, captures information associated to the course of the trajectory until point i: normalised time (between 0 and 1), normalised sums of step sizes to the powers 1, 2 and 4, up to point i, normalised distance to origin and maximal distance to origin up to point i. The normalization terms were chosen so that features do not directly account for the trajectory’s scale (a global scale is computed instead and intervenes downstream of graph convolutions). The edges incoming at each node originate only from nodes in the past (respecting causality): node i receives edges from nodes i − ⌊γ⁰⌋, i − ⌊γ¹⌋, …, i − ⌊γ^k−1⌋, with γ = (i − 1)^1/(k−1) and where k is the fixed maximal in-degree of nodes. We remove eventual repeated edges. For trajectories shorter than k, we connected each node to all its predecessors. In the following, we chose k = 10.

Each edge feature vector , of size n_y = 6, contains information about the trajectory’s course between the two nodes it connects: normalised time difference, normalised distance, normalised sums of step sizes to the powers 1, 2 and 4, between points i and j, correlation of the two displacements. The graph construction is illustrated in Fig 2A and exact expressions of edge and node features are provided in S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Model architecture.

(A) Construction of a graph from a single trajectory (left). Positions, colored according to time, are treated as nodes for which features X are computed (middle). Nodes are then connected by edges (grey lines) following a given wiring scheme, and edge features Y are computed. Edges terminating at the trajectory’s last point are shown in red. Feature matrices for both nodes (X) and edges (Y) are depicted in small insets with color coded values. (B) Graph neural network. 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.

https://doi.org/10.1371/journal.pcbi.1010088.g002

After initialization, the nodes and edge features are embedded into a space of higher dimension (16 for edges, 10 for nodes) using a multi-layer perceptron. Then, nodes vectors are sequentially updated by a series of graph convolution layers, as illustrated in Fig 2B. More precisely, we used three successive GIN convolution layers introduced in [50], with 32 convolution filters each. A pooling layer follows these convolutions, aggregating all node features vectors of a trajectory using an attention mechanism, such that in the second part of the network, each trajectory is represented by a single vector whose dimension is independent of the trajectory length. A feature accounting for the scale of the trajectory is concatenated to this vector, which is then passed to a multi-layer perceptron. Its output is a 16-dimensional vector of learnt summary statistics, which we designate in the following as the “latent representation”, or “latent vector” obtained from GRATIN. We purposely chose to use a relatively high number of latent dimensions to benefit from the effects of over-parameterization [51]. During the training phase, the latent vector is fed to two separate MLPs predicting the anomalous exponent and the underlying model of the random walk, respectively. The loss minimised during the network’s training is the sum of two task-specific losses: the mean squared error of the prediction of the anomalous exponent (excluding OU trajectories which are not anomalous random walks) and the cross-entropy of the predicted and true model classes. Training the network on such physically informed tasks makes it build a relevant latent representation of trajectories, as shown in our previous work [40]—the architecture presented here is an improved version of the previously published one.

Latent representation of trajectories.

As mentioned above, once processed by the encoder, each trajectory is projected into a 16-dimensional space. Maintaining a relatively high dimension for the latent space helps the network’s convergence. Nonetheless, such high-dimensional statistics vectors are not convenient for visual interpretation of the subsequent results. We thus chose, in the rest of the process, to resort to a compressed 2D version of this latent vector, projecting it on a plane using parametric-UMAP [52]. This variation of UMAP allows us to learn the transformation projecting the data from 16 to 2 dimensions on simulated trajectories, so that it is independent of the experimental trajectories and is only trained once. The GNN was trained first, then the parametric-UMAP projection was learnt and both their sets of weights were frozen. We designate as “GRATIN 2D latent representation” the two-dimensional vector, output by the parametric-UMAP for each trajectory. We compared the power of MMD tests based on the 16- and 2-dimensional latent representations of trajectories and found no significant difference of performance (see S4 Fig and S1 Text).

By comparing Fig 3A, 3B and 3C with respectively Fig 3D, 3E and 3F, we first see that latent representations of simulated and experimental data largely overlap, and that the experimental trajectories fall within the region covered by the simulated ones. We furthermore observe that the random walk model, the diffusivity and the anomalous diffusion exponent are prominent determinants of the latent space structure. Finally, the fact that parameter values predicted on experimental trajectories match the true values of simulated trajectories lying in the same region of the latent space demonstrates the robustness of the encoding. Fig 3G, 3H and 3I illustrate the diversity of αSyn:Eos4 trajectories that can be found in a presynaptic bouton, and how this diversity is captured by the latent representation: trajectories whose latent vectors are located in remote regions of the latent space exhibit distinct dynamics. This diversity of dynamics within a given synapse suggests that α-synuclein molecules can transition between various dynamic modes and that its state is not exclusively determined by its location.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. 2D Latent space of trajectories.

(A, B, C) Latent representations of simulated trajectories, colored according to the true log-diffusivity, the true random walk model (fBM: fractional Brownian motion; LW: Levy walk; sBM: scaled Brownian motion; OU: Ornstein-Uhlenbeck process; CTRW: continuous-time random walk) and the true anomalous diffusion exponent. (D, E, F) Latent representations of recorded trajectories, colored according to the same variables as above, this time using inferred values. (G) Recorded trajectories at synapses and in the axon. (G) Zoom on a presynaptic bouton, delimited by the red square in panel D. Three individual trajectories are highlighted. (I) Latent representation of the trajectories at this synapse, with colored dots corresponding to the three trajectories highlighted in H. (J) Examples of acquired trajectories, located according to their position in the latent space. Each square has a side length of 1.5 micrometer, trajectories length vary from 14 to 16.

https://doi.org/10.1371/journal.pcbi.1010088.g003

Using the approach described in this subsection, we can associate to any set of trajectories, a set of constant-sized vectors characterising their dynamics. Each microscope recording, or organelle within it, can thus be characterised by a set of N 2-dimensional feature vectors, N being the number of trajectories. The conditions are thus met to perform statistical testing.

As already stated in introduction, we have shown that this description of trajectories is robust to noise and short trajectory length, and accounts for several aspects of the dynamics, but it is entirely possible to follow the process described in the following paragraphs with another set of descriptive statistics if the features on which the analysis should be focused are well circumscribed. We verified that latent spaces obtained with models trained on different tasks and random walk models were exploitable as well, and hence that the results presented hereafter are not too dependent on this precise method of training. Results are shown in S5 Fig.

Maximum mean discrepancy.

Maximum mean discrepancy (MMD), introduced in [14], is a measure of distance between distributions. It was developed to perform statistical testing between two sets of independent observations lying in a metric space , X = {x₁, …, x_m} drawn from probability measure p and Y = {y₁, …, y_n} drawn from q, with the the goal of assessing whether or not p and q are different.

Given a class of functions from to , the MMD between two probability measures p and q is defined as (1) where E_x and E_y denote expectation w.r.t. p and q, respectively.

If the function class is the unit ball in a Reproducing Kernel Hilbert Space (RKHS) [14] , the square of the MMD can directly be estimated from data samples. Denoting k the kernel operator such that , an unbiased estimator of the square of the MMD between X and Y is given by: (2)

In our case, , and we used the classical Gaussian kernel . We set the kernel bandwidth σ either to the median of the pairwise Euclidian distances between samples from X and Y or we optimised it in specific conditions.

The MMD is capable of detecting subtle differences such as the ones between data generated by generative adversarial networks (GANs) and real data [53]. It has also proved efficient in discovering which variables exhibit the greatest difference between data sets [14, 54].

Statistical test.

We adapted the bootstrap test described in [14] to assess whether dynamics of two sets of trajectories exhibit significant differences. The sets of trajectories are represented by their two sets X and Y of descriptive vectors, respectively drawn from unknown probability densities p and q. For simplicity, we assume here that X and Y have the same number of elements, m = n, using the same notation as in the last section. In practice, the number of observed trajectories varies significantly across experimental replicates. Hence, to ensure that all replicates have an equal importance when the two sets do not have the same number of trajectories we randomly sub-sampled the larger of the two sets to equalise their sizes. The null hypothesis H₀ of the statistical test is that p = q, i.e. the two conditions lead to the same distribution of random walks. Under H₀, we approximated the distribution of by bootstrapping, i.e. we drew random samples from the union of X and Y and distributed them in two groups X′ and Y′ (whose sizes respectively match those of X and Y), on which we computed . We repeated this procedure a sufficient number of times to obtain an estimation of the distribution of under the assumption that X′ and Y′ are drawn from the same distribution. Then, if the original was greater than the 1 − β quantile of this distribution, we rejected H₀. This test is said to be of level β, because with probability β, we will reject the null hypothesis when it is actually true.

Some factors unrelated to the dynamics of a biomolecule might have an influence on its descriptive statistics: this is notably the case of the trajectory length, or of the localization uncertainty. To compensate these undesired differences, it is possible to “marginalize” for these factors prior to comparing two sets X and Y of descriptive vectors, by making sure that the distribution of such factors is equal in both sets. In the following, we always marginalized comparisons by trajectory length. In cases when a large number of conditions are analyzed, and if one is worried that p-values of comparisons might be close to the significance threshold (because of a low number of trajectories, or when the difference between conditions is expected to be subtle), we advise to carefully choose the comparisons to perform, in accordance with pre-existing biological knowledge, so as to limit the amplitude of the correction to apply to the significance threshold (e.g. using the Bonferroni method).