## Abstract

Superresolution optical fluctuation imaging (SOFI) is a simple and affordable super-resolution imaging technique, and attracted a growing community over the past decade. However, the theoretical resolution enhancement of high order SOFI is still not fulfilled. In this study, we identify “cusp artifacts” in high order SOFI images, and show that the high-order cumulants, odd-order moments and balanced-cumulants (bSOFI) are highly vulnerable to cusp artifacts. Our study provides guidelines for developing and screening for fluorescence probes, and improving data acquisition for SOFI. The new insight is important to inspire positive utilization of the cusp artifacts.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Superresolution (SR) imaging techniques such as stimulated emission depletion microscopy (STED) [1], photo activated localization microscopy (PALM) [2], structured illumination microscopy (SIM) [3], stochastic optical reconstruction microscopy (STORM) [4] and their many derivatives have gained prominence in recent years [5–8] by providing imaging below the diffraction limit of light. Superresolution optical fluctuation imaging (SOFI) [9] is an affordable alternative to these methods. In SOFI, consecutive frames are acquired to form a movie of the imaging sample, which is labeled with stochastically blinking probes. The auto- and cross-correlations of the time trajectories of the pixel intensities are then calculated and subsequently used to construct the different-order cumulants in order to obtain high-order SOFI images. Since SOFI does not require any special hardware and is based on a simple mathematical algorithm, it has the potential to “democratize” SR imaging. The only requirement for SOFI is that the fluorescence probes used should exhibit stochastic blinking at a rate that can be captured by a camera. Quantum dots (QDs) [9], organic fluorophores (dyes) [10], fluorescence proteins [11,12], carbon nanodots [13], and Raman probes coupled to plasmonic nanoparticles [14] have all been used for SOFI. Other forms of optical fluctuations have also been exploited for SR imaging using SOFI, such as those related to diffusion-assisted Forster resonance energy transfer [15], protein-protein interactions [16], and the diffusion of nonblinking probes [17]. The large variety of probes available and the various implementations of SOFI suggest that it may be useful in a variety of applications.

The resolution enhancement of SOFI is manifested by the reduced width of the point spread function (PSF) in the reconstructed SOFI image. Theoretically, the PSF width for a *n ^{th}*-order SOFI image is reduced by a factor of

*n*as compared to that of the PSF in the original acquisition. When deconvolution or Fourier reweighting [18] is combined with the estimation of the system’s optical PSF [18], an additional enhancement of

^{1/2}*n*can be realized, bringing the total theoretical resolution enhancement factor to

^{1/2}*n*. Moreover, in principle, the resolution could be improved even further by increasing the SOFI order,

*n*.

However, in practice, the resolution enhancement of high-order SOFI is limited by two fundamental issues. The first issue is the nonlinear expansion of the dynamic range of the pixel intensities in high-order (order > 2) SOFI images. This increase in the dynamic range renders the fine details in SOFI images imperceptible. The balanced SOFI (bSOFI) [19] method has been developed and adopted widely [20] for solving the first problem. This method involves calculating the balanced cumulants, which compensate for the expanded dynamic range of the pixel intensities in high-order SOFI images. However, the application of bSOFI in the case of high-order cumulants (> 4) has rarely been reported and generally results in artifacts in practice. The second issue, which was mostly ignored in the past but was identified in this study, is attributable to the positive and negative oscillations of the cumulants (Fig. 1). More specifically, the boundaries between the negative and positive regions result in artifacts (we have labeled these as “cusp artifacts”). We found that cusp artifacts are intrinsic to high-order SOFI cumulants due to finite and inhomogeneous blinking profiles with limited signal to noise ratio. Subsequent post-processing steps that rely on the high-order SOFI cumulants would probably be adversely affected by these artifacts, such as deconvolution algorithms with positivity constraints (e.g., MATLAB’s “deconvlucy” and “deconvblind” functions), and high-order bSOFI reconstructions, which assumes perfect deconvolution.

We identified cusp artifacts in experimental data, and analyzed them both theoretically and through simulations of various conditions. Our study shows that the cusp artifacts are originated from the nature of a mixture of positive and negative “virtual emitters” (explained in section 3) in the high order cumulant image, and the insights gained regarding such mixed negative and positive values can serve as guidelines for avoiding the cusp artifacts, and inspire new applications.

Concerning avoiding the cusp artifacts, our study suggests emitters with uniform photophysical parameters (blinking and bleaching) should be used; long movies should be acquired to ensure sufficient statistics in the acquisition, and the blinking statistics of the fluorophores should be tuned to be free of cusp-artifacts (see Fig. 1). Besides, bleaching should be minimized and bleaching correction [14] should be performed on the data set. Cusp artifacts are more pronounced for high-order cumulant reconstructions, where the expansion of the dynamic range degrades the quality of the reconstructed image. We used a cusp-artifact-independent method to compress the dynamic range of the pixel intensities (labeled as *ldrc*; the details are provided in [21]). We compared the obtained results with those of bSOFI reconstructions. The analyses of the theory, simulations and the experimental data all showed that the bSOFI algorithm fails to faithfully reconstruct the true image at high orders when influenced by cusp artifacts. We also examined and compared a moments-reconstruction approach to cumulants-reconstruction approach. While more details on moments reconstruction are provided in [21], we show here that even-order moments are immune to cusp artifacts and tend to smoothen the features-of-interest, while still allowing for some degree of resolution enhancement (over the diffraction limit, with the limit being 2^{1/2}-fold for the *n*^{th}-order moment and an additional *n*^{1/2}-fold when deconvolution is performed, resulting in an overall theoretical resolution enhancement of (2*n*)^{1/2}). Therefore, the proposed method has significant potential for use as a practical (but mathematically nonrigorous) solution for avoiding cusp artifacts. Compared to the case for bSOFI, the even-order moments combined with *ldrc* yield significantly more faithful images up to the 6^{th} order. (Details are presented in the accompanying manuscript [21]). We showed both theoretically and experimentally that cusp artifacts could be avoided by using even-order-moment reconstruction [21]. With the provided insights about the nature of mixed positive and negative values in the cumulant image, new methods could be developed to decipher the underlying physics through the statistics revealed by the nature of the cumulants by combining multiple orders of cumulants and solve the underlying blinking statistics altogether as a global inverse problem.

