Interactive simulation and visualization of point spread functions in single molecule imaging

The point spread function (PSF) is fundamental to any type of microscopy, most importantly so for single-molecule localization techniques, where the exact PSF shape is crucial for precise molecule localization at the nanoscale. However, optical aberrations and fixed fluorophore dipoles can lead to non-isotropic and distorted PSFs, thereby complicating and biasing conventional fitting approaches. In addition, some researchers deliberately modify the PSF by introducing specific phase shifts in order to provide improved sensitivity, e.g., for localizing molecules in 3D, or for determining the dipole orientation. For devising an experimental approach, but also for interpreting obtained data it would be helpful to have a simple visualization tool which calculates the expected PSF for the experiment in mind. To address this need, we have developed a comprehensive and accessible computer application that allows for the simulation of realistic PSFs based on the full vectorial PSF model. It incorporates a wide range of microscope and fluorophore parameters, enabling an accurate representation of various imaging conditions. Further, our app directly provides the Cramer-Rao bound for assessing the best achievable localization precision under given conditions. In addition to facilitating the simulation of PSFs of isotropic emitters, our application provides simulations of fixed dipole orientations as encountered, e.g., in cryogenic single-molecule localization microscopy applications. Moreover, it supports the incorporation of optical aberrations and phase manipulations for PSF engineering, as well as the simulation of crowded environments with overlapping molecules. Importantly, our software allows for the fitting of custom aberrations directly from experimental data, effectively bridging the gap between simulated and experimental scenarios, and enhancing experimental design and result validation.


Introduction
Single-molecule localization microscopy (SMLM) techniques offer a powerful approach to discern molecular structure and dynamics of biological samples below the diffraction limit [1].The overall idea is to virtually dilute single molecule signals in time, e.g., by stochastic switching of fluorophores between a bright and a dark state [2,3], by transient binding of fluorescent ligands [4,5], or by massive underlabeling in single molecule tracking [6].All these methods yield image stacks with very low densities of single molecule signals per frame, which ideally show no overlap.From such images it is possible to estimate the emitter positions to a precision that is mainly limited by the signal to noise ratio of the single molecule images [7].
Crucially, the reliability of the final superresolution images or acquired tracking data critically depends on the quality of the positions obtained from the fitting procedure.In the simplest case, single molecule signals are fitted using a Gaussian function [8].For this, it is assumed that a fluorophore's emission yields an isotropic point spread function (PSF).However, various factors such as microscope aberrations and fluorophore characteristics alter the shape of the PSF.The mismatch of the fitted simple model of a Gaussian function and the intricate shape of the true PSF may well lead to biases in the estimated positions of up to tens of nanometers, thus highly distorting the obtained results [9].Hence, it is vital to incorporate a realistic model in the fitting procedure [10].
In response to this challenge, we have developed a comprehensive and accessible computer application that allows for interactive simulation and visualization of realistic PSFs under a wide range of microscope and fluorophore parameters, enabling an accurate representation of various imaging conditions.In contrast to similar previous tools that allow for calculating and displaying PSFs [11], we implement the full vectorial PSF model [12][13][14] and allow to set a fixed fluorophore dipole orientation; the application is hence not limited to PSFs of isotropic emitters.
In practice, fluorophores are dipole emitters, and their rotation is often restricted by steric hindrances [15].Recently, researchers became interested in performing SMLM under cryogenic temperatures [16,17].Under these conditions, the dipole orientation of a fluorophore is fixed, which leads to generally anisotropic emission patterns and altered PSFs [18].In addition to facilitating the simulation of isotropic PSFs, our tool therefore allows to particularly simulate and visualize the PSF of non-rotating anisotropic emitters with arbitrary dipole orientations.
The best precision that can be achieved in localization techniques depends on the signal to noise ratio, the PSF shape, and the chosen fitting procedure.Ultimately, the localization precision is limited by the Cramér-Rao bound (CRB) [19], which is a theoretical limit for the precision any unbiased estimator can possibly achieve under the given conditions.Notably, the CRB depends on the shape of the PSF.This fact can be taken advantage of by shaping the PSF via manipulations in the back focal plane, often referred to as PSF engineering [1,20].First, this allows to shape the PSF in a way to achieve best lateral position estimate.Second, this breaks the ground for estimating not only lateral position, but also encoding additional parameters including axial position and dipole orientation in the PSF [15].To easily assess the effect of different PSF engineering approaches, we included a feature for adding manipulations of the emitted light in the back focal plane in our app.Eventually, we provide the option to directly determine the CRB for any given PSF shape.
In practice, the shape of the PSF is often affected and degraded by various types of aberrations, for example imperfections of the optical setup such as coma, spherical aberrations and astigmatism.These aberrations negatively affect the fitting procedure and lead to decreased quality of the position estimates and all other estimated parameters.In our app, arbitrary aberrations can be easily included using Zernike coefficients.Moreover, we provide an extended feature that allows for fitting a specific microscope's aberrations from data recorded from a calibration sample.The obtained coefficients can be directly loaded back into the PSF simulation app, yielding simulation results tailored to the user's setup.
As a final feature, we implemented the option to simulate multiple fluorophores in the same region of interest.Overlapping PSFs often occur in cryoSMLM applications, as switching of fluorophores between bright and dark states is decelerated under these conditions [17,21].In addition, overlapping PSFs occur in step-wise photobleaching methods [22], where the signals of multiple fluorophores initially overlap, complicating the localization procedure.
Our application allows to visualize and examine the single fluorophore images, as they would be obtained for specific conditions.It further allows to export the obtained images for further analysis.It could be used, e.g., to provide ground truth information for more advanced fitting procedures.
It supports both the import and export of data, encompassing input parameters and the calculated PSF images.This facilitates a streamlined and accurate approach for the assessment of PSFs under various conditions and PSF engineering approaches.

