Evaluation of sted super-resolution image quality by image correlation spectroscopy (QuICS)

Quantifying the imaging performances in an unbiased way is of outmost importance in super-resolution microscopy. Here, we describe an algorithm based on image correlation spectroscopy (ICS) that can be used to assess the quality of super-resolution images. The algorithm is based on the calculation of an autocorrelation function and provides three different parameters: the width of the autocorrelation function, related to the spatial resolution; the brightness, related to the image contrast; the relative noise variance, related to the signal-to-noise ratio of the image. We use this algorithm to evaluate the quality of stimulated emission depletion (STED) images of DNA replication foci in U937 cells acquired under different imaging conditions. Increasing the STED depletion power improves the resolution but may reduce the image contrast. Increasing the number of line averages improves the signal-to-noise ratio but facilitates the onset of photobleaching and subsequent reduction of the image contrast. Finally, we evaluate the performances of two different separation of photons by lifetime tuning (SPLIT) approaches: the method of tunable STED depletion power and the commercially available Leica Tau-STED. We find that SPLIT provides an efficient way to improve the resolution and contrast in STED microscopy.

Here, we introduce a simple algorithm to evaluate systematically and in an unbiased way the quality of STED images by image correlation spectroscopy (QuICS). Image correlation spectroscopy (ICS) is a general and versatile method to quantitatively analyze fluorophore distribution in microscopy images 16 . ICS can be used to extract parameters such as size 17 , distances 18,19 and aggregation state 20 from static images and dynamic parameters such as diffusion coefficient 21 and velocity 22 from time-resolved images. In this work, we focus only on the analysis of static images. We apply ICS to extract three quantities that are related to the quality of the super-resolved image: the width of the autocorrelation function, related to the spatial resolution; the brightness, related to the image contrast; the relative noise variance, related to the signal-to-noise ratio of the image. Within this study, we describe how the modulation of image acquisition parameters can influence STED efficiency and the image formation of DNA replication sites in U937-PR9 cells, an in vitro model of leukemia 23 . Our study reveals that to optimize the imaging conditions for a given sample, a balance between different parameters must be found. We found a valid solution to this elusive balance by applying the method of Separation of Photons by LIfetime Tuning (SPLIT) 24,25 to STED microscopy. In particular, we show SPLIT images obtained using the method of tunable STED depletion power 25,26 , or acquired by a recently developed, lifetime-based commercial setup (the Leica Tau-STED microscope). QuICS analysis reveals that SPLIT images have higher resolution and non-reduced brightness and noise parameters, compared to their counterpart STED images.
We developed the QuICS algorithm based on the growing need for analysis of nuclear processes performed at the level of individual cells, also taking into account that certain events typically occur in a relatively small fraction of cells in the population at any given time (i.e., events taking place in a specific phase of the cell cycle). Recent advances observed a considerable variability and heterogeneity in genome organization at the single-cell level 27 . Imaging and super-resolution can thus provide a unique view of nuclear organization and functions in intact cell nuclei.

