2.1.

Image Formation Process

The MC simulation modeled our current 2D-DWP-based sSMLM system as previously reported.¹³ We used two excitation lasers (532 nm, Exlsrone-532-200, Spectra-Physics; and 647 nm, 2RU-VFL-P-2000-647-B1R, MPB Communications) and an inverted microscope body (Ti2-E, Nikon) equipped with a 100× total-internal reflection fluorescence (TIRF) objective (CFI Apochromat TIRF 100XC Oil, Nikon). Figure 1(a) shows the photon flow in the detection path. The emitted photons from the fluorescent sample pass through a DWP module, consisting of a BS prism (BS010 for R50:T50; and BS049 for R30:T70, Thorlabs), a right-angle prism (RAP) (84-514, Edmund Optics), and a custom-made DWP.¹³ The BS divides the photons into a spatial imaging (zeroth order) path and a spectral imaging (first order) path. The RAP then reflects photons in the spatial imaging path, and the custom-made DWP disperses the photons in the spectral imaging path. Finally, an electron-multiplying charge-coupled device (EMCCD) (iXon Ultra 897, Oxford Instruments) simultaneously detects all of the photons in the zeroth order and first order paths for further processing.

Fig. 1

Illustration of the 2D-DWP image formation and MC simulation steps. (a) Schematic of the 2D-DWP image formation; (b) photon emission simulation process; (c) a simulated image of the photon distribution on the EMCCD following a Gaussian distribution of photon wavelengths; (d) corresponding EMCCD-detected image for photon distribution shown in panel (c); (e) simulated EMCCD background, including background and readout noises; (f) combined image from panels (d) and (e); (g) 1D fitting the spectral peak and with from spectral profile, where the zeroth order image provides the reference needed to obtain the corresponding first order spectral information; (h) measured photon count distribution of AF647 with a 10 ms exposure time; (i) measured bulk emission spectrum of AF647; (j) measured $x$ and $y$ variations of the individual photons due to the diffraction limit, with $x$- and $y$-summed histograms shown; (k) measured spectral deflection characteristics of the DWP module, and (l) measured ratios of first order/zeroth order efficiency in the DWP module when different BSs are used. RAP, right-angle prism; BS, beamsplitter.

2.2.

MC Simulation Process

Figure 1(b) shows our MC simulation process, from the generation of a blinking event by individual fluorophores to image detection. We start from the ground truth of $N$ photons emitted in a single frame from one photoswitching event with parameters $x$, $y$, and $\lambda$, where $N$ photons are either sampled from a lognormal function or assumed to be a constant value, $x$ and $y$ are the spatial positions of each photon, and $\lambda$ is the wavelength of each photon sampled from the emission spectra. Because BSs often have wavelength-dependent splitting ratios and introduce photon losses, we pass the $N$ photons through an efficiency function $\epsilon (\lambda )$ to determine whether each photon goes to the zeroth order path or first order path or is lost. Photons in the zeroth order path will be detected directly by the EMCCD, whereas photons in the first order path will pass through a dispersion function $f(\lambda )$ to determine their wavelength-dependent shift along the $x$-axis ($x_{\mathrm{disp}}$), mimicking the behavior of the DWP. Next, we combine the locations of the photons, $X=x+x_{\mathrm{disp}}$ and $Y=y$, to determine the exact locations where each photon will be incident on the EMCCD array in Fig. 1(c). This signal is then discretized by individual EMCCD elements due to the finite pixel size ($a$) with varying intensities to give a signal ($\mathbf{S}$) illustrated pictorially in Fig. 1(d). We further add a background signal consisting of Poisson-distributed background noises ($\mathbf{B}$) and Gaussian-distributed readout noises²¹ ($\mathbf{R}$), shown in Fig. 1(e). By summing the digitized signals from Figs. 1(d)–1(e), we apply signal amplification using a Gamma distribution, which is known to effectively model the behavior of an EMCCD²¹. This process yields $\mathbf{Z}$, the simulated sSMLM frame shown in Fig. 1(f).