The rest of the manuscript is organized as follows: in Section 2, we briefly review the underlying theory of SOFI. In Section 3, we introduce the mathematical concept of “virtual emitters” and “virtual PSF” for high-order SOFI images (to be used in the subsequent sections). In Section 4, we present a theoretical explanation of cusp artifacts. In Section 5, we examine the conditions that result in cusp artifacts. Next, in Section 6, we evaluate the adverse effects of cusp artifacts on balanced cumulants and post-processing deconvolution algorithms. Further, we show that cusp artifacts can be eliminated completely by utilizing even-order moments (instead of cumulants) for image reconstruction. In Section 7, we compare the performances of the various algorithms using real data. Finally, we conclude the manuscript by discussing the implications of our findings in Sections 8.

## 2. Review of SOFI theory

A brief review of SOFI theory is given below. For SOFI reconstruction, a stack of frames (a movie) is obtained using a simple wide-field imaging system. The sample is labeled with stochastically blinking probes. Each point emitter (probe) in the sample plane is imaged onto the camera plane via the optical imaging system. Further, owing to the diffraction limit of light, the intensity distribution of imaging system takes the shape of the PSF. The signal captured at a given camera pixel located at $\vec{r}$ can be expressed as follows (excluding the binning effects due to pixilation):

where $\vec{r}$ is the location of the pixel in the imaging plane,*N*is the total number of emitters,

*k*is the emitter index, ${\vec{r}_k}$ is the location of the

*k*

^{th}emitter, ɛ

_{k}is the brightness of the

*k*

^{th}emitter when it is in the “on” state, and

*b*(

_{k}*t*) is the blinking time trajectory of the

*k*

^{th}emitter.

*b*(

_{k}*t*) = 1 when the emitter is in the “on” state, and

*b*(

_{k}*t*) = 0 when the emitter is in the “off” state. $U(\vec{r})$ is the PSF of the imaging system, which is determined by the optical setup as well as the emission wavelength of the emitters.

In SOFI, the temporal average of each pixel’s time trajectory is subtracted from the signal, such that only the fluctuations (around zero) are considered:

_{t}represents the time-average operator, and δ

*b*(

_{k}*t*) represents the fluctuations in the blinking time trajectory

*b*(

_{k}*t*). The cumulants of $\delta F(\vec{r}, t)$ can then be calculated. In the case of a 2

^{nd}-order cumulant, the cumulant (C

_{2}) is equivalent to the correlation function:

^{nd}-order cumulant of δ

*b*(

_{k}*t*). The derivation above can be extended to higher-order cumulants as well [22]. Notice that the cumulants are additive [22]. Hence, the n

^{th}-order cumulant of $\delta F(\vec{r}, t)$ can be expressed as the sum of the cumulants of the individual emitters:

^{th}-order cumulant of δ

*b*(

_{k}*t*) and can be simplified as ${C_n}[\delta {b_k}(t)]$ when the time lags are not specified.

## 3. The “virtual-emitter” interpretation for high-order SOFI

In order to better facilitate the analysis and discussions about cusp-artifacts, we propose a conceptual metaphor to interpret the physical meaning of high order SOFI cumulants based on its mathematical description (Eq. (2.6)). We note that Eq. (2.1) (which describes the direct observation image) and Eq. (2.6) (which describes the n^{th}-order SOFI cumulant) share a common form. The original image (Eq. (2.1)) is formed by the summation of weighted signal from all the emitters, where the weighting factor $U(\vec{r} - {\vec{r}_k})$ describes the weight determined by the PSF $U(\vec{r})$ and the distance between the pixel location ($\vec{r}$) and the emitter location (${\vec{r}_k}$) in the sample plane, and the term ${{\epsilon}_k}{b_k}(t)$ describes the time dependent brightness profile of the *k*^{th} emitter. Equation (2.6) is of the same form as Eq. (2.1) with two terms replaced: The point spread function term $U(\vec{r} - {\vec{r}_k})$ in Eq. (2.1) is replaced by ${U^n}(\vec{r} - {\vec{r}_k})$ in Eq. (2.6), and the time dependent brightness term δ*b _{k}*(

*t*) in Eq. (2.1) is replaced by ${\epsilon}_k^n{\omega _{n, k}}({\tau _1},{\tau _2},\ldots ,{\tau _{n - 1}})$ in Eq. (2.6). Therefore, Eq. (2.6) can be perceived as describing an image acquired by an optical system with a PSF ${U^n}$.with signals generated by sample with ‘virtual’ emitters located at the same locations of the real emitters, but with emitter brightness described by ${\epsilon}_k^n{\omega _{n, k}}({\tau _1},{\tau _2},\ldots ,{\tau _{n - 1}})$. We address the term ${\epsilon}_k^n{\omega _{n, k}}({\tau _1},{\tau _2},\ldots ,{\tau _{n - 1}})$ in Eq. (2.6) (which corresponds to the brightness term in Eq. (2.1)) as the “virtual brightness”, and the term $U(\vec{r} - {\vec{r}_k})$ (which corresponds to the PSF term in Eq. (2.1)) as the ‘virtual PSF’. Notice that the emitter locations ${\vec{r}_k}$ are not changed, i.e. virtual emitters share the same locations as the real emitters in the sample. Additionally, the virtual PSF is the original PSF raised to the power of