Results
Autocorrelation function as a source of information about image quality. In order to extract a series of parameters associated to an image I(x,y), we start by calculating a radial autocorrelation function (ACF) G(ρ), where ρ represents the spatial lag in the radial direction. This function is calculated by performing an angular average on the two-dimensional ACF G(δ x ,δ y ), where δ x and δ y are the spatial lag variables (see "Methods"). In general, the function G(ρ) contains information on all the intensity variations in the image, including fluctuations due to statistical noise. Let's call G NF (ρ) the noise-free correlation function, i.e. the corresponding function in the absence of noise. By fitting G NF (ρ) to a Gaussian model (see Eq. (5)), we extract the amplitude G NF (0) and the width parameter w. We define the following three quantities ( Fig. 1):  www.nature.com/scientificreports/ where we have indicated I av as the average intensity value over all the pixels of the image. In order to understand the physical meaning of R, B and N, let's assume, for simplicity, that the sample contains randomly distributed point-like fluorescent particles so that the corresponding image is the convolution of the emitters and the Point Spread Function (PSF) of the optical system. The PSF in the lateral direction is assumed to be a Gaussian function of the form exp(− 2r 2 /w PSF 2 ), where w PSF represents the PSF waist (1/e 2 width) and is related to the Full Width at Half Maximum (FWHM) of the same PSF by the relationship FWHM PSF = (2ln2) ½ w PSF . The resulting autocorrelation function will be a Gaussian of the form exp(− ρ 2 /w PSF 2 ) and can be fitted by Eq. (5). In this case, R corresponds to the resolution of the optical system expressed in terms of the FWHM of the PSF, R = FWHM PSF . More in general, since the sample may contain features of finite size, it will be R ≥ FWHM PSF . Thus, the estimated resolution of the optical system is at least equal to R. Note that, even if the PSF of the optical system is not a Gaussian, Eq. (5) is a convenient approximation to extract the parameter R from the ACF.
The quantity B in Eq. (2) is called brightness 28,29 and is equal to σ p 2 /I av where σ p 2 is the variance of the intensity due to the particles. The brightness of the particles depends on the number of fluorophores per particle and on the actual brightness of the fluorophores at the specific imaging settings (e.g. excitation intensity level, detector gain, pixel dwell time). Let's assume that, in addition to the signal from the particles, whose average intensity is I av,p , there is a uniform background signal, with average intensity I av,bkgd , so that I av = I av,p + I av,bkgd . The brightness is given by B = σ p 2 /(I av,p + I av,bkgd ). Thus, a reduction of the brightness parameter B is related to a decrease in the contrast of the particles in the image.
Finally, we note that G NF (0) = σ p 2 /I av 2 and G(0) = (σ p 2 + σ noise 2 )/I av 2 , where σ noise 2 is the variance of the intensity due to noise. Thus, the quantity N in Eq. (3), N = σ noise 2 /σ p 2 , represents the variance of the noise normalized to the variance of the particles (relative noise variance). We use N to quantify the noise level in the image, where the limits N = 0 indicates no noise, and N = 1 indicates that intensity fluctuations due to noise are comparable to those due to the particles.
The noise-free correlation function, G NF (ρ), required for this analysis, can be obtained by cross-correlating two statistically independent realizations of the same image, in analogy to what is done in Fourier Ring Correlation (FRC) methods [12][13][14] . The acquisition of two statistically independent images is straightforward in singlemolecule localization microscopy but can be a more difficult task with other super-resolution techniques 15 . In order to estimate G NF (ρ) from a single acquired image, we propose performing a fit of the autocorrelation function G(ρ) either skipping the zero lag point (which contains the contribution of intensity fluctuations due to statistical noise), or cross-correlating two statistically independent images obtained by down-sampling the original image 30 . In our experimental data, the two approaches provided similar results (Supplementary Figure 1).