To obtain the subpixel peak location from the zeroth order spatial image, we fit the zeroth order image with a bivariate Gaussian function using the maximum likelihood estimation (MLE) method,¹¹ giving us an estimate of the expected values of the peak locations $x$ and $y$ as $\hat{x}$ and $\widehat{y}$, respectively. To obtain the spectral peak from the first order spectral image, we first calibrate the $x$-axis of the spectral graph with a spectral calibration procedure (detailed in Sec. 2.4), which allows us to calculate the wavelength values along the $x$-axis for each zeroth order localization. We then fit the spectrally calibrated image with either a univariate or bivariate Gaussian function to estimate the spectral peak (${\widehat{\mu }}_{\lambda }$) and spectral width (${\widehat{\sigma }}_{\lambda }$) together with the spatial localization. The procedure is illustrated with a one-dimensional (1D) spectral graph in Fig. 1(g), where we integrated the signal intensity along the $y$-axis for clarity. An alternative method commonly used for spectral characterization is the spectral centroid method, which takes the weighted average of the different wavelengths detected to calculate a spectral centroid associated to each localization.^12,13,22 In this study, we fitted the resulting image with a bivariate Gaussian function. This approach enhances noise rejection by taking advantage of the larger number of samples available in 2D data, as opposed to summing it along the $y$-axis in 1D. As the 2D Gaussian fit utilizes information from both dimensions, it provides a more accurate representation of the underlying distribution and consequently improves the robustness against noise.

To illustrate the functions $N$, $\lambda$, $f(\lambda )$, and $\epsilon (\lambda )$ in Table 1 used in the simulation, we extracted the distribution of photon count and spectral signature of the commonly used SMLM dye Alexa Fluor 647 (AF647). Figure 1(h) shows the lognormal fitted distribution of photon counts $N$ of AF647 extracted using the process described in Sec. 2.3.1. Figure 1(i) shows the spectral distribution $\lambda$ of AF647 obtained with the procedure described in Sec. 2.3.2, and Fig. 1(j) shows a simulated point spread function (PSF) using a 2D-Gaussian function, which is a close approximation to the Airy disc.²³ Figure 1(k) shows the spectral dispersion function, $f(\lambda )$, experimentally measured using the procedure described in Sec. 2.3.3. Figure 1(l) shows the efficiency function, $\epsilon (\lambda )$, for BSs with reflectance (R):transmittance (T) of R30:T70 and R50:T50 experimentally measured using the procedure described in Sec. 2.3.4. Finally, we repeated the process 20,000 times to obtain the spatial and spectral precisions, defined as the standard deviation of $\hat{x}$ and ${\widehat{\mu }}_{\lambda }$, respectively. We summarize the parameters used in the simulation in Table 1.

Table 1

Summary of simulation parameters.

2.3.

Extraction of Parameters from Experimental Data

2.4.

Post-Processing of EMCCD Images

2.4.2.

Correcting for nonlinearity in spectral fits in the DWP module

To accurately fit the spectra obtained from the DWP module, it is crucial to consider the inherent nonlinearity in spectral dispersion. As shown in Fig. 2(a), a hypothetical light source with a uniform spectral power density undergoes spectral dispersion after passing through the DWP module, resulting in non-uniform intensity across the image sensor, as shown in Fig. 2(b). This effect is particularly noticeable in regions with lower spectral dispersion, where the incident power is greater. Consequently, fitting a Gaussian curve to the spectrum would result in a bias toward the longer wavelengths due to the higher intensity in the region. To account for this spectral bias, we normalized the spectra with the spectral dispersion, as shown in Fig. 2(c), which results in a corrected image with uniform intensity throughout the image sensor. This normalization process corrected the nonlinear spectral dispersion in the DWP module and ensured accurate spectral fitting.

Fig. 2

Illustration of the correction process of nonlinearity in the spectral dispersion with DWP. (a) Illustration of the wavelength distribution of an ideal uniformly distributed spectral density light source; (b) power density detected by the EMCCD after the DWP; and (c) corrected power density.

2.5.

Validation with Fluorescence Microspheres

To validate our MC model, we imaged microspheres with emission peaks at 560, 580, 605, 645, and 720 nm (F8800, F8794, F8801, F8806, and T8870, Invitrogen) deposited on #1.5-thick glass coverslips (12-541-B, Fisher Scientific). The other physical properties of these microspheres are detailed in Table S2 in the Supplementary Material. We first washed the glass coverslips with phosphate-buffered saline (PBS) and coated the coverslips with a 0.001% poly-L-lysine solution for 5 min. We then washed the coverslips with PBS and stained the glass coverslips with a diluted sample of each microsphere for 1 h. Next, we washed off the unbounded microspheres with PBS three times and added a single drop of antifade buffer (P36961, Invitrogen). Finally, we sealed the coverslips with black nail polish (114-8, Ted Pella) and left them overnight to dry.

