Langevin and Smoluchowski equations as a model of motion

Statistical methods to extract information from SPTs are based on stochastic models, such as the Langevin equation or its Smoluchowski's limit and associated models that account for additional localization point identification noise or memory kernel.^[14] The Langevin equation describes a stochastic particle driven by a Brownian force ${\displaystyle \Xi }$ and a field of force (e.g., electrostatic, mechanical, etc.) with an expression ${\displaystyle F(x,t)}$:

${\displaystyle m{\ddot {x}}+\Gamma {\dot {x}}-F(x,t)=\Xi ,}$

where m is the mass of the particle and ${\displaystyle \Gamma =6\pi a\rho }$is the friction coefficient of a diffusing particle, ${\displaystyle \rho }$ the viscosity. Here ${\displaystyle \Xi }$ is the ${\displaystyle \delta }$-correlated Gaussian white noise. The force can derived from a potential well U so that ${\displaystyle F(x,t)=-U'(x)}$ and in that case, the equation takes the form

${\displaystyle m{\frac {d^{2}x}{dt^{2}}}+\Gamma {\frac {dx}{dt}}+\nabla U(x)={\sqrt {2\varepsilon \gamma }}\,{\frac {d\eta }{dt}},}$

where ${\displaystyle \varepsilon =k_{\text{B}}T,}$ is the energy and ${\displaystyle k_{\text{B}}}$ the Boltzmann constant and T the temperature. Langevin's equation is used to describe trajectories where inertia or acceleration matters. For example, at very short timescales, when a molecule unbinds from a binding site or escapes from a potential well ^[15] and the inertia term allows the particles to move away from the attractor and thus prevents immediate rebinding that could plague numerical simulations.

In the large friction limit ${\displaystyle \gamma \to \infty }$ the trajectories ${\displaystyle x(t)}$ of the Langevin equation converges in probability to those of the Smoluchowski's equation

${\displaystyle \gamma {\dot {x}}+U^{\prime }(x)={\sqrt {2\varepsilon \gamma }}\,{\dot {w}},}$

where ${\displaystyle {\dot {w}}(t)}$ is ${\displaystyle \delta }$-correlated. This equation is obtained when the diffusion coefficient is constant in space. When this is not case, coarse grained equations (at a coarse spatial resolution) should be derived from molecular considerations. Interpretation of the physical forces are not resolved by Ito's vs Stratonovich integral representations or any others.

General model equations

For a timescale much longer than the elementary molecular collision, the position of a tracked particle is described by a more general overdamped limit of the Langevin stochastic model. Indeed, if the acquisition timescale of empirical recorded trajectories is much lower compared to the thermal fluctuations, rapid events are not resolved in the data. Thus at this coarser spatiotemporal scale, the motion description is replaced by an effective stochastic equation

${\displaystyle {\dot {X}}(t)={b}(X(t))+{\sqrt {2}}{B}_{e}(X(t)){\dot {w}}(t),\qquad \qquad (1)}$

where ${\displaystyle {b}(X)}$ is the drift field and ${\displaystyle {B}_{e}}$the diffusion matrix. The effective diffusion tensor can vary in space ${\displaystyle D(X)={\frac {1}{2}}B(X)B^{T}X^{T}}$ (${\textstyle X^{T}}$ denotes the transpose of ${\textstyle X}$). This equation is not derived but assumed. However the diffusion coefficient should be smooth enough as any discontinuity in D should be resolved by a spatial scaling to analyse the source of discontinuity (usually inert obstacles or transitions between two medias). The observed effective diffusion tensor is not necessarily isotropic and can be state-dependent, whereas the friction coefficient ${\displaystyle \gamma }$ remains constant as long as the medium stays the same and the microscopic diffusion coefficient (or tensor) could remain isotropic.

Empirical estimators for the drift and diffusion tensor of a stochastic process

Several empirical estimators have been proposed to recover the local diffusion coefficient, vector field and even organized patterns in the drift, such as potential wells.^[17] The construction of empirical estimators that serve to recover physical properties from parametric and non-parametric statistics. Retrieving statistical parameters of a diffusion process from one-dimensional time series statistics use the first moment estimator or Bayesian inference.

The models and the analysis assume that processes are stationary, so that the statistical properties of trajectories do not change over time. In practice, this assumption is satisfied when trajectories are acquired for less than a minute, where only few slow changes may occur on the surface of a neuron for example. Non stationary behavior are observed using a time-lapse analysis, with a delay of tens of minutes between successive acquisitions.

The coarse-grained model Eq. 1 is recovered from the conditional moments of the trajectory by computing the increments ${\displaystyle \Delta X=X(t+\Delta t)-X(t)}$:

${\displaystyle a(x)=\lim _{\Delta t\rightarrow 0}{\frac {E[\Delta X(t)\mid X(t)=x]}{\Delta t}},}$

${\displaystyle D(x)=\lim _{\Delta t\rightarrow 0}{\frac {E[\Delta X(t)^{T}\,\Delta X(t)\mid X(t)=x]}{2\,\Delta t}}.}$

Here the notation ${\displaystyle E[\cdot \,|\,X(t)=x]}$means averaging over all trajectories that are at point x at time t. The coefficients of the Smoluchowski equation can be statistically estimated at each point x from an infinitely large sample of its trajectories in the neighborhood of the point x at time t.

Empirical estimation

In practice, the expectations for a and D are estimated by finite sample averages and${\displaystyle \Delta t}$ is the time-resolution of the recorded trajectories. Formulas for a and D are approximated at the time step ${\displaystyle \Delta t}$, where for tens to hundreds of points falling in any bin. This is usually enough for the estimation.

To estimate the local drift and diffusion coefficients, trajectories are first grouped within a small neighbourhood. The field of observation is partitioned into square bins ${\displaystyle S(x_{k},r)}$of side r and centre ${\displaystyle x_{k}}$ and the local drift and diffusion are estimated for each of the square. Considering a sample with ${\displaystyle N_{t}}$ trajectories ${\displaystyle \{x^{i}(t_{1}),\dots ,x^{i}(t_{N_{s}})\},}$ where ${\displaystyle t_{j}}$ are the sampling times, the discretization of equation for the drift ${\displaystyle a(x_{k})=(a_{x}(x_{k}),a_{y}(x_{k}))}$at position ${\displaystyle x_{k}}$ is given for each spatial projection on the x and y axis by

${\displaystyle a_{x}(x_{k})\approx {\frac {1}{N_{k}}}\sum _{j=1}^{N_{t}}\sum _{i=0,{\tilde {x}}_{i}^{j}\in S(x_{k},r)}^{N_{s}-1}\left({\frac {x_{i+1}^{j}-x_{i}^{j}}{\Delta t}}\right)}$

${\displaystyle a_{y}(x_{k})\approx {\frac {1}{N_{k}}}\sum _{j=1}^{N_{t}}\sum _{i=0,{\tilde {x}}_{i}^{j}\in S(x_{k},r)}^{N_{s}-1}\left({\frac {y_{i+1}^{j}-y_{i}^{j}}{\Delta t}}\right),}$

where ${\displaystyle N_{k}}$is the number of points of trajectory that fall in the square ${\displaystyle S(x_{k},r)}$. Similarly, the components of the effective diffusion tensor ${\displaystyle D(x_{k})}$ are approximated by the empirical sums

${\displaystyle D_{xx}(x_{k})\approx {\frac {1}{N_{k}}}\sum _{j=1}^{N_{t}}\sum _{i=0,x_{i}\in S(x_{k},r)}^{N_{s}-1}{\frac {(x_{i+1}^{j}-x_{i}^{j})^{2}}{2\,\Delta t}},}$

${\displaystyle D_{yy}(x_{k})\approx {\frac {1}{N_{k}}}\sum _{j=1}^{N_{t}}\sum _{i=0,x_{i}\in S(x_{k},r)}^{N_{s}-1}{\frac {(y_{i+1}^{j}-y_{i}^{j})^{2}}{2\,\Delta t}},}$

${\displaystyle D_{xy}(x_{k})\approx {\frac {1}{N_{k}}}\sum _{j=1}^{N_{t}}\sum _{i=0,x_{i}\in S(x_{k},r)}^{N_{s}-1}{\frac {(x_{i+1}^{j}-x_{i}^{j})(y_{i+1}^{j}-y_{i}^{j})}{2\,\Delta t}}.}$

The moment estimation requires a large number of trajectories passing through each point, which agrees precisely with the massive data generated by the a certain types of super-resolution data such as those acquired by sptPALM technique on biological samples. The exact inversion of Lagenvin's equation demands in theory an infinite number of trajectories passing through any point x of interest. In practice, the recovery of the drift and diffusion tensor is obtained after a region is subdivided by a square grid of radius r or by moving sliding windows (of the order of 50 to 100 nm).

Automated recovery of the boundary of a nanodomain

Algorithms based on mapping the density of points extracted from trajectories allow to reveal local binding and trafficking interactions and organization of dynamic subcellular sites. The algorithms can be applied to study regions of high density, revealved by SPTs. Examples are organelles such as endoplasmic reticulum or cell membranes. The method is based on spatiotemporal segmentation to detect local architecture and boundaries of high-density regions for domains measuring hundreds of nanometers.^[18]