Biophysics of high density nanometer regions extracted from super-resolution single particle trajectories: application to voltage-gated calcium channels and phospholipids

The cellular membrane is very heterogenous and enriched with high-density regions forming microdomains, as revealed by single particle tracking experiments. However the organization of these regions remain unexplained. We determine here the biophysical properties of these regions, when described as a basin of attraction. We develop two methods to recover the dynamics and local potential wells (field of force and boundary). The first method is based on the local density of points distribution of trajectories, which differs inside and outside the wells. The second method focuses on recovering the drift field that is convergent inside wells and uses the transient field to determine the boundary. Finally, we apply these two methods to the distribution of trajectories recorded from voltage gated calcium channels and phospholipid anchored GFP in the cell membrane of hippocampal neurons and obtain the size and energy of high-density regions with a nanometer precision.


MLE estimator for a potential well
We modified the Maximum likelihood Estimator (MLE) procedure to reconstruct from SPTs, the geometrical parameters (center and covariance matrix) of a well. Using a Ornstein-Ulhenbeck process, we apply the MLE procedure to the points of the trajectories falling inside the ensemble where ρ is the steady-state probability density function of the OU process. The advantage of the Maximum-likelihood approach is that no spatial discretization is needed. We recall that the transition probability density of an OU process centered at µ, with diffusion coefficient D and spring constant λ The transition probability is where a trajectory is discretized in y 1 , . . . , y M , with a fixed time step ∆t. The log likelihood is The maximum-likelihood approach consists in estimating λ that maximizes the log-likelihood l(y 1 , . . . , y M ). We change variables x = e −λ∆t and v = λ 2πD(1−e −2λ∆t ) so that l(y 1 , . . . , y M |x, At the maximum, leads to the coupled equationŝ This procedure can be generalized in two dimensions and in addition, we apply the estimator of eq. (9) to the ensemble Γ α , defined by (1). We thus obtain the following estimatorŝ where M α is the number of points y i ∈ Γ α . We apply this estimator to numerical simulations in Figs. S1 and S2 and to CaV2.2 and GPI-GFP SPTs in Figs. 7 and 8 (main text) respectively. To recover the potential well from the drift distribution, we use a least square estimator where N is the number of points X i = (x i , y i ) and the potential well is so that The minimizers are given by from which, we obtain the center and the eigenvalues of the covariance matrix: In practice, we computed the center µ x , µ y and the eigenvalues λ x , λ y over the points X i falling inside the well.

Center location and semi-axes from optimal fit
We derive here a close formula for the eigenvalues and the center associated to the optimal estimators of the drift. For an OU process, we recall that the eigenvalues are given by and eq. (20) in eq. (21), we obtained where µ y is solution of the quadratic equation Similarly, we get Relations (24) and (25) lead to a close expression of the eigenvalues (eq. (19)).

Comparing MLE with density estimators
To recover the center and eigenvalues of an Ornstein-Uhlenbeck process with the potential well given in eq. (14), we apply the MLE procedure (section 1) for various values of the parameter α. We compare the results with the density (section 3 main text) and the LSQ (section 2) methods. We find that the MLE and density approaches are quite robust and give similar results, as shown in Fig. S1: Interestingly, all three estimators allow to recover the center (µ x , µ y ) with high accuracy for the disk and the ellipse, when 100% to 50% of the points are taken into account α ∈ [0; 0.5] (Fig. S1). However, the estimation of the eigenvalues is acceptable for the MLE only, in the range α ∈ [0; 0.5], because it diverges in the two other cases, except when α 0.1 (90% of the distribution is used). When the time step δt of the numerical simulations of eq. (2) and the sampling time ∆t are equal, the three methods lead to a good recovery of the centerμ x ,μ y (Fig. S2), but differ for recovering the eigenvalues (λ x ,λ y ). The least square approach is less dependent on the parameter α than the two others. Probably because the distribution was generated with a large time step so that the statistics are calculated on trajectories far from the equilibrium. This result shows that the least square approach does not require to sample over a steady state distribution and thus recovering the parameters from the drift is possible for a large range of the parameter α. In section 4, we will estimate the effect of changing the time steps.

Influence of the time and spatial discretizations on the Least Square Estimation
In this section, we shall estimate the impact of the initial points distribution of the trajectories on the estimation of the drift. For the stochastic equation the optimal estimator for the drift A at a time resolution ∆t is obtained by the formula [8,5], where X n = X(n∆t) and E [.] in eq. (27) is the expectation. When a grid of size ∆x is used to estimate the drift map, all points in bin k leads to the same drift a ∆t (x k ), where x k is the center of the bin. For long trajectories, the stochastic process samples the steady-state distribution p(x). Such distribution might influence the computation of the drift inside bin k. To estimate this contribution, we normalize the steady-state distribution q(x) of the stationary process by The estimated drift a ∆t (x) depends on the distribution of points falling into the bin ∆(x) as follows where p ∆t (y|x) is the pdf to find X(t) at the point y at time t + ∆t when it started at point x at time t. In the small time limit, We shall now obtain a further approximation by using a Taylor's expansion of the function F (x) = ∫ x 0 p(s)a(s)ds. We obtain that ∫ which leads to the approximation: ) .
The formula in higher dimension is given by ) .
We conclude that a discretization in bins of size ∆x perturbs the drift recovery by a term (∆x) 2 .