Features
Our application is structured into several subwindows (Fig. 1).The main window allows to set all simulation parameters and visualization options.Fig. 2 shows a schematic overview of the implemented setup and all parameters that can be varied.We simulate the emission of a fluorophore as a dipole emitter and its PSF as observed in a conventional wide-field microscope setup [14].Here, we will present a short overview of the most important features and illustrative examples of obtained PSFs.For an exhaustive documentation detailing all parameter settings and screenshots of all tabs and subwindows of the app we refer to the app manual in the Supplementary Material.

Fluorophore and microscope parameters
First, the app allows to configure the sample and microscope.The sample is assumed to be a fluorophore, modeled as dipole emitter [14].The fluorophore parameters that can be set include the emission wavelength, its lateral and axial position, and the number of observed photons.Further, the fluorophore can be assumed to be either freely rotating, yielding an isotropic PSF, or fixed.In the latter case, the dipole orientation of the fluorophore can be specified via its inclination angle θ and azimuthal angle ϕ.In addition, the app allows to account for reduced excitation due to dipole inclination by automatically reducing the number of observed photons.For simplicity, the excitation dipole is assumed to be aligned with the emission dipole here.Further, the actual number of photons in the pixels of the image can be subjected to photon shot noise.
The configurable parameters for the microscope include settings of the tube lens (focal length) and objective (numerical aperture, magnification, focal length).The focus position of the objective is set to the transition between objective and intermediate layer.The focus position can be adjusted by changing the defocus value.In addition, an intermediate layer can be added, representing, e.g., a layer on the coverslip with a different refractive index than the immersion medium.For the camera, the pixel size can be specified either by specifying the physical camera pixel size or the pixel size in object space.Further, background noise can be added to the image.The background noise is assumed to follow a Poisson distribution; the parameter for the background noise specifies its standard deviation.

Multiple fluorophores
Ideally, the signals of active emitters in SMLM are well separated.However, this might not always be the case, in particular for high-density SMLM [23] and stepwise photobleaching methods [22], where the PSFs of several fluorophores may overlap.In addition to single emitters, we allow to configure multiple fluorophores with different positions (lateral and axial) and dipole orientations in the same image.The resulting PSF is the superposition of all individual PSFs and can show a distinctively different shape than the individual PSFs as shown in Fig. 4.