*n*. In cases when the PSF can be estimated by a three-dimensional Gaussian function, the width of the PSF is reduced by a factor of

*n*

^{1/2}, as shown below:

_{xy}and σ

_{z}are the respective widths of the original PSF in the xy- and z-planes. The reduction of the virtual PSF widths provides the basis for resolution enhancement in the SOFI cumulants (in all three dimensions). It is worth noticing that the virtual PSF is the n

^{th}power of the real PSF regardless of the profile of the real PSF.

In summary, a SOFI cumulant image is equivalent to an image captured with a virtual microscope that has a virtual PSF being the nth power of the original PSF (i.e., reduced width and increased resolution with Gaussian PSF estimation), and the captured signal is created by virtual emitters that are located at exactly the same locations as the real emitters, but with altered brightnesses (i.e. virtual brightnesses) that are the cumulant of the blinking profile of each emitter (and can exhibit either positive or negative values).

We will demonstrate in the following sections that ${\omega _{n, k}}({\tau _1},\;{\tau _2},\;\ldots ,{\tau _{n - 1}})$ exhibits negative and positive values due to two main causes: blinking behavior and finite acquisition time.

## 4. Origin of cusp artifacts

For simplicity, we discuss here a noise free condition without bleaching, and assume a two-state blinking profile *b*(*t*) for the real emitters, with *b*(*t*) = 1 indicating the on state and *b*(*t*) = 0 indicating the off state. Further, we define ρ* _{k}* as the fraction of time the emitter spends in the on state during the signal acquisition step (ρ

*is the on-time ratio). Adjacent emitters could have different ρ*

_{k}*values (in principle) owing to the following reasons: (i) finite acquisition time (i.e., not reaching statistical significance to represent the underlying statistical process); (ii) subtle changes in the local microenvironments; and (iii) inhomogeneity in the photophysical properties of the emitters (for example, quantum dots). And the cumulants of*

_{k}*b*

_{k}(

*t*) can be expressed as a function of ρ

*. If the acquisition time is long enough, and if the blinking behavior can be statistically described by a Poisson process (or Bernoulli statistics), ρ*

_{k}*converges to*

_{k}_{k}does not converge [23], or with insufficient acquisition time to ensure statistical significance, or when the blinking behavior follows a different model, different formulas of on-time ratio ρ

_{k}can be deduced. Regardless, 0≤ρ

_{k}≤1 holds for any two-state blinking model. Therefore, for generality, we present the following discussions in the form of functions of ρ

_{k}. If, for simplicity, we set all the time lags equal to zero in the cumulant calculations, we find that the cumulants with orders higher than two will exhibit positive-negative oscillations as a function of ρ

_{k}. The expressions for the first six cumulants (2

^{nd}order to 7

^{th}order) as functions of ρ

_{k}are given as follows (the derivations are given in Appendix 1 [28]):

*k*, and the time lags (τ

_{1}, …, τ

_{n−1}) have been eliminated to simplify the notation.

Figure 1(v) shows cumulants of different orders as functions of the on-time ratio, ρ. The higher-order cumulants (> 2^{nd} order) clearly exhibit positive-negative oscillations as a function of ρ. If the virtual brightnesses of two nearby virtual emitters have opposite signs (Fig. 1(i) and Fig. 1(ii)), the corresponding amplitude cross-section of the SOFI amplitude image (consisting of the convolution of the virtual emitters with the virtual PSF) will exhibit positive and negative lobes (Fig. 1(iii)). The absolute value of this cross-section (determining this is a common step in displaying an image) exhibits a cusp (Fig. 1(iv), arrow). We have therefore labeled such artifacts as “cusp artifacts.”

In order to better demonstrate cusp artifacts, we performed simulations for three adjacent blinking emitters with the same on-state brightness, but with on-time ratios of 0.831, 0.416, and 0.103 respectively (Fig. 2). The simulations show clearly that high-order cumulants can take negative values that coexist with the positive values, leading to cusp artifacts.

Moreover, when displaying gamma-corrected [24] experimental high-order SOFI-processed images (before performing deconvolution or Fourier reweighting), cusp artifacts can be observed readily for orders greater than two, as can be seen from the gray/red boundaries in Fig. 3.

## 5. Exploration of parameter space of cusp artifacts