Influence of the time discretization ∆t on the drift estimation
To study the consequences of a discrete sampling time on the reconstruction of the drift from SPTs, we focus on the one-dimensional OU-process where λ, µ are fixed. A direct integration of equation (34) for s ≤ t leads to and for two consecutive points x(t) and x(t + ∆t), we have Thus, when ∆t is small, the drift at position x and resolution ∆t is We conclude that at resolution ∆t, the approximation error is suggesting that the drift of an OU process is always under-estimated using the displacement estimator.

Time and space discretization for Ornstein-Ulhenbeck process
We shall now evaluate the cumulative effect of a temporal ∆t and spatial ∆x discretization on the recovery of an OU-process. The spatial grid of size ∆x and the drift at position x are estimated empirically using the points falling in the bin ∆(x) = [x − ∆x/2, x + ∆x/2]. We start with the conditional steady-state distribution q ∆ (x) of points falling in ∆(x), which is linked to the pdf p(x) of the OU-stationary process by The drift term from eq. (37) can be approximated as where µ is the center of the OU-process and the stationary pdf is given by To estimate eq. (40), we use eq. (31). For this computation, we set µ = 0.
In that case, we have 11 In the small ∆x approximation, we have Using eqs. (40) and (43), we obtain an approximation for the drift at finite time step ∆t and grid size ∆x To conclude, relation (44) reveals that the empirical displacements x(t + ∆t) − x(t) collected over trajectories for an OU, can be used to recover the drift, with an additional exponential order correction in ∆t and a second order in ∆x.

Empirical estimations of the drift
The empirical estimatorã of the drift at position x for finite time ∆t and spatial steps ∆x, is defined where N is the number of points x i (t j ) located in the square bin of infinitesimal surface (∆x) 2 around x. Using eq. (44), the approximation at second order in ∆x gives that To recover the parameter λ at various order of ∆x, we can expressλ using a regular series expansioñ The first term is obtained by setting ∆x = 0 in eq. (46), leading for any bin centered around In that case, a linear regression method can be used using the two coordinates x k − µ andã and to fit the distribution with a line and invert eq. (48). If necessary, the next in the expansion can be found. Note that a formal inversion of eq. (48) shows that for each k, we can obtain an estimation for λ 0 : so thatλ showing that numerical fluctuations inã(x k ) for |x k − µ| small can drastically affect the estimation. We use this result to study the recovery of the parameters in Fig. 5C (main text) and Fig. S3C, where we indeed observe larger errors near the center of the well than inside. We refer to Fig. 6 (Main text) for the estimation of the eigenvalue with and without the center bin.