Tuning of STED depletion power.
During STED image acquisition, the most intuitive way to improve the resolution is to increase the STED beam's intensity. To validate our method as a function of the STED beam intensity, we acquired images of U937-PR9 cells samples in which we were able to visualize the DNA replication thanks to the incorporation of the nucleoside Ethynyl deoxyuridine (EdU) to the newly replicated DNA strand. We then coupled the EdU molecule to an azide molecule carrying the Alexa fluorophore, taking advantage of a Cu-catalyzed Click-iT reaction. During microscopy acquisition, we were careful to select actively replicating cells with DNA replication sites spread all over the nucleus and thus more suitable for resolution evaluation analysis.
First, we acquired a confocal and a STED image of a cell nucleus (Fig. 2a, upper row) by applying a relatively small depletion beam power (9 mW) and the two line profiles of the same structure were plotted to compare the achieved resolution ( Fig. 2b, upper panel). Then, we acquired a confocal and a STED image of a second nucleus, doubling the power of the depletion beam ( Fig. 2a, lower row), and we compared the line profiles of the same structure from the confocal and the STED image. The line profile analysis yielded FWHM = 189 nm and FWHM = 212 nm for the two peaks detected at 9 mW, and FWHM = 160 nm and FWHM = 178 nm for the two peaks detected at 18 mW. However, this result depends on the specific structures selected for the line profile analysis and does not take into account the totality of labeled sites in the whole cell nucleus. Therefore, we acquired at least ten STED images with each of the two different depletion beam powers, and we calculated the autocorrelation function in order to obtain the average resolution R related to the entire nuclei. As a result, we obtained that the doubling of the STED depletion beam power from 9 to 18mW led to an improvement of spatial resolution from R = 234 ± 3 nm to R = 213 ± 3 nm (mean ± s.d., Fig. 2c). This result is in keeping with the line profile analysis, as expected, and represents an average of the whole nucleus structures. The obtained values of R strongly depend on the average apparent size of the structures (i.e. replication foci) in the images, meaning the molecular volume plus the resolution, and therefore their values are larger than the maximum resolution of the optical system.
In contrast, we observed that the image brightness B was significantly lower for the images acquired with the higher STED depletion power (Fig. 2d). We interpreted this reduction of B as a reduction in the image contrast. In fact, we compared images acquired exactly with all the same instrumental settings (e.g. same excitation power, same detector gain, same pixel dwell-time) other than the STED depletion power. In this case, the action of STED reduces the average intensity per pixel due to the particles but does not decrease the average intensity of the background signal (originating, for instance, from undepleted out-of-focus fluorescence signal 24,31,32 ). As explained in the previous section, this causes a reduction of the parameter B.
Finally, we observed that the relative noise variance N was higher at the higher STED depletion power (Fig. 2e). This is in line with the expected reduction of signal-to-noise ratio at increasing STED depletion power. To evaluate how increasing averages would influence the quality of the image, we acquired sequential STED images of the same cell, with a depletion beam power of 18 mW (see "Methods" for a detailed description of the sequential acquisition settings). From each sequential acquisition, we generated STED images with a different number of line-averages (Fig. 3a) and applied the QuICS algorithm. As expected, the image's noise significantly decreased following doubling of the number of lines averaged (Fig. 3b). On the other hand, the average resolution R did not improve with an increasing number of lines in the average (Fig. 3c). The brightness B decreased as a function of the number of averages (Fig. 3d) indicating a reduction of the contrast.  www.nature.com/scientificreports/ To interpret these results, we monitored photobleaching as a function of the number of averages of the STED image (Fig. 3e). Photobleaching was calculated as the percentage reduction of average fluorescence intensity with respect to the initial value. We observed that each line averages-acquisition step induced a significant increase in photobleaching of the sample's fluorophores (Fig. 3e). Consequently, the image contrast decreased, thus leading to a brightness reduction as a function of the number of averages (Fig. 3d,f). These data also show that, in our samples, for photobleaching levels above about 40%, there is no significant improvement in the signal-to-noise ratio of the images (Fig. 3b,f). Thus, in case it cannot be avoided, photobleaching should be at least kept below this level.

Comparison between SPLIT and STED imaging.
To increase the spatial resolution of a STED microscope, the most straightforward way is to increase the depletion beam's intensity. However, as we have seen, this may reduce the contrast and signal-to-noise of the images, quantified in QuICS via the brightness and noise parameters. Here, we examine the advantages of increasing spatial resolution via application of Separation of Photons by LIfetime Tuning (SPLIT) 24 . The SPLIT method provides an increase in spatial resolution by decoding the spatial information encoded into an additional channel. The first reported SPLIT configuration exploited, as an additional channel, the fluorescence lifetime gradient induced by a continuous-wave STED beam 24,33 . Subsequent studies demonstrated that SPLIT is not limited to analysis of fluorescence lifetimes. SPLIT could also be applied to stacks of STED images obtained with tunable depletion power, with the depletion power used as the additional channel for SPLIT 25,26 , or even to structured illumination microscopy images 34 .
As described in Fig. 4a, we first applied the SPLIT method to stacks consisting of two STED images at different depletion power: a confocal (0 mW STED depletion power) and a STED image (18 mW STED depletion power). In this case, the fluorescence intensity variations due to the tuning of the STED depletion power, allow the separation of the contributions from fluorophores in the center or the periphery of the PSF (Fig. 4a). Since the excitation intensity can also be easily tuned along the stack 26 , we set the excitation level of the confocal image so that it induced negligible photobleaching. In this way, the data acquisition for SPLIT was straightforward and did not induce more photobleaching than the acquisition of the STED image alone. Figure 4b shows application of this approach to imaging of replication foci in a U937-PR9 cell in S phase. Shown are the confocal, the STED image and the resulting SPLIT image. We compared the line profiles of the same structure and we observed a resolution improvement from FWHM = 147 nm and FWHM = 154 nm, for the two peaks detected in the STED line profile, to FWHM = 135 nm and FWHM = 107 nm, for the two peaks detected in the SPLIT line profile (Fig. 4c). QuICS analysis of at least ten samples, revealed a significant improvement in the average resolution of the SPLIT image (R = 129 ± 9 nm, mean ± s.d.) compared to the STED image (R = 213 ± 9 nm, mean ± s.d.) (Fig. 4d). Notably, this improvement in resolution is not achieved at the expense of the image brightness (Fig. 4e) or the signal-to-noise ratio (Fig. 4f). Figure 4g shows an image of replication foci in a U937-PR9 cell in late S-phase acquired with the Leica Tau-STED microscope. Here the SPLIT image (i.e. the Tau-STED image) is compared with the STED image and with a time-gated STED image (time gate = 1-8 ns). QuICS analysis indicates an improvement of resolution from R = 231 nm (STED image) and R = 225 nm (gated-STED image) to R = 181 nm (Tau-STED image) (Fig. 4h). The brightness is significantly higher in the Tau-STED image than in the STED and gated-STED images (Fig. 4i). This increase in brightness is probably due to the improvement of contrast provided by SPLIT, which has the capability of filtering out background signal originating for instance by direct excitation from the STED beam 24 . The SPLIT image has higher SNR than the STED and gated-STED images, in line with the overall reduction of background in the image (Fig. 4j) and in keeping with previous studies 35 . The large variation among the QuICS parameters can be due to the different stage of the cell cycle, and thus the different EdU distribution pattern, of cells acquired for this measure.