In order to better understand the prevalence and significance of cusp artifacts in SOFI reconstructions, we performed a series of realistic simulations to explore the relevant parameter space. A simulator was developed that propagates the emissions of point emitters in the sample plane onto the detector (camera) plane using the Gibson-Lanni PSF model [25]. The point emitters’ time blinking trajectories were simulated to blink according to Poisson statistics. Poisson noise was simulated as described in the accompanying manuscript [21]. Actual experimental background noise was recorded using an electron-multiplying charge-coupled device (EMCCD) camera and added to the simulated movie. Next, the SOFI cumulants of up to the 7^{th} order were calculated and analyzed. Long (20,000 frames) simulations were performed to ensure that the blinking trajectories exhibited suitably high statistical significance. Details of the simulator are given in the accompanying manuscript [21] and posted on a public GitHub account as SR_Simu3D [26].

In the first set of simulations (“Simulation-1”), we simulated four different populations of emitters (P1, P2, P3, and P4) with different ρ distributions (Fig. 4(i), dashed red, blue, black, and green curves), whose ranges were as follows: 0.49 ≤ ρ ≤ 0.51 for P1, 0.53 ≤ 0.87 for P2, 0.11 ≤ ρ ≤ 0.93 for P3, and 0.05 ≤ ρ ≤ 0.95 for P4 (Fig. 4(i)). Comparing these distributions to Fig. 1(v) allowed us to predict the signs of the resultant virtual brightnesses (Fig. 4(ii)). A simulated filamentous morphology was then populated with emitters from either P1, P2, P3, or P4. The resulting simulated movies were SOFI-processed up to the 7^{th} order. The signs of the virtual brightnesses of the virtual emitters (defined in Section 3) were mostly in keeping with the predictions in Fig. 4(ii) for the different orders, with the exception of the out-of-focus P2 emitter regions (details given in Fig. S1) and the 4^{th}-order cumulant for P3. As can be seen from Fig. 4(i), for P1, ρ is distributed in a region with positive lobes for the 2^{nd}- and 6^{th}-order cumulants, a negative lobe for the 4^{th}-order cumulant, and positive/negative transition regions for the 3^{rd}-, 5^{th} -, and 7^{th}-order cumulants. This means that all P1 virtual emitters will exhibit positive virtual brightnesses for the 2^{nd}- and 6^{th}-order cumulants, negative virtual brightnesses for the 4^{th}-order cumulant, and both negative and positive virtual brightnesses for the 3^{rd}-, 5^{th}-, and 7^{th}-order cumulants. We could therefore predict that P1 will exhibit cusp artifacts for the 3^{rd}-, 5^{th}-, and 7^{th}-order cumulants. Small discrepancies could be introduced by the Poisson noise in the signal as shown by the small green region for 6^{th}-order cumulant image for P1. Similarly, for P2, ρ is distributed in a region with positive lobes for the 2^{nd}- and 5^{th}-order cumulants, negative lobes for the 3^{rd}-, 4^{th}-, and 7^{th}-order cumulants, and positive/negative transition regions for the 6^{th}-order cumulant. We could therefore predict cusp artifacts for the 6^{th}-order cumulant of P2. Discrepancies can be introduced by Poisson noise in the signal (5^{th}-order image for P2). Additionally, for P3 and P4, ρ is broadly distributed in the positive/negative transition regions for all the cumulants with orders higher than two and is purely positive only for the 2^{nd}-order cumulant. Therefore, we predicted that P3 and P4 would exhibit cusp artifacts for all the cumulants with orders higher than 2. However, as shown in Fig. 4(iii), the 4^{th}-order of P3 exhibits pure negative virtual brightnesses, which can be attributed to the fact that the percentile population with positive virtual brightnesses in P3 for the 4^{th}-order cumulant was too small. Therefore, the signal is canceled by the large population with negative brightnesses. This is confirmed in the case of P4, where we increased the population with positive virtual brightnesses for the 4^{th}-order cumulant, and the cusp artifacts is as predicted.

In the simulations, we assumed that the blinking phenomenon is governed by a Poisson process, that is, the photons emitted from an individual emitter exhibit a telegraph-noise-like time trajectory. If enough data (i.e., a large number of movie frames) with sufficient statistical significance were to be collected, it would be seen that the estimator for the “on”-time ratio, ρ, would converge to ${\tau _{on}}/({\tau _{on}} + {\tau _{off}})$ (Eq. (4.1)). As the total number of frames in the simulation is reduced, the actual (calculated) value of ρ starts to deviate from its estimated value. This can lead to unexpected cusp artifacts for cumulants of certain orders in the regions that are predicted to be free of these artifacts. Appendix 2 [28] shows that the higher the SOFI order, the higher the number of frames needed for SOFI processing to realize the theoretically predicted cusp-artifact-free images. Moreover, since the high-order virtual brightnesses exhibit greater oscillations, they are more susceptible to heterogeneities in the photophysical properties of the emitters and hence more vulnerable to cusp artifacts.