Effect of the grids intersecting the boundary in the estimation of the drift
The recovery of a truncated OU involves estimating several parameters that depend on the accurate detection of the boundary. We focus here on the drift estimation for the part of the square grid that intersects the boundary (green bins in Fig. 6, main text). In that case, for the interior part that intersects the elliptic domain, the empirical estimation recovers the local vector, while outside, it fluctuates around zero, due to the nature of the Brownian motion (no drift). Thus the error of the drift estimation at the boundary increases with the area fraction of the bin falling outside the domain. To estimate this error, we recall that the truncated OU-process is defined bẏ where Since the drift is zero for a diffusion process, located outside Γ E , we have where ∆ 1 (x) = ∆(x) ∩ Γ E is the part of the grid interior to the ellipse. Indeed, the drift a diffusion process is zero. In addition, we suppose that the sample is made according to a normalized distribution where p(x) is any distribution that could be the steady-state distribution of a truncated OU inside the well and is uniform outside. From eq. (53), we get In first approximation, where ∆ 2 (x) = ∆(x) ∩ Γ E ∆(x). We shall now estimate ∆ 1 (x). We first note that the conservation of surfaces: ∆ 1 (x) + ∆ 2 (x) = (∆x) 2 . For a square centered at a boundary of the ellipse (x, y) we consider the square grid with integer coordinates k = [ x ∆x ] and q = [ y ∆x ]. To compute ∆ 2 (x), we subtract the total area of the rectangle to the surface underneath the ellipse between the point k∆x and ∆x + k∆x: .
Thus the computations lead to and for x > ∆x, we get A similar computation leads for 0 < x ≤ ∆x to To conclude, except for the bin at the four extreme positions of the ellipse, the error is of order O((∆x) 2 ) for each grid bin, leading to a cumulative error along the total length of O((∆x)). The error contribution is shown in Fig.  6E.

Influence of the time steps in stochastic simulations
In the Smoluchowski's limit of the Langevin equation, the first order stochastic equation from which trajectories are generated is obtained by choosing a time step δt. This time step should not be smaller than the reciprocal of the friction coefficient γ so that the successive points X(δt), X(2δt), ..X(nδt), ... should approximate the physical trajectory. When the sampling rate is such that ∆t ≫ δt, we can compare the drift estimation in that case and also study the extreme case when the sampling and simulation time steps are identical, leading to a jump process. When ∆t ≫ δt, the drift is computed after n steps. Using the empirical estimator, in the limit of a large number of trajectories N , we get where p s (y|x) is the pdf of the process X(t) starting at x at time 0 and ending at y at time s. This result is quite different from the classical estimator of eq. (29). We note that the result of eq. (64) is very different from the estimation from a single observation time step ∆t. Using the pdf an Onrstein-Uhlenbeck process we get Using a Taylor's expansion in the drift term: Thus the estimator of the drift with many time steps is at first order in ∆t given by We now evaluate the consequences of simulating a process with many time steps δt = 10 −4 s whereas the sampling time was ∆t = 0.02s in Fig. 5 (main text). We show in Fig. S3 how the drift field can be recovered. This situation corresponds to large jumps of the underlying physical process. To conclude, we find that the center and peripheral grid bins are generating most of the error, especially for large grid sizes (∆x = 50 and 90 nm).

Conditional drift estimation
In this last section, we discuss the effect of estimating the drift by conditioning the end points of a displacement to stay inside the potential well. Computing the displacement ∆X by selecting only trajectories that stay inside the well gives a bias estimator of the drift. This situation appears when the trajectories never reach the boundary. The drift estimator is computed from the conditional pdf p * of the process that stays inside the potential well. To find such a drift, we introduce the probability that a stochastic particle hits a ball of radius ϵ centered on the well before escaping from the well [9], then where p is the pdf in the entire space. q is solution of  Figure S3: Recovering the vector field with an equal simulated and sampled time step δt = ∆t A Recovering the local drift field inside a circular well for different grid sizes (10 nm, 50 nm, 90 nm) using numerical simulations with a sampling ∆t = 20 ms, with the constraints that at least 10 points falls inside a bin. B Error between the true and observed fields averaged over all the square bins inside the well versus the time step ∆t. C Error between the true and observed fields averaged over the radial angle versus the distance r to the center for various timestep (see color code). L * is the backward Fokker-Planck equation associated the process X(t) (p.77 [9]), defined by In that case, To conclude, by restricting the computation of the displacements to empirical trajectories that only remain in the well, an additional term has to be accounted for, which diverges as the distance from the point x to the boundary tends to zero (Fig. S4). We estimated the drift for displacements that do not exit the well (assuming the boundary is known).