We first imaged the microspheres without the DWP module to obtain the photon count of each nanosphere and then imaged the same microspheres with the DWP module for 100 frames. To obtain the spatial localization precision, we processed the 100 frames of microspheres using ThunderSTORM²⁴ and treated the standard deviation from $\hat{x}$ calculated from the 100 frames as the localization precision. To obtain the spectral precision, we identified a region of interest from 450 to 800 nm relative to the zeroth order image after spectral calibration using the method previously described.^5,13,22 This region of interest gives us a spectral profile similar to the one seen in Fig. 1(f), which we use to analyze the spectral signature of the fluorophore further. Next, we subtracted the background using averaged neighboring windows of the first order image. Then, we performed a least-squares fit of the spectral profile with a bivariate Gaussian function to identify the spectrum corresponding to the zeroth order blinking.

3.1.

Validating MC Simulation with Theoretical Analyses

To validate our MC model, we first benchmarked the precisions generated from our model against theoretical models of spatial localization precision $\mathrm{\Delta }x$ defined as¹¹

Equation (5) describes the performance limit of the MLE when fitting the PSF into an approximate 2D Gaussian curve, also referred to as the Cramér–Rao lower bound.¹¹ Equation (6) assumes (1) a uniform spectral dispersion, (2) a Gaussian spectral shape with a spectral peak of ${\mu }_{\lambda }$ and a spectral FWHM of ${\sigma }_{\lambda }$, and (3) spectral centroid in the calculation of the spectral precision.¹² The spectral centroid method takes the weighted average of the intensities captured in the pixels associated with different wavelengths, yielding a weighted representation of the spectra recorded.

To provide a fair comparison to validate our simulation, we matched the simulation based on the assumptions of Eq. (6) listed above. We computed the spatial localization precision, defined as the standard deviation of the localized spot from the true location, and the spectral precision, defined as the standard deviation of the spectral centroid or the fitted spectral peak. We performed MC simulations under the following conditions: $N=700$, $b_{\mathrm{bg}}=15$, spectral $\text{peak}=680\mathrm{nm}$, and spectral $\mathrm{FWHM}=41\mathrm{nm}$, using a linear spectral dispersion of 5 nm/px, unless otherwise stated.

Figure 3 shows the influence of photon counts $N$, background photons $b_{\mathrm{bg}}$, and emission spectral peak ${\mu }_{\lambda }$ on lateral localization and spectral precisions. In Figs. 3(a)–3(c), we compared the theoretical localization precisions calculated with Eq. (5) with our MC simulation, with mean absolute errors of 2.5, 0.8, and 0.3 nm relative to the theoretical values, respectively. Similarly, in Figs. 3(d)–3(g), we compared the theoretical spectral precisions calculated with Eq. (6) with our MC simulation using the spectral centroid method, with mean absolute errors of 2.1, 7.6, and 2.2 nm, respectively. For the spectral fitting method, the simulated precisions are consistently better, with mean absolute errors of 6.3, 4.6, and 3.6 nm for Figs. 3(d)–3(g), respectively.

Fig. 3

Comparison of localization precisions between theoretically predicted results (dashed line) and simulation results (dots) with respect to (a) photon count, (b) background photons, and (c) emission spectral peak. Comparison of spectral precisions between theoretically predicted results (dashed lines), simulated results calculated with spectral centroids (crosses), and simulated results calculated with a bivariate Gaussian fit (dots) with respect to (d) photon count, (e) background photons, and (f) spectral peak.

In Figs. 3(a) and 3(d), we observe that increasing the number of emitted photons from 1500 to 5000 improved the localization precision from 8.1 to 4.0 nm. Similarly, the spectral precision improved from 4.4 to 2.0 nm with the spectral centroid method and 3.2 to 1.6 nm with the spectral fitting method. With a limited photon budget, splitting the photon budget between zeroth order and first order images can affect the tradeoff between the localization precision and spectral precision.