Aberrations, transmission and phase retrieval
Up to now, we have assumed an ideal PSF only affected by noise.However, imperfections in the optical path or inhomogeneous refractive indices in the sample can lead to aberrated PSF shapes [24].We allow to include such aberrations in the simulation by introducing a phase shift in the back focal plane that is expanded into a linear combination of Zernike polynomials (see Eq. 5).The most common aberrations, including spherical aberra-  tions, astigmatism and coma can be directly selected and their coefficients adjusted.Alternatively, an arbitrary set of Zernike modes and corresponding coefficients can be specified.Fig. 5 shows an illustrative example how aberrations can affect the PSF shape.
In addition, the PSF shape can be affected by apodization, i.e., nonhomogeneous transmission of the emitted light through the objective.In particular, towards the outer rim of the objective, light transmission is re-duced [25].In order to model this attenuation, the app allows to load a custom transmission mask.
As an additional feature, we provide a subroutine that allows for retrieving the aberrations and transmission of a specific setup from experimental data.For this, a stack of images at various defocus positions from a calibration sample (a fluorescent bead) is required.For details on the recording of the calibration data see 4.4.An illustrative example of experimental data and the fitted model is shown in Fig. 6.

Phase masks
In contrast to undesired aberrations, phase manipulations can be used deliberately in PSF engineering approaches.Here, phase shifts are exploited for shaping the PSF in a way that allows to encode more information.For example, the double helix PSF has been used to allow for determination of the axial position [26], and a vortex phase mask has been shown to allow for retrieving information about the lateral and axial position, as well as the emitter's dipole orientation [15].Any such phase manipulation can be introduced by adding an additional phase factor in the back focal plane (see Eq. 4).Our app offers the feature to select from a variety of commonly used phase masks, including the vortex, double helix, and pyramid phase masks.In addition, a custom phase mask can be loaded.Further, the selected phase masks can be altered by cutting out an inner disk or rotation of the phase mask.A selection of phase masks and their influence on the PSF shape is given in Fig. 7.

Visualization options
Our app allows to visualize the calculated PSF for a set of given parameters in several ways.The default visualization option is a 2D lateral view of the PSF.In addition to the full PSF, the emission can be split into x-and y-polarization channels that can be viewed separately.Further, a full 3D model of the PSF can be calculated.For visualization, we show an xz-projection along with an isosurface plot.The value of the isosurface can be adjusted to show various isosurfaces of the 3D PSF.
The visualization of the plots can be adjusted in several ways as depicted in Fig. 9. First, the simulated region of interest can be specified by setting either the side length of the desired region of interest or the number of pixels per lateral axis (panel a).Second, in case of high background noise, adjusting the contrast of the image may help to better discern the PSF shape (panel b).Third, the colormap of the images can be set by selecting from a few options including the viridis, parula, hot and gray colormaps (panel c).

Import and export options
As noted in the previous subsections, our application allows to import custom data for aberrations, phase masks and transmission.In the subwindow for fitting of the PSF model to experimental data, the fitting results (including aberrations and transmission) can be exported.The saved aberrations and transmission can then be imported into the main window via the aberration and transmission tabs.This allows to incorporate the fit results directly into the simulated model.The calculated simulation data for 2D and 3D PSFs can be exported under the tab Options for further processing or benchmarking of fitting algorithms.

Cramér-Rao bound
Finally, the theoretically achievable localization precision for lateral and axial position can be directly calculated in our app under the tab Options.We calculate the CRB as the diagonal elements of the inverse Fisher information matrix (see Eq. 6).The CRB is affected by the shape of the PSF, the number of observed photons and the background noise.Note that in some cases the PSF does not provide enough information to estimate a particular parameter.For example, the axial position of a molecule that is in focus cannot be accurately estimated from an isotropic PSF in the absence of PSF engineering.In this case, numerical errors lead to unstable values of the CRB.Thus, we do not provide an exact value for the CRB if the calculated number is greater than 1000 nm, i.e., much greater than the diffraction limit of around 200 nm.