Discussion
Applications of super-resolution microscopy to biology are increasing. However, despite the availability of several types of commercial setups, optimization of the conditions of imaging still requires some degree of expertise. It is important to find the conditions that maximize the quality of the image, paying attention to the onset of potentially degrading effects such as fluorophore photobleaching. Our approach provides an unbiased measurement of the super-resolution image quality based on the three parameters R, B, N defined by Eqs. (1), (2), (3). We note that R and N can be readily used to compare the resolution and signal-to-noise ratio of images acquired under different conditions. On the other hand, B depends on the image contrast, but also on many instrumental factors (e.g., the excitation intensity level, the detector gain, and the pixel dwell-time) that should be considered when performing any comparison.
There is an important difference between QuICS and the FRC method. The FRC resolution merges into a single parameter information about both the relevant spatial frequencies and the noise content of an image. In other words, the FRC resolution describes the length scale below which the image lacks signal content 13 . In QuICS, the resolution parameter R contains average information on the characteristic size (e.g. specimen features, PSF of the optical system) whereas the parameter N contains information on the noise content of the image. In the limit of infinitely high signal-to-noise ratio, the two values of resolution extracted by QuICS and FRC are the same. Conversely, for low signal-to-noise ratio, we expect the FRC value to increase whereas the QuICS resolution to remain constant, since it represents the average apparent size of particles in the image (or the size of the PSF, in the limit of point-like particles). Thus, an advantage of QuICS is that the same algorithm can be used not only to evaluate the image quality but also to quantify biophysical parameters such as the size and the molecular brightness, important in many biophysical applications 16,17,28,29,36 .  www.nature.com/scientificreports/ The combination of super-resolution microscopy with the correlation spectroscopy toolbox undoubtedly offers several advantages 37 . Here, we have shown how analysis of an angle-averaged, image correlation function can provide useful hints on the optimization of the imaging conditions. As a case study, we have focused our attention on STED imaging of DNA replication foci in fixed U937 cells. However, even if not demonstrated, we expect that the approach can be adapted to images containing arbitrary features (for instance cytoskeletal structures). Similarly, we expect that it can be used to evaluate the quality of images acquired in confocal microscopy or other types of super-resolution techniques.

Methods
Cell culture and treatments. U937-PR9 cells were cultured in RPMI-1640 medium (Sigma Aldrich R7388) supplemented with 1% penicillin/streptomycin (Sigma-Aldrich P4333) and 10% fetal bovine serum (Sigma-Aldrich F9665) and maintained at 37 °C and 5% CO 2 . U937-PR9 were seeded on poly-l-lysine (Sigma-Aldrich P8920) coated glass coverslips immediately before experiments. Cells were incubated with 10 µM of the synthetic nucleoside 5-Ethynyl-2′-deoxyuridine (EdU) (Thermo Fisher Scientific) for 25 min at 37 °C and 5% CO 2 . Image acquisition. Images of Figs. 2, 3 and 4b were acquired on a Leica TCS SP5 gated-STED microscope, using an HCX PL APO 100× 100/1.40/0.70 oil immersion objective lens (Leica Microsystems, Mannheim, Germany). Emission depletion was accomplished with a 592 nm STED laser and depletion beam power was measured after the objective. Excitation was provided by a white laser at the desired wavelength for each sample. Alexa488 was excited at 488 nm and its fluorescence emission detected at 500-560 nm, with 1.5-5 ns time gating using a hybrid detector (Leica Microsystems). 512 × 512 pixel images were acquired with a pixel size of 20 nm.