In Figs. 3(b) and 3(e), increasing background photons from 0 to 30 per pixel worsens the localization precision from 10.6 to 15.0 nm. Similarly, spectral precision worsens from 4.4 to 13.6 nm with the spectral centroid method and from 2.7 to 4.3 nm with spectral fitting. The spectral fitting method is more robust against noise because it estimates the background term by fitting for and averaging over all pixels available in the spectral window, reducing the variation in the estimated spectral peak ${\hat{\mu }}_{\lambda }$. By contrast, the spectral centroid method takes the weighted average of the spectral window, which is more susceptible to noise variations. In addition, accurately estimating the background level can be challenging,²⁸ which can result in biases to the spectral centroid value. A positive bias can result from the presence of background contribution to the spectral window, which can shift the spectral centroid value toward the mean of the spectral window and reduce the spectral precision. Hence, we prefer the bivariate Gaussian fitting method due to its robustness in estimating the spectral peak ${\hat{\mu }}_{\lambda }$ and calculating for spectral precisions.

In Figs. 3(c) and 3(f), we observe that, as the center wavelength of the emission peak increases from 500 to 740 nm, the localization precision deteriorates from 9.6 to 15.0 nm. Simultaneously, the spectral precision is stable at $\sim 9.2\mathrm{nm}$ with the spectral centroid method and at 3.5 nm with the spectral fitting method. This effect is more pronounced in Fig. 3(c), where the wider FWHM of the PSF results in improved localization precision. However, in Fig. 3(f), a spectral emission FWHM of 47 nm plays a more significant role in determining the overall width of the first order image with added spectral dispersion. Overall, our MC model agrees well with the theoretical predictions, with a worse-case average absolute error of 2.5 nm for lateral precision and 2.1 nm for the spectral precision with the spectral centroid method across all cases.

3.2.

MC Simulation of Different BSs

After validation, we used the MC model to predict the performance of different BSs used in the DWP module. Figure 4 shows how photon counts $N$, background photons $b_{\mathrm{bg}}$, and emission spectral peak ${\mu }_{\lambda }$ affect spatial localization and spectral precisions with different BSs. We conducted MC simulations under the following conditions: $N=700$, $b_{\mathrm{bg}}=15$, ${\mu }_{\lambda }=680\mathrm{nm}$, and emission spectral $\mathrm{FWHM}=41\mathrm{nm}$, and we used the DWP system with different BS ratios unless otherwise stated.

Fig. 4

Comparison of localization precisions between an R50:T50 (orange dots) and R30:T70 (blue squares) BS with respect to (a) photons count, (b) background photons, (c) spectral peak, and (d) spectral FWHM. Comparison of spectral precisions between an R50:T50 (orange dots) and R30:T70 (blue squares) BS with respect to (e) photons count, (f) background photons, (g) spectral peak, and (h) spectral FWHM.

Figures 4(a)–4(d) show that using an R30:T70 BS results in a 35% worse localization precision on average compared with an R50:T50 BS. However, Figs. 4(e)–4(h) demonstrate that an R30:T70 BS provides 8% better spectral precision on average compared with an R50:T50 BS. This is because the R30:T70 BS redirects more of the photon budget toward the spectral channel, which leads to lower spectral precisions, as shown in Fig. 3(d). The tradeoff between spectral and spatial localization precisions is especially prominent in photon-limited photoswitching events when $N<1000$, as shown in Figs. 4(a) and 4(e), where we observed up to a 12.3 nm worsening in the localization precision and up to 2.7 nm improvement in the spectral precision for $N=320$. From Fig. 1(h), many photoswitching events of AF647 have $N<1000$, which suggests that the change of BS can play a significant role in optimizing the accuracy in the classification of different fluorophores.

Increasing the number of background photons from 0 to 50 per pixel leads to an increase in the localization and a minimal increase in the spectral precision in both BSs. In Fig. 4(b), the localization precision worsens from 11.4 to 17.8 nm for the R50:T50 BS and from 16.1 to 24.7 nm for the R30:T70 BS. In Fig. 4(f), the spectral precision worsens slightly from 3.2 to 4.5 nm for the R50:T50 BS and from 2.9 to 4.0 nm for the R30:T70 BS. Overall, switching from an R50:T50 BS to an R30:T70 BS led to an average worsening in the localization precision by 5.6 nm and an improvement in the spectral precision by 0.3 nm when the number of background photons per pixel is varied.