Conclusion and outlook
We have presented an easy-to-use MATLAB application that enables the simulation of point spread functions as they appear in SMLM applications.Reflecting the variability in experimental setups, the application offers a wide range of customizable parameters.This allows to tailor the simulation to specific microscope setups, ensuring that the simulated data closely aligns with the experimental reality.The key features of the application include • the simulation of fixed and rotating dipole emitters • the simulation of aberrations modelled by Zernike polynomials in the objective pupil • the addition of optical elements such as phase masks and polarizers in the emission path • retrieval of Zernike aberrations from an experimentally recorded PSF stack • the calculation of the Cramér Rao Bound to assess the theoretically achievable localization precision under specific conditions • flexible visualization options.
The results can be exported in various formats, allowing users to easily generate simulated data for an array of purposes.Most of the parameters can be adjusted via sliders or numerical input fields with a near real-time calculation and visualization of the PSF.While we strived to offer computational efficiency for this first release of our app, a further speed up of calculations, in particular for 3D PSFs and small pixel size, can be achieved in future releases by leveraging GPU capabilities.The features of this app and its interactivity allows users to observe the effect of specific parameters on the shape of the PSF, which, combined with the Cramér-Rao bound may assist, e.g., in designing novel PSF engineering approaches.

Mathematical model
We start by stating the PSF model for a fixed dipole emitter situated on the optical axis in an aberration-free optical system as illustrated in Fig. 2. We assume that the emitter is embedded in a sample medium (e.g., water), followed by an intermediate layer (e.g., a transparent coating) and an objective layer (e.g., immersion oil).The interfaces are assumed to be planes orthogonal to the optical axis.
The angles (θ, ϕ) represent the inclination and azimuth angles characterizing the dipole orientation.We denote by λ the emission wavelength in vacuum and by ⃗ x b ∈ R 2 and ⃗ x f ∈ R 2 coordinates in the back focal plane and image plane, respectively.
The dipole emission pattern, interaction with the different layers of medium and subsequent passage through the infinity-corrected objective is comprehensively described in [14].The starting point of our model is the electric field vector E BFP = E BFP (⃗ x b ; θ, ϕ), expressed in Cartesian coordinates.This field in the back focal plane is defined by [14,Eq. (18)].
A tube lens with focal length f is positioned between the objective and the camera to produce the image.The passage of the unaberrated field through this tube lens is modelled by the Fourier transform, Integration is performed over the circular pupil area.The intensity distribution in the focal plane of the tube lens is given by the absolute value of the electric field, (2) The coordinate system in the back focal plane can in principle be chosen arbitrarily.If one uses emission polarizers, the directions of ⃗ x b are chosen to lie along the orthogonal directions of those polarizers.Then the field resulting from a linear polarizer oriented to transmit x-polarization (y-polarization) only is given by the first (second) component of E(⃗ x f ) [14].

Discretization and noise
The BFP field E BFP is calculated as implementation of [14, Eqs. ( 10)-( 18)] on a discrete grid.To simulate the pixelated grid of the camera, we consider the integrated intensity over the k-th pixel □ k , We call the amount of supporting points used in each dimension for the calculation of (3) the oversampling factor.
The isotropic PSF resulting from a freely rotating emitter is modelled as the superposition of three fluorophores with pairwise orthogonal orientations.
The intensity pattern resulting from multiple nearby emitters is calculated as the superposition of the respective individual intensities.
Camera shot noise is modelled via the realization of a Poissonian random variable, with the calculated noise-free PSF as the mean.A constant background may be added before the application of the Poissonian random variable to model the fluorescence emission of a homogeneous background.