The next set of simulations (“Simulation-2”) were performed to examine the effects of photobleaching and noise on cusp artifacts (Fig. 5). The time trajectories simulated during Simulation-1 for P1 were stochastically truncated (using Poisson bleaching statistics) to simulate the bleaching events (see Appendix 3 [28] and Fig. 5(i) – (iii)). The predicted signs of the virtual brightnesses of the emitters in P1 are shown on top in the first row of Fig. 5(iv). As can be seen from Fig. 5(iv), when there is no bleaching correction (see below), the signs of the virtual brightnesses displayed in the SOFI images deviate from those predicted. This is because bleaching causes ρ to deviate from its estimated value (Eq. (4.1)), as the “bleached” state is equivalent to an ultralong “off” state, causing a decrease in the ρ values of the emitters and vary from the theoretical predicted ρ value shown in Eq. (4.1). The “bleached” state also affects the assumption of independence of the blinking trajectories of the different emitters, rendering the predictions less reliable. Next, we used a bleaching correction algorithm [14] with minor modifications (Appendix 3 [28]) on the simulated data by dividing the movie into individual blocks of frames, wherein each block had a signal decrease (from the beginning to the end of the block) of fbc = 1% of the overall signal decrease (from the beginning to the end of the whole movie; fbc is called the “bleaching correction factor”). The final cumulant reconstruction was computed as the average of the cumulants of all the time blocks of the same cumulant order. It can be seen that the bleaching-corrected reconstructions (second row in Fig. 5(iv)) exhibit the predicted virtual brightness signs. Next, we examined the effects of the signal-to-noise ratio (SNR) on the bleaching correction. The noise of the EMCCD camera was recorded and added to the simulated bleaching data with altered signal levels (as described in Appendix 3.2 [28]) to ensure SNR levels of 1.47 dB (third row) and −1.33 dB (bottom row). The results of the simulations with bleaching and background noise were SOFI-processed up to the 7^{th} order. As can be seen in Fig. 5, the noise severely degrades the quality of the images in most of the cases, especially for the cumulant orders with a large proportion of positive virtual emitters. However, if the emitters’ blinking parameters were to be tuned to have either purely or primarily negative virtual brightnesses (as for the 4^{th}-order cumulant in the third row, which has pure negative virtual brightnesses, and the 7^{th}-order cumulant in the third row, which primarily has negative virtual brightnesses), a significant enhancement in contrast would be achieved (as shown for this set of simulations). If the majority of the virtual emitters are negative, the negative signal would still create a negative contrast against the positive noise background (as can be seen for the 7^{th}-order cumulant in the third and fourth rows). We would like to note here that the cumulants of the background noise were positive for all orders higher than two; this was because the higher-order (>2) cumulants of a random variable that follows a Poisson distribution remain constant (Appendix 4 [28]).

Here, we would like to highlight the trade-offs between the bleaching correction factor (f_{bc}), noise level, and statistical significance of the blinking trajectories (details available in Appendix 3 [28]). When we decreased f_{bc}, the block size decreased, and the total number of bleaching events occurring within each block decreased as well. Thus, the bleaching effect was suppressed to a greater degree. However, because the smaller f_{bc} led to a decrease in the degree of distinctiveness between the blinking events and noise within each block, the construction of the cumulants became more vulnerable to noise. In addition, a smaller f_{bc} reduces the statistical significance of the blinking trajectories within each block and increases the total number of blocks for constructing the final cumulant. Thus, these two factors counterbalance each other in terms of the effects of f_{bc} on the overall statistical significance of the blinking trajectories.

In the simulations discussed above (Simulation-1 and Simulation-2), the spatial distribution of the ρ value of the different emitters was random. In Fig. 6, we show the results of another set of simulations (“Simulation-3”) in which the ρ values were varied slowly over the range of 0.01 to 0.99 in increments of 0.01 across the spatial field of view, and the SOFI cumulants of the 2^{nd} to 7^{th} orders were calculated. As can be seen from the figure, the number of zero-crossing nodes increased with the order of the cumulants, with the positive and negative domains being in good agreement with the ground-truth predictions. Further, the spatial variations in the blinking statistics can be seen visually as green/red segments (with the cusps concentrated at the segments boundaries).

To summarize this part: mixed populations of positive and negative virtual brightnesses cause cusp artifacts. Knowing the distribution of the virtual brightnesses based on the blinking statistics can allow one to predict cusp artifacts. SOFI cumulants that intrinsically have cusp artifacts are more vulnerable to background noise because of the lowering of the amplitude of the virtual signal. Such signal attenuation is owing to the neutralization of the positive and negative virtual brightnesses as well as the attenuation caused by the small amplitudes of the cumulant wn. Bleaching affects the ability to predict cusp artifacts; however, performing bleaching correction can help overcome this problem, making cusp artifact predictions valid again. Too small of an fbc value can make the SOFI cumulants more vulnerable to background noise (especially when the cumulant order is high) but does not have a significant effect on the statistical significance of the blinking trajectories. Lastly, and more interestingly, in the ideal cases, where the spatial variation of the blinking statistics is very small, cusps can serve as boundaries to segments with similar blinking statistics within the segment, but different blinking statistics between different segments.

## 6. Effects of cusp artifacts on image deconvolution and dynamic-range expansion correction and moment alternative

We have shown in the previous sections that spatial variations in the “on”-time ratio, ρ, lead to positive/negative-value variations in the cumulants, which, in turn, produce cusp artifacts once the absolute value operator has been used with the calculated cumulant images. The subsequent steps in the SOFI protocol for obtaining the final SOFI image include deconvolution (and/or Fourier reweighting) [18] and dynamic-range expansion correction (bSOFI [19], or ldrc as introduced in the accompanying manuscript [21]). In this section, we examine the effects of cusp artifacts on these two final (“post-processing”) steps. We argue that if cusp artifacts are not taken care of properly, their adverse effects could be amplified by the post-processing steps.

In order to determine the factors that contribute to the imperfections in deconvolution and Fourier reweighting, simulations were performed for an ideal collection of positive and negative virtual emitters (Fig. 7(i)) corresponding to a 3^{rd}-order cumulant with zero background and noise; the results are shown in Fig. 7. The blurred image (Fig. 7(ii)) was generated by convolving the ground-truth image with the PSF. In the amplitude display (absolute value) of this image (Fig. 7(iii)), cusps can be seen clearly between the positive/negative regions. The images were generated on fine grids without noise, eliminating the imperfections which could arise from noise or pixilation, allowing us to analyze how the deconvolution step could be influenced by the cusp artifacts and the nature of mixed positive/negative virtual brightnesses. Notice that the imperfections caused by the discrete Fourier transform are still present in the image. As can be seen in the ideal deconvolution panel (Fig. 7(iv)), where under this ideal case, a perfect deconvolution is performed by direct division by the optical transfer function (OTF) in the Fourier space and the subsequent use of the inverse Fourier transform. In this case, we can see that even with a perfect imaging condition, we still expect minor imperfections (shown as the ringing artifacts in Fig. 7(iv)) caused by the discrete Fourier transform with finite boundaries. However, the virtual emitters are recovered in Fig. 7(iv), with their virtual brightnesses remaining unchanged (including the negative ones), suggesting that the influence of discrete Fourier transform is negligible. In the case of Fourier reweighting (mimicking the case of 3^{rd}-order cumulants), the transformed image was divided by a modified OTF with a damping factor (we used the machine epsilon variable in MATLAB, which represents a small number). The Fourier spectrum was subsequently multiplied by an OTF with a wider support. This was equivalent to convolving the ground truth with a smaller Gaussian PSF (Fig. 7(v)). Since the signs of the virtual brightnesses were conserved, the cusp artifacts remained in the high-density regions. On the other hand, when we performed deconvolution with solvers that imposed positivity constraints (for example, MATLAB’s “*deconvlucy*” or “*deconvblind*” functions), deconvolution failed (Fig. 7(vi)). Such failure is due to the presence of negative virtual brightnesses and cusps, which conflict with the positivity constraints used in the deconvolution solver. When the ground-truth is a mixture of positive/negative signal sources, the positivity constraint forces the pixel values to be positive either by taking only the absolute values or by rejecting all the negative ones during the recovery iterations, causing the image no longer to be the result of a convolution between a PSF and a ground-truth image, therefore deconvolution becomes fundamentally invalid.

In cases where a subsequent dynamic range compression step (i.e., the “balancing” step in bSOFI) is performed after deconvolution, the final SOFI image will be even more distorted. Note that the current open source bSOFI package (MATLAB version [27]) uses “*deconvlucy*,” which imposes positivity constraints and therefore improperly handles negative virtual brightnesses. The reliance of bSOFI on positivity constraint and the assumption of perfect deconvolution should be removed. In contrast, Fourier reweighting does not alter the sign of the virtual brightnesses but is very sensitive to noise and does not eliminate cusps.

As shown in Eq. (2.6) and discussed in Section 2, cumulants are additive. Therefore, the *n*^{th}-order cumulant can be expressed as a sum of the cumulants of the individual emitters. This property is, in fact, the reason why cumulant reconstruction was originally chosen [9]. Reconstruction by moments (rather than by cumulants) was not considered owing to the presence of mixed terms that contain signal contributions from multiple individual emitters; such a reconstruction would be mathematically non-rigorous and the resulting physical meaning uninterpretable. Nonetheless, since even-order moments are intrinsically free of cusp artifacts, we pursued moments-based reconstruction as a ** practical** approach and have shown that moments do result in an enhancement in the resolution as compared to the diffraction limit (as shown in the accompanying manuscript [21] and Section 8).

## 7. Comparison of performances

As shown in Fig. 8, we next compared the performances of SOFI, bSOFI, and 6^{th}-order moments (M6) + ldrc for 6^{th}-order reconstructions in terms of cusp artifacts.

All reconstructions were performed on the same data set (movie) of an α-tubulin network in a fixed HeLa cell that was labeled with blinking QD800 QDs (ThermoFisher Scientific, Carlsbad, CA). A total of 2000 movie frames were collected at 33 Hz using an EMCCD camera (Andor USA, South Windsor, CT). The results are shown in Fig. 8. Note that neither deconvolution nor Fourier reweighting was performed for M6 in order to isolate the factors responsible for the enhancement in the resolution and to assess the degree to which cusp artifacts degrade the overall performance. Average, XC2, M6 + ldrc, and bSOFI all show faithful reconstructions in the regions where the filament density is low and the α-tubulins are well separated. On the other hand, XC2 and bSOFI perform better than M6 + ldrc in terms of feature visibility. However, in the regions where the filament density is high, as shown in the boxed region in (i) as well as in (ii), (iv), (vi), and (viii), M6-ldrc out-performs Average, XC2, and bSOFI. This is further confirmed by a comparison of the cross-sectional intensity profiles in (ix). In some cases, for M6 + ldrc, two distinct peaks are observed, in contrast to the case for XC2 (see left-most panel in (ix)), suggesting that the degree of resolution enhancement was higher than that for XC2. However, we would still like to set the lower limit of resolution enhancement for M6 + ldrc to be 2^{1/2}, meaning that the resolution enhancement in this case should at least be similar to that for XC2. An additional theoretical *n*^{1/2}-fold enhancement in the resolution for the *n*^{th}-order moment can be realized when deconvolution is also used. A more detailed performance analysis of the moment-reconstruction method is presented in the accompanying manuscript [21]. Note that, for the bSOFI reconstruction, a deconvolution step was included because balanced cumulant reconstruction is a post-processing step performed after deconvolution.