Figure 4(c) shows that, as the spectral peak changes from 500 to 750 nm, the localization precision worsens from 8.6 to 11.8 nm for the R50:T50 BS and from 10.9 to 16.6 nm for the R30:T70 BS. This worsening is attributed to the change in PSF’s FWHM, as observed in Fig. 3(c). However, Fig. 4(g) reveals a nonlinear relationship between the spectral precision and spectral peak, unlike in Fig. 3(f), with mean values of 3.0 nm for the R50:T50 BS and 2.8 nm for the R30:T70 BS. This nonlinearity results from the nonlinear spectral dispersion of the DWP system, as shown in Fig. 1(k), causing spectral precision to vary based on the emission wavelength. Overall, switching from an R50:T50 BS to an R30:T70 BS led to an average worsening in the localization precision of 3.5 nm, an average improvement in the spectral precision of 0.2 nm, and up to a 0.4 nm improvement in the spectral precision at a spectral peak of $\sim 540\mathrm{nm}$.

Figure 4(d) shows that the localization precision remains relatively stable when the FWHM spectral bandwidth increases from 10 to 80 nm. Figure 4(h) shows that the spectral precision worsens from 2.3 to 6.4 nm for the R50:T50 BS and from 2.3 to 5.5 nm for the R30:T70 BS. This increase is due to the broader spread of the first order image on the EMCCD, with results in a lower signal-to-noise ratio in the detected spectra for a fixed photon count $N$, and, consequently, poorer precision. However, the localization precision remains mostly unaffected as the PSF does not change significantly compared with the spectral shape. Overall, switching from an R50:T50 BS to an R30:T70 BS led to an average worsening in the localization precision by 5.2 nm and an average improvement of the spectral precision by 0.4 nm.

3.3.

Testing MC Simulation Against Experimental Data

Figure 5 tests our MC simulation against experimental results obtained by imaging microspheres. The results demonstrate a good match between the MC-DWP simulation and our experimental data for the localization precision, as shown in Figs. 5(a)–5(r), with an average absolute error of 2.7 nm for the R30:T70 BS and 1.5 nm for the R50:T50 BS. However, for the spectral precision shown in Figs. 5(c)–5(t), we observed a slightly higher error, with an average of 0.9 nm for the R30:T70 BS and 1.0 nm for the R50:T50 BS, compared with the simulation.

Fig. 5

Localization and spectral precisions of different microspheres under varying numbers of emitted photons with DWP-based sSMLM system: (a)–(r) localization precisions and (c)–(t) spectral precisions. Each row shows the localization and spectral precisions of a single species of microspheres. The BS used in each case is indicated at the top of each column.

The difference in the spectral precision can be attributed to two factors: (1) spectral variations within individual microspheres and (2) imperfect background estimations in the experimental results. Although the simulation results used the bulk spectrum of the microspheres to generate the data, the experimental results measured the spectra of individual microspheres. At the molecular level, individual spectra are not readily available, and it is known that, when adding two unequal Gaussian probability density functions (PDFs), the result is a PDF with a broader width. Because the individual spectral signatures of the microspheres differ from the bulk phase due to heterogeneity,²² their narrower widths likely collectively combine to produce a broader bulk phase spectral signature. As shown in Fig. 4(h), increasing the spectral FWHM by 10 nm results in $\sim 10\%$ degradation in the spectral precision. The influence of spectral heterogeneity present in our experimental findings explains the slightly improved calculated spectral precision in the experiments compared with the simulation results.

Furthermore, background estimation and subtraction can be challenging and biased because the microspheres are not photoswitching. During data processing, background subtraction may have caused some microspheres to have consistently lower spectral precisions than the simulated values, resulting in systematic errors in the calculated spectral precision.²⁸ Nonetheless, our MC-DWP simulation provides a reasonable estimate of the performance of the different nanospheres.

Finally, our analysis does not consider changing the spectral dispersion of the system because the DWP module was designed using an optimal spectral dispersion range of 4 to 6, which minimizes the spectral uncertainty when $N$ ranges from 1000 to 2000.¹² This $N$ range is the photon budget available for most fluorophores commonly used in biological imaging applications. Therefore, the beamsplitting ratio and fluorophore combinations are key to maximizing system performance. In particular applications, such as the spectroscopic analysis of dyes,²² in which the photon budget is fixed, we may omit the BS and use the entire available photon budget for spectroscopic analysis. In this case, the only way to tune the spectral performance of the system is to change the spectral dispersion.