EdU
For the experiment reported in Fig. 3, the first four STED images were acquired with 16 scans average per pixel line, while the fifth image was acquired with 64 scans average per line. We designed this experiment intending to mimic an acquisition with 128 scans average per line and to be able to monitor the trend of the resolution, noise, and brightness after each 16-averages step. Besides, before and after each STED image, we also acquire a confocal image, in order to monitor the trend of the photobleaching after each STED acquisition. To do so we carefully choose the confocal acquisition parameters in order to induce a confocal-related negligible photobleaching. The final images were then obtained by combining and averaging frame of STED images after each acquisition step so that the resulting image had 16 + 16 averages, 16 + 16 + 16 averages and so on (Fig. 3a).
Images of Fig. 4g were acquired on a Leica Stellaris 8 Tau-STED microscope, using an HC PL APO CS2 100×/1.40 oil immersion objective lens (Leica Microsystems, Mannheim, Germany). Emission depletion was accomplished with a 775 nm STED laser. Excitation was provided by a white light laser at the desired wavelength for each sample. Alexa594 was excited at 561 nm and its fluorescence emission detected at 570-620 nm using a hybrid detector (Leica Microsystems). 1024 × 1024 pixel images were acquired with a pixel size of 14 nm.
Generation of SPLIT images. Separation of photon by lifetime tuning (SPLIT) images in Fig. 4b were generated using the method of tunable depletion power 25,26 . A simplified version of the algorithm described in 24 was implemented in Matlab and applied to two-frame stacks consisting of a confocal and a STED image.
QuICS algorithm. The QuICS analysis was performed in MATLAB (The MathWorks) using a custom code.
Given an image I(x,y), a two-dimensional (2D) image correlation function G 2D (δ x ,δ y ) was calculated as: where δ x and δ y are the spatial lag variables, I(x,y) is the fluorescence intensity detected at pixel (x,y), the angle brackets indicate averaging over all the selected pixels of the image. The numerator in Eq. (4) was calculated by a 2D fast Fourier transform algorithm. Before calculation, a region of interest (ROI) corresponding to the nucleus was defined using the counterstain signal and the corresponding mask has been applied to the image as described previously 18,38 . This step is useful to minimize the effects of nuclear borders on the correlation functions.
The 2D correlation function was then converted into a one-dimensional radial correlation function, G(ρ), by performing an angular mean 17 . To this aim, the spatial lag variables δ x and δ y were expressed via polar coordinates, defined by the relations δ x = ρ•cosθ and δ y = ρ•sinθ, where ρ is the radial lag and θ is the angle lag. Then, the radial correlation function G(ρ) was calculated as G(ρ) = <G 2D (ρ,θ)> θ , where G 2D (ρ,θ) is the 2D correlation function expressed via polar coordinates and the angle brackets indicate averaging over the angular coordinate.
To estimate the noise-free correlation function from a single image, we performed a Gaussian fit of the radial correlation function G(ρ) by skipping the zero lag point: www.nature.com/scientificreports/ where the width parameter w corresponds to the 1/e 2 of a Gaussian function and it is related to the Full Width Half Maximum (FWHM) by the relationship w = FWHM/(2ln2) 1/2 ; G NF (0) represents the amplitude; G NF (∞) represents an offset value. The fitting range was determined as follows. The value ρ min was set to one pixel, in order to skip the zero-lag point which contains contribution from statistical noise. The value ρ max was set, by visual inspection of the data, in such a way to fit a single Gaussian component (Supplementary Figure 1). As an alternative approach, we generated two independent images I′(x,y) and I″(x,y) by downsampling the image I(x,y) to half the size. The image I′(x,y) was obtained by averaging the intensity of pixel (i,j) with that of pixel (i + 1, j + 1), with i + j even. The image I″(x,y) was obtained by averaging the intensity of pixel (i,j) with that of pixel (i + 1, j + 1), with i + j odd. The images I′(x,y) and I″(x,y) were then resampled back to the original size. The 2D cross-correlation function was calculated as: The 2D cross-correlation function was then converted into a one-dimensional radial cross-correlation function, G cc (ρ) and fitted with Eq. (5) by setting ρ min = 0. The two approaches yielded similar results in our experimental data (Supplementary Figure 1).
A user-friendly version of the Matlab code is available at https:// github. com/ llanz ano/ QuICS.