Aberrations and phase shifts
The infinity-corrected optical system provides a space between the objective and the tube lens, where additional optical components such as phase plates can be placed.Any wavefront deformation in this space, either caused by aberration or deliberate distortion, can be modeled by introducing additional phase factors.
We split the phase into two separate parts, φ = φ z + φ r , where φ z is a Zernike term and φ r models the addition of further optical elements, e.g., phase masks in the back focal plane, see Section 2. We expand φ z into a linear combination of orthonormal Zernike polynomials, i.e., where Z i , denotes the i-th Zernike polynomial (using Noll's indices [27]) and w i is the corresponding Zernike coefficient.In particular, we use tip and tilt to model the PSF of an emitter that is laterally displaced from the optical axis.Setting w 1,2 = 1λ produces a lateral shift of 2 N A λ in horizontal and vertical direction, respectively.

Cramér-Rao Bound
We calculate the Cramér Rao Bound (CRB), which is a tool from estimation theory that provides a benchmark for the achievable precision of an estimator [19].The CRB is given by the diagonal elements of the inverse Fisher information matrix of the likelihood function.The likelihood function can be constructed from the forward model (3) and an appropriate noise model.As noise model we choose Poissonian noise, which is a good approximation for camera shot noise.Following [28], we can then calculate the Fisher information as where the summation is over the pixels of the image.The parameter vector ξ denotes the parameters that one wishes to estimate, which is typically just the lateral position.However, they could also include axial position or orientation.The Fisher information matrix is always a symmetric matrix with as many rows as the amount of parameters that are estimated.

Phase retrieval
The PSFs of real microscopes usually differ from the theoretical model.This is in one part explained by the design of the optical elements, in particular the objective lens.Although a modern microscope objective consists of many individual lenses, truly isoplanatic imaging performance cannot be obtained and significant amounts of astigmatism and coma appear at increasing distance from the optical axis.Other effects are known to introduce spherical aberrations, predominantly a refractive index mismatch between the sample buffer solution and objective immersion medium, but also too high or low environmental temperatures or age-related refractive index changes of the immersion oil.Even when these aberrations are small, they can cause systematic errors in the molecule position estimates on the order of several tens to hundreds of microns.In order to avoid these errors, we include aberrations in the model via the Zernike term φ r .The necessary Zernike coefficients are estimated by a phase retrieval algorithm [29] which operates on an experimental 3D image (z-stack ) acquired from a single small fluorescent bead.The algorithm finds Zernike coefficients that minimize the squared L2-Norm of a vectorial error metric ϵ, which is defined as: Here, E and S represent the experimentally recorded and simulated threedimensional bead intensity images and k the voxel index.The quantity γ is a user-definable scalar value between 0 and 1 that influences the fit performance.Smaller values of γ assign increased weight to voxels of lower intensity, e.g., those in out-of-focus planes.We discovered that a γ value of 0.5 was effective in the tested cases, whereas a value of 1 led to inaccuracies as the algorithm stopped at non-ideal local minima.Appropriate z-stacks should be recorded at a maximum possible signal to noise ratio and contain about 10 widefield images covering an axial range from about −1 to 1 µm around the axial bead position.Ideally, the bead is immersed in a mounting medium with a refractive index higher than the NA of the objective.This avoids the formation of a supercritical angle fluorescence (SAF) zone in the objective pupil, and also avoids complications caused by fluorophores located in different regions within the bead emitting different amounts of SAF.The diameter of the bead should not be larger than 200 nm.

Figure 1 :
Figure 1: Overview of app windows.The main app window (top left) comprises several tabs and allows to set all simulation parameters and visualization options.The PSF is visualized in a separate window (top right).The bottom row shows the applied phase mask (left), the aberration (middle) and transmission mask (right).

Figure 2 :
Figure 2: Schematic of implemented microscope setup.(a) On the left, the components of the microscope are listed; on the right, all the parameters that can be varied in the simulations.EBFP and E are the electric fields in the back focal plane of the objective and in the focal plane of the tube lens, respectively.The path of the ray is simplified and intended to show the infinity space between the objective and the tube lens where phase manipulations can be easily modeled.(b) Inclination angle θ and azimuthal angle ϕ defining the dipole orientation

Figure 7 :
Figure 7: Phase masks.The top row shows the phase shift introduced by various phase masks and the bottom row the resulting PSF shape.

Figure 8 :
Figure 8: PSF visualization.The top row shows visualization options in 2D, including the 2D PSF (a) and the split of the emission into polarized channels, (b) and (c).The bottom row shows visualization options for the 3D PSF.(d) xz-projection of the 3D PSF.(e) Isosurface of 3D PSF.Scale bars: 200 nm.