Fokker-Planck analysis of superresolution microscopy images

A method for the analysis of superresolution microscopy images is presented. This method is based on the analysis of stochastic trajectories of particles moving on the membrane of a cell with the assumption that this motion is determined by the properties of this membrane. Thus, the purpose of this method is to recover the structural properties of the membrane by solving an inverse problem governed by the Fokker-Planck equation related to the stochastic trajectories. Results of numerical experiments demonstrate the ability of the proposed method to reconstruct the potential of a cell membrane by using synthetic data similar those captured by superresolution microscopy of luminescent activated proteins.


Introduction
The pioneering works [5,11,21] mark the development of the revolutionary superresolution microscopy (SRM) that allows to go beyond the Abbe limit for conventional light microscopy [1].The SRM method consists in labelling the molecules moving on a biological support with fluorophores and then in sampling the microscopic images of the activated fluorescent molecules.
Observation of the frames of the sampled SRM microscopic images have suggested that the motion of the molecules could be modelled by a stochastics Langevin equation [24,26].Specifically, it appears that an adequate model of the observed trajectories of 2-dimensional images is given by the following stochastic differential equation (SDE) [12]: where b represents the drift, σ is the dispersion coefficient, and X t ∈ R 2 denotes the position of the observed molecule at time t.In this framework, it is well-known that the drift and dispersion coefficients satisfy where the expected values are computed with respect to the process having value z at t = s; the operator E[• | X s = z] denotes averaging w.r.t the measure of the trajectories conditioned to be at z at time s.
The formulas above suggest that suitable approximations to b and σ can be obtained by tracking single molecules; see, e.g., [9,12].However, this approach may suffer of the highly fluctuating values of the trajectories and the difficulty of discerning between different molecules that come closer than the resolution limit.
For this reason, already in [2] the authors have pursued an alternative strategy that allows to build a robust methodology for estimation of the drift based on the observation of ensemble of trajectories.Our approach is built upon the assumption that this ensemble is driven by a velocity field (the drift), given by a potential velocity field U (x), as follows: b(x; U ) = −∇U (x). ( Moreover, one assumes a constant diffusion coefficient whose value is chosen consistently with estimates of laboratory measurement [20].Further, in agreement with our statistical approach based on ensembles, we focus on the evolution of the probability density function (PDF) of the positions of the molecules (not on the single trajectories) whose evolution is governed by the Fokker-Planck (FP) problem given by [23,26]: where Q = Ω × (0, T ) and Σ = ∂Ω × (0, T ).In this formulation, f (x, t) represents the PDF of a particle at x ∈ R 2 at time t, ∇U (x) is the Cartesian gradient of the potential U , f 0 is the initial density, and ∆ is the twodimensional Laplace operator.Notice that we require zero-flux boundary conditions, where F (f )(x, t) is the following flux of probability We choose zero-flux boundary conditions since they reasonably model the situation where a similar number of particles enters and exits the domain; see, e.g., [4,26].
Our proposal is to construct a FP-based imaging modality that is based on the formulation of an inverse problem for U and the observation of a time sequence, in a time interval [0, T ], of numerical PDFs (two-dimensional histograms), which are obtained from a uniform binning of SRM particles' positions.We denote this input data with f d (x, t) which is a piecewise constant function.In this setting the initial condition is given by f This proposal is similar to that in our previous work [2].However, in [2] the assumption of interacting particles was made that resulted in a nonlinear FP model and very involved and CPU time demanding calculations.
At the continuous level, our FP-based imaging tool is formulated as the following inverse problem: with the given initial-and boundary conditions for the FP equation, and α, ξ > 0.
In this problem, the objective functional J is defined as the weighted sum of a space-time bestfit term Ω T 0 (f (x, t) − f d (x, t)) 2 dx dt, and at final time Ω (f (x, T ) − f d (x, T )) 2 dx, and of a suitable 'energy' of the potential ∥U ∥ 2 = Ω (|U (x)| 2 + |∇U (x)| 2 ) dx, which corresponds to the square of the H 1 (Ω) norm of U .Notice that this formulation allows to avoid any differentiation of the data and makes possible to choose the binning size and, in general, the measurement setting, independently of any choice of parameters that are required in the numerical solution of the optimization problem.
Our second main concern in determining the potential U is to provide a measure of uncertainty, and thus of reliability, of its reconstruction.Statistically, this is achieved by many repetition of the same experiment, that could not be feasible for (short) living cells.However, inspired by the so-called model predictive control (MPC) scheme [10] already used for optimal control problems [3,4], we propose a novel procedure to quantify the uncertainty of the estimation of U by using the data of a single experiment.
Our methodology is to consider a sequence of non-parametric inverse problems like (8) defined on time windows (t k , t k+1 ), k = 0, . . ., K − 1, that represent an uniform partition in K subintervals of the time interval [0, T ].Therefore a statistical analysis can be performed on the set of the corresponding K solutions for U that are obtained in the subintervals.
For development and validation, we consider images of cell's membrane structures (actin, cytoskeleton), as expression of potentials, that is pixel grey values, with which we generate our synthetic data.In particular, we use an image of actin from a cytoskeleton obtained with a Platinium-replica electron microscopy [8,17].
With this images taken as gray-level representation of potential functions, we perform Monte Carlo simulation of motion of particles to generate images of molecules at different time instants, thus constructing the datasets representing the output of measurements.This setting is illustrated in Figure 1, where the image of actin [14] and a plot of a few trajectories of the corresponding stochastic motion of the particles in this potential are shown.
Once the synthetic measurement data is constructed, we perform a preprocessing step on this data to construct the numerical PDF required in our method and solve our inverse problem to find the estimated-measured potential U .The latter is compared with that one used in the MC simulation, by a measure of similarity based on the pixel cross-correlation between the two images.
In Section 2, we discuss a numerical methodology for solving our FPbased reconstruction method for the potential U .In Section 3, we provide all details of our experimental setting and introduce some analysis tools for determining the accuracy of the proposed reconstruction.In Section 4, we validate our reconstruction method, and use our uncertainty quantification procedure.In Section 5, we investigate the resolution of the proposed FP image reconstruction as optical instrument.A section of conclusion and acknowledgements completes this work.

Numerical methodology
Our aim is to reconstruct the potential U from the data consisting of a temporally sampled SRM images of the positions of particles subject to this potential; see Figure 2 (left), for a schematic snap-shot of this data.This image is subject to a pre-processing binning procedure in order to construct histograms by counting the number of particles in a regular square partition of Ω.The height of an histogram is proportional to the number of particles in a bin of the domain.This procedure for the image at time t defines the histogram function f d (•, t); see Figure 2 (right).
In order to illustrate our numerical framework, we introduce the potential- to-state map U → f = S(U ), that is, the map that associates to a given U ∈ H 1 (Ω) the unique solution to our FP problem ( 4)-( 6), with given initial condition f 0 .For the analysis of well-posedness and regularity of the map S we refer to [13].
Next, we remark that with the map S, we can define the reduced objective functional Ĵ(U ) := J(S(U ), U ) and consider the equivalent formulation of (8) given by min which has the structure of an unconstrained optimization problem.Thanks to the regularity of S and the quadratic structure of J, existence of an optimal U can be stated by well known techniques; see, e.g., [18].Further, since S and J are Fréchet differentiable, it is possible to characterize an optimal U as the solution to the following first-order optimality condition where ∇ U Ĵ(U ) denotes the so-called reduced gradient [6].
In the Lagrange framework, this condition results in the following opti-mality system: where p denotes the adjoint variable, which is governed by a backward adjoint FP equation.The last equation in (10) is the so-called the optimality condition equation, and the Neumann boundary condition ∂ nU = 0 is our modelling choice.One can show that its left-hand side represents the L 2 gradient along the FP differential constraint with respect to U of the objective functional.We have Our approach for solving our FP optimization problem ( 8) is based on the nonlinear conjugate gradient (NCG) method; see, e.g., [6].This is an iterative method that resembles the standard CG scheme and requires to estimate the reduced gradient ∇ U Ĵ(U ) at each iteration.
In order to illustrate the NCG method, we start with a discussion on the construction of the gradient.For a given U n obtained after n iterations, we solve the FP equation and its adjoint, and use (11) to assemble the L 2 gradient.However, since the potential is sought in H 1 (Ω), we need to obtain the H 1 gradient that satisfies the following relation where (•, •) denotes the L 2 (Ω) scalar product.Now, using the definition of the H 1 inner product, we obtain which must hold for all the test functions δU ∈ H 1 (Ω).Therefore we obtain with the boundary conditions ∂ ∂ n ∇ U Ĵ(U )| H 1 = 0 on ∂Ω; see [2] for more details.In Algorithm 1 our procedure for computing the gradient is given.
Solve forward the FP equation with inputs: f 0 (x), U (x) Solve backward the adoint FP equation with inputs: In this algorithm, the FP problem and its optimization FP adjoint are approximated by the exponential Chang-Cooper scheme and the implicit BDF2 method; see [4].Now, we can discuss the NCG method.The NCG iterative procedure is initialized with U 0 (x) = 0. We denote the optimization directions with d n .In the first update, we have where α 0 is obtained by a backtracking linesearch procedure.After the first step, in the NCG method the descent direction is defined as a linear combination of the new gradient and the past direction as follows: where β −1 = 0, and that is, the Dai-Yuan formula; see, e.g., [6].
The tolerance tol and the maximum number of iterations n max are used for termination criteria.Summarizing, in Algorithm 2 we present the NCG procedure.For our numerical experiments these algorithms have been implemented with object oriented programming in C++, by using the numerical libraries Armadillo [25], Lapack [15] and SuperLU [7].

Experimental design and analysis tools
In a real application, the input data f d is given by frames of a recorded sequence of a SRM experiment, and the desired output is the reconstructed potential denoted with U r .In our case, we construct this data based on a Once we have computed U r with our optimization procedure, we aim at providing a quantification of its uncertainty.Thus, we compute the following normalized cross-correlation factor between the reconstructed potential U r and the one used to generate the synthetic data U s .We have In this formula, U r and U s are considered as vectors and • represents the scalar product.Therefore if cc = 1 we have that the two potentials match perfectly, whereas if its value is close to 0 the two potentials are dissimilar, and 0.5 it is poor.Notice that cross-correlation is commonly used in medical imaging and biology; see, e.g., [16,22,27].
Clearly, one could consider many repetition of the simulation of the motion of the particles with the same initial condition and make the final binning on the average of the resulting frames.This procedure would result in a less fluctuating f d (x, t ℓ ) that allows a better reconstruction.However, this scenario seems difficult to realize in the real laboratory setting of a living cell.On the other hand, in SRM imaging is possible to visualize the motion of the particles on a cell membrane for a relatively long time (T ≫ 1 in our setting) and our approach exploits this possibility considering a subdivision of the time interval in a number K of time windows, and solving our optimization problem in each of these windows almost independently.This approach allows to improve the reconstruction U r and makes possible to quantify the uncertainty of the reconstruction.Now, to illustrate our approach, consider a uniform partition of [0, T ] in time windows of size ∆t = T /K with K a positive integer.Let t k = k∆t, k = 0, 1, . . ., K, denote the starting-and end-points of the windows.At time t 0 , we have the initial PDF f 0 , and we solve our optimization problem (8) in the interval [t 0 , t 1 ].This means that the final time is t 1 and the terminal condition for the adjoint variable is given by p(x, t 1 ) = −ξ (f (x, t 1 ) − f d (x, t 1 )).The resulting potential is denoted with U 1 .Thus, the solution obtained in this window provides also the PDF at t = t 1 .
Clearly, we can repeat this procedure in the interval (t 1 , t 2 ) with the computed PDF at t = t 1 as the initial condition and t 2 as final time, to compute U 2 .This procedure is recursive and can be repeated for k = 1, . . ., K, thus obtaining U k , k = 1, . . ., K.
Notice that small values of K in relation to L produce a rough estimate of the average potential and its standard deviation due to statistical fluctuations of the Monte Carlo experiments.On the other hand, for greater values of K, the number of frames for each window of our approach is reduced when L is kept fixed, thus resulting in a worsening of the reconstruction procedure.
For the purpose of our analysis, we apply a scaling of these potentials so that their point-wise values are in the interval [0, 1].This scaling is performed as follows: Thereafter, we the reconstructed potential by pixel-wise average of the U k is given by Moreover, we can also compute the following pixel-wise standard deviation Next, we provide conversion formulas for our parameters in order to accommodate data from real laboratory experiments.We introduce a unit of length u such that the side length l of our square domain Ω is given by l = 6 u, and the unit of the noise amplitude σ is given by √ u/s.In real biological experiments, the typical measure of the length l of a cell membrane is given in µm.Further, the particle's diffusion constant D = σ 2 /2 is given in µm 2 /s, hence we have the correspondence σ = l/ l√ 2D in unit √ u/s, whose value is used for MC simulations.
The depth of the potential Ũ is expressed in unit of K B T , where K B is the Boltzmann constant and T the absolute temperature.In experimental papers, the equation ( 3) is written with the diffusion constant D and K B T , i.e.D Ũ /(K B T ).As above, we obtain the relationship between the values of a potential U and the scaled Ũ , as Ũ = U ( l/l) 2 /D in the unit of K B T .
As an illustration of the setting above, we see that in an experiment, the superresolution of an acquired image frame can reach the value of 0.02 µm/pixel.With an image of 500 × 500 pixels, we have l = 10 µm.The average diffusion coefficient of particles (protein molecules) observed in SRM imaging is estimated with D = 0.1 µm 2 /s [20].By superresolution techniques, it is possible to activate a density of 0.5 ÷ 2/µm 2 visible particles, which in terms of image pixels corresponds to 0.5 ÷ 2 particles in a square of 50 pixels of side.Each frame is usually sampled at time intervals of δt = 30 ms.
In order to set up a consistent MC simulation of a real experiment, by mapping an image of a square of side 10 µm on our domain Ω, we get from the above mentioned formula: σ = 6/10 √ 0.2 ≈ 0.268 u/ √ s.

Numerical validation
In this section, we discuss results of experiments in a setting that is close to real laboratory experiments involving SRM imaging.The results of these experiments demonstrate the ability of our methodology to reconstruct the potential from the simulations of the SRM measurements of the motion of particles on a cell's membrane.
We consider a potential that corresponds to a portion of cytoskeleton as depicted in Figure 3, with 200 × 200 pixels.We assume that the pixel is 50 nm, which corresponds to an area of 100 µm 2 .In the figure, the white regions represent the structure of the cytoskeleton; the black ones are the 'valleys' where the proteins are supposed to be attracted.
For the MC simulations for generating the synthetic data, we choose σ = 0.268 u/ √ s.This value of σ corresponds to a diffusion constant of D ≃ 0.1 µm 2 /s.We consider N p = 1000 particles, i.e. an average density of 10 particles per µm 2 .In this case, we consider a sequence of L = 3000 frames and T = 90, obtained by the numerical integration of the stochastic differential equation with an integration step τ = 30 • 10 −3 s.The frames have δt = 30 ms, similar to a real experiment.The resulting (single run) particle trajectories are collected in a binning process based on a mesh Ω of 50 × 50 bins.
For our reconstruction method, we choose a numerical partition of Ω of 100 × 100 subdivisions, corresponding to a mesh size of 100 nm.The time integration step coincides with that of the frames.For the tracking functional, we set α = 10 −4 and ξ = 1.Further, in the FP setting, we have σ = 0.7 u/ √ s.Notice that σ in the FP model is chosen larger than the one used in the MC simulations.This choice is dictated by numerical convenience and it appears that it does not affect the quality of the reconstruction.The calculations are performed according to the MPC procedure with K = 5 time windows.
With this setting, we obtain the reconstructed potential shown in Figure 3 (right).We see that the reconstruction is less sharp as we expected considering the much finer structure of the cytoskeleton and the small number of particles involved.
Further, in Figure 4, we depict the potentials obtained on each time window of the MPC procedure and the values of the corresponding cc.With these results, we have obtained the reconstructed potential U r in Figure 3, which we re-plot in level-set format in Figure 5 for comparison.In Figure 5, we also depict the standard deviation that suggests that we have obtained a reliable reconstruction with small uncertainty.

Resolution of FP-based image reconstruction
In this section, we investigate the optical resolution of our reconstruction method, that is, try to determine a confidence value related to the scale at which our method can resolve variations of the potential.As a guideline, we remark that single molecule localization microscopy (SMLM) can distinguish distances of molecules of approximately 20 nm resolution.Therefore we assume this resolution range of the fluorescently labelled particles images, and attempt to quantify the smallest scale at which geometric features of the reconstructed potential U can be distinguished.
For our purpose, we consider the following 'target potential', appearing as an alternating sequence of black and white circles (likewise those in test targets used for the resolution measurement of optical instruments), to synthetically generate the motion data of particles.We have where A denotes the semi-amplitude of the variation between the minimum and the maximum of the potential, l is the length of the side of the domain, d is the distance between two peaks of the potential as a fraction of l.Now, we consider a single MC simulation of 500 particles with the setting: σ = 0.5, T = 90 and L = 3000 frames, integrated with the time step τ = 0.03.In Figure 6, we show (left) the given potential with A = 0.05 and d = 1/20, with a gray-scale value representation conveniently adjusted for illustration pourpose.According to the above working hypothesis, we suppose that the pixel's width of the image is 20 nm In Figure 6 (left), we depict U in a square of side of 500 pixels, corresponding to l = 10 µm.Hence, the distance between two peaks is λ = 10/20 = 500 nm.Further, the particle's density is 5 particles per µm 2 , the diffusion coefficient D ≃ 0.3472 µm 2 /s, and the potential depth, i.e. the difference between the maximum and the minimum, is Ũ = 0.8 K B T .For the reconstruction process, we use a binning of 50 × 50, α = 10 −4 , ξ = 1.In the numerical setting, we use a grid of 100 × 100 points, and K = 5.Also in Figure 6 (right), we show the reconstructed potential ⟨U r ⟩ and notice its high accuracy that is also confirmed by the high value 0.82 of the cross correlation.Notice, that the quality of the reconstruction can be further improved by using post-processing techniques of images.Now, with the aim to define a criteria to establish the resolution measure for the potential, we introduce a confidence level for the quality of the reconstructed potential by setting a threshold for the calculated cross-correlation.This approach has been adopted in [16] for the detection of cellular objects from images acquired from electron tomography.For that purpose the authors used the threshold value of 0.5, whereas in our case, we set a more strict threshold-cc level equal to 0.8.With this threshold, we can state that the test pattern depicted in Figure 6 (left) is satisfactorily reconstructed and determine that the resolution measure associated to our 'imaging instrument' is 500 nm.Notice that this value is affected by the value of the potential U and the diffusion D, and it can be further improved by changing the other parameters of the experiment, such as K or the time T of the motion sampling.

Conclusion
A novel method for the analysis of superresolution microscopy images was presented, and applied to the reconstruction of the structure of a cell membrane potential based on observation of the motion of particles on the membrane.
The working principle of this method is the modeling with the Fokker-Planck equation of the ensemble of the stochastic trajectories of particles moving on the membrane of a cell, and the solution of an optimization problem governed by this equation, where the purpose of the optimization is to find a potential such that a least-squares bestfit term of the computed and ob-served particles' density and a Tikhonov regularization term are minimized.
Results of numerical experiments were presented that demonstrated the ability of the proposed method to reconstruct the potential of a cell membrane by using data of superresolution microscopy of luminescent activated proteins.

Figure 1 :
Figure 1: A picture of actin from a cytoskeleton as cell membrane potential (close up) (Courtesy of Koch Institute [14]) (left); a few simulated trajectories of particles (black dots) on the membrane (in reverse colors).

Figure 2 :
Figure 2: A frame of particles (left) and the corresponding histogram f d (x, t) on a mesh of 40x40 bins for a fixed time, from simulated data.

Figure 6 :
Figure 6: Left: the potential (19) with A = 0.05 and d = 1/20 (the gray scale levels spans from U = 0 to 0.1).Right: result of the reconstruction with the gray levels expanded to the min/max of ⟨U r ⟩.The cross correlation between the two images is 0.82.