## 8. Discussion and conclusions

Cusp artifacts had not been identified previously for several reasons: (i) 2^{nd}-order cumulants are always positive and therefore are not susceptible to these artifacts; hence, many studies on SOFI stopped at the 2^{nd} order; (ii) post-processing steps such as balancing and deconvolution mask the existence of cusp artifacts and their origin; and (iii) SR methods in general, and SOFI and SOFI-derivatives reconstructions in particular, are an attempt to solve an ill-defined inverse problem (the location of the emitters in the sample). It is very hard, if not impossible, to identify cusp artifacts from experimental data. It was only after we had performed theoretical analysis of high order cumulants, realistic simulations (including those that considered dark noise, read-out noise, background noise, out-of-focus lighting, spatial variability in the photophysical properties with either limited or sufficient statistics, and the variation of the ground truth feature) and compared the results of the reconstruction algorithms to the ground truths were we able to successfully identify and characterize cusp artifacts. As shown in Fig. 3, once cusps had been identified, a careful reexamination of the experimental data did indeed confirm their existence.

Cusp artifacts are, in fact, hard to avoid. Even if the photophysical properties (blinking and photobleaching) of the emitters are more or less uniform across the sample, the finite acquisition time of a SOFI experiment (usually ∼2,000 to 20,000 frames) is often not long enough to ensure that the statistical significance of the blinking behavior will be reached. Hence, the ρ values exhibit a broad distribution. This, in turn, leads to positive and negative higher-order (>2) cumulant values (Fig. 1). It is only when all the emitters in the sample exhibit a narrow ρ distribution during the data acquisition stage that cusp artifacts can, in principle, be avoided (as shown in Fig. 4 for P1 with cumulants of 2^{nd}, 4^{th}, and 6^{th} orders.). However, a narrow ρ distribution could still be positioned close to a zero crossing of one (or more) of the high-order cumulants, leading to the coexistence of both positive and negative cumulant values.

To minimize the adverse effects of cusp artifacts, the following guidelines should be considered: (1) emitters with uniform photophysical parameters (blinking and bleaching) should be used to ensure an intrinsically narrow ρ distribution; (2) long movies (with a large number of frames) should be acquired (to narrow down the experimental ρ distribution); (3) to the extent possible, the ρ distribution should be tuned to a zero-crossing-free zone (see Fig. 1); and (4) bleaching should be minimized and bleaching correction [14] should be performed on the data set, while following guidelines (1)-(3) for each bleaching correction block.

In summary, we successfully identified cusp artifacts in SOFI reconstructions with orders greater than two. These artifacts have either missed completely or interpreted erroneously in previous studies. In this work, we proposed a virtual-emitter interpretation of high-order SOFI cumulants, using which we could determine the theoretical origins of cusp artifacts. We performed a series of realistic simulations that provided insights into the origin of cusp artifacts and have proposed guidelines on how to minimize their adverse effects. We were able to demonstrate that moment-based reconstructions can improve the resolution and sometimes even eliminate these artifacts. We also suggest guidelines on how to screen for improved probes that could minimize cusp artifacts in high-order SOFI cumulants. In addition, based on the novel theoretical framework of virtual-emitter interpretation, our work provides insights into high-order cumulants as well as cusp artifacts that could potentially inspire positive utilization of the abnormal virtual brightness distribution and cusps where ρ can serve as an indicator.

## Funding

National Science Foundation (DMR-1548924); Dean Willard Chair Fund; Lawrence Livermore National Laboratory.

## Acknowledgements

The first author thanks Dr. Jianmin Xu for advises and training on cell culturing and immunostaining. We also thank Ms. Yingyi Lin, Mr. Xi Lin, and Mr. Sungho Son for their help with the project as undergraduate student researchers. Finally, we acknowledge the computational and storage services provided by the Hoffman2 Shared Cluster maintained by the Research Technology Group of the UCLA Institute for Digital Research and Education. This work was partially performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344.

## Author Contributions

S. W. and X. Y. designed the study; X. Y. performed the experiments, performed the simulations and analyses and wrote the associated codes; and S. W. and X. Y. wrote the manuscript. The authors declare no conflict of interests.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **S. W. Hell and J. Wichmann, “Breaking the diffraction resolution limit by stimulated emission: stimulated-emission-depletion fluorescence microscopy,” Opt. Lett. **19**(11), 780–782 (1994). [CrossRef]

**2. **E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacino, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science **313**(5793), 1642–1645 (2006). [CrossRef]

**3. **M. G. Gustafsson, “Surpassing the lateral resolution limit by a factor of two using structured illumination microscopy,” J. Microsc. **198**(2), 82–87 (2000). [CrossRef]

**4. **M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nat. Methods **3**(10), 793–796 (2006). [CrossRef]

**5. **A. M. Sydor, K. J. Czymmek, E. M. Puchner, and V. Mennella, “Super-resolution microscopy: from single molecules to supramolecular assemblies,” Trends Cell Biol. **25**(12), 730–748 (2015). [CrossRef]

**6. **E. Betzig, “Single Molecules, Cells, and Super-Resolution Optics (Nobel Lecture),” Angew. Chem., Int. Ed. **54**(28), 8034–8053 (2015). [CrossRef]

**7. **W. E. Moerner, “Single-Molecule Spectroscopy, Imaging, and Photocontrol: Foundations for Super-Resolution Microscopy (Nobel Lecture),” Angew. Chem., Int. Ed. **54**(28), 8067–8093 (2015). [CrossRef]

**8. **S. W. Hell, “Nanoscopy with focused light,” Ann. Phys. **527**(7-8), 423–445 (2015). [CrossRef]

**9. **T. Dertinger, R. Colyer, G. Iyer, S. Weiss, and J. Enderlein, “Fast, background-free, 3D super-resolution optical fluctuation imaging (SOFI),” Proc. Natl. Acad. Sci. **106**(52), 22287–22292 (2009). [CrossRef]

**10. **T. Dertinger, M. Heilemann, R. Vogel, M. Sauer, and S. Weiss, “Superresolution optical fluctuation imaging with organic dyes,” Angew. Chem. **122**(49), 9631–9633 (2010). [CrossRef]

**11. **P. Dedecker, G. C. Mo, T. Dertinger, and J. Zhang, “Widely accessible method for superresolution fluorescence imaging of living systems,” Proc. Natl. Acad. Sci. **109**(27), 10909–10914 (2012). [CrossRef]

**12. **X. Zhang, X. Chen, Z. Zeng, M. Zhang, Y. Sun, P. Xi, J. Peng, and P. Xu, “Development of a reversibly switchable fluorescent protein for super-resolution optical fluctuation imaging (SOFI),” ACS Nano **9**(3), 2659–2667 (2015). [CrossRef]

**13. **A. M. Chizhik, S. Stein, M. O. Dekaliuk, C. Battle, W. Li, A. Huss, M. Platen, I. A. Schaap, I. Gregor, and A. P. Demchenko, “Super-resolution optical fluctuation bio-imaging with dual-color carbon nanodots,” Nano Lett. **16**(1), 237–242 (2016). [CrossRef]

**14. **T. Dertinger, A. Pallaoro, G. Braun, S. Ly, T. A. Laurence, and S. Weiss, “Advances in superresolution optical fluctuation imaging (SOFI),” Q. Rev. Biophys. **46**(2), 210–221 (2013). [CrossRef]

**15. **S. Cho, J. Jang, C. Song, H. Lee, P. Ganesan, T.-Y. Yoon, M. W. Kim, M. C. Choi, H. Ihee, and W. Do Heo, “Simple super-resolution live-cell imaging based on diffusion-assisted Förster resonance energy transfer,” Sci. Rep. **3**(1), 1208 (2013). [CrossRef]

**16. **F. Hertel, G. C. Mo, S. Duwé, P. Dedecker, and J. Zhang, “RefSOFI for mapping nanoscale organization of protein-protein interactions in living cells,” Cell Rep. **14**(2), 390–400 (2016). [CrossRef]

**17. **L. Kisley, R. Brunetti, L. J. Tauzin, B. Shuang, X. Yi, A. W. Kirkeminde, D. A. Higgins, S. Weiss, and C. F. Landes, “Characterization of porous materials by fluorescence correlation spectroscopy super-resolution optical fluctuation imaging,” ACS Nano **9**(9), 9158–9166 (2015). [CrossRef]

**18. **T. Dertinger, R. Colyer, R. Vogel, J. Enderlein, and S. Weiss, “Achieving increased resolution and more pixels with Superresolution Optical Fluctuation Imaging (SOFI),” Opt. Express **18**(18), 18875–18885 (2010). [CrossRef]

**19. **S. Geissbuehler, N. L. Bocchio, C. Dellagiacoma, C. Berclaz, M. Leutenegger, and T. Lasser, “Mapping molecular statistics with balanced super-resolution optical fluctuation imaging (bSOFI),” Opt. Nanosc. **1**(1), 4 (2012). [CrossRef]

**20. **S. Jiang, Y. Zhang, H. Yang, Y. Xiao, X. Miao, R. Li, Y. Xu, and X. Zhang, “Enhanced SOFI algorithm achieved with modified optical fluctuating signal extraction,” Opt. Express **24**(3), 3037–3045 (2016). [CrossRef]

**21. **X. Yi, S. Son, R. Ando, A. Miyawaki, and S. Weiss, “Moments reconstruction and local dynamic range compression of high order Superresolution Optical Fluctuation Imaging,” Biomed. Opt. Express **10**(5), 2430–2445 (2019). [CrossRef]

**22. **M. G. Kendall, Alan Stuart, and J. K. Ord, * The Advanced Theory of Statistics* (Wiley, 1968), Vol. 3.

**23. **K. T. Shimizu, R. G. Neuhauser, C. A. Leatherdale, S. A. Empedocles, W. Woo, and M. G. Bawendi, “Blinking statistics in single semiconductor nanocrystal quantum dots,” Phys. Rev. B **63**(20), 205316 (2001). [CrossRef]

**24. **P. Charles, *Digital Video and HDTV Algorithms and Interfaces *(Morgan Kaufmann Publishers, 2003), 260, 630.

**25. **S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A **9**(1), 154–166 (1992). [CrossRef]

**26. **X. Yi, “SR_Simu3D” (2018), retrieved https://xiyuyi.github.io/SR_simu3D/.

**27. **e. a. Marcel Leutenegger, Balanced super-resolution optical fluctuation imaging, 2012.

**28. **X. Yi and S. Weiss. “Cusp-artifacts in high order superresolution optical fluctuation imaging.” bioRxiv: 545574 (2019).