Spectral cross-cumulants for multicolor super-resolved SOFI imaging

Super-resolution optical fluctuation imaging provides a resolution beyond the diffraction limit by analysing stochastic fluorescence fluctuations with higher-order statistics. Using nth order spatio-temporal cross-cumulants the spatial resolution and the sampling can be increased up to n-fold in all spatial dimensions. In this study, we extend the cumulant analysis into the spectral domain and propose a multicolor super-resolution scheme. The simultaneous acquisition of two spectral channels followed by spectral cross-cumulant analysis and unmixing increases the spectral sampling. The number of discriminable fluorophore species is thus not limited to the number of physical detection channels. Using two color channels, we demonstrate spectral unmixing of three fluorophore species in simulations and experiments in fixed and live cells. Based on an eigenvalue/vector analysis, we propose a scheme for an optimized spectral filter choice. Overall, our methodology provides a route for easy-to-implement multicolor sub-diffraction imaging using standard microscopes while conserving the spatial super-resolution property.


Supplementary Note 1: Theory of spectral unmixing using spectral cross-cumulants n th order spectral cross-cumulant analysis between two adjacent physical spectral channels
For the case of 2 physical spectral channels (here: transmission channel T and reflection channel R, collecting fluorescence light in a specific wavelength range. A according to the combined spectral response of all the filters implemented in the microscope, with corresponding transmission and reflection coefficients per fluorophore species i) and a total of Nc fluorophore species we write: with ( ) being the intensity distribution of species i measured on detector pixel r. This linear system cannot be inverted to solve for the images of the isolated fluorophore species for more unknowns 1 ( ),

General case of n th -order spectral cross-cumulant analysis and Np physical color channels
For the most general case of Np physical spectral channels (psc) and Nc fluorescent species i, we can define a corresponding proportion , of the intensity that is directed into the specific spectral channel. With ( 1 , … , ) the set of physical spectral channels ∈ �1, … , p � denoting the cross-cumulant that is computed using pixels from the physical spectral channel . { ( , )} = ; denotes the n th -order cumulant of the different fluorophore species i. This computation of additional channels is the key to enable unmixing by inversion of the linear system of equations.

Single-color single-species cumulant for m emitters:
For m fluorescent emitters of a single fluorophore species i recorded in a single color channel, the n thorder cumulant can be written as where ( ) is the the spatial distribution of the molecular brightness, , � on, ; � is the n th -order cumulant of a Bernoulli distribution with on-time ratio on, = on, on, + off, and ( ) is the system PSF for the fluorophore species i 2 .

Flattening:
Using cross cumulants, virtual pixels are calculated in between the physical pixels acquired by the camera.
Subsequently, proper weights are assigned to these virtual pixels in the so-called flattening operation assuming a known PSF 3 (see Equation S6) or optimal weights are calculated using a computationally demanding approach based on jackknife resampling 4 .
In this study, we used the simple yet effective approach of weighing to the same mean within sub-grids of the image. Calculating the n th order cross-cumulant, we obtain n-1 virtual pixels in between each pair of pixels of the original pixel grid. For example, in the case of the 2 nd order cumulant, we generate 3 virtual pixels for each physical pixel (i.e. there are 4 "pixel types" in the new, finer grid). This new grid can be divided into 4 mutually shifted sub-grids (each composed out of pixels of the same "pixel type"). These sub-grids represent the same image shifted by ps/n, where ps is the projected pixel size of the original image and n is the cumulant order. Assuming that these mutually shifted subsampled versions of the full image are supposed to have the same mean, the flattening can be performed by simply normalizing the sub-grids to the same mean value as the mean of the original image (i.e. sub-grid composed of the physical pixels).

Supplementary Note 2: Transmission and reflection coefficients
From spectral data of the used fluorophores and filters The known fluorophore emission spectra are weighted with the spectral response curves of the reflection and transmission channels obtained from transmission data of the different (dichroic) filters that are

Experimental determination
Cells are labelled with a single fluorophore species i and widefield images across the same filter combination as in the multicolour experiments are obtained in the reflection and emission channel.
Subsequently, the transmission and reflection coefficients for fluorophore species i can be calculated by summing the background corrected transmitted and reflected intensity, respectively, and normalizing by the total emission.

Supplementary Note 3: Optimization of multicolor SOFI via eigenvalue and eigenvector analysis
The selection of the best combination of dyes and adequate filter sets can be an overwhelming and challenging task. In this section, we discuss a systematic approach to guide potential users in this process.
Each fluorophore is characterized by its emission spectrum ( ). If we assume that we use a perfect dichroic filter with a transmission function where ( ) is the Heaviside step function and is the characteristic wavelength of the dichroic, we can express the transmission and reflection coefficients , as where we assume that the spectrum ( ) is already normalized. Similarly, we can express the unmixing matrix as a function of The matrix ( ) is not invertible if one or more of its eigenvalues are equal to 0. In practice, the linear system is degraded by several noise sources and a matrix with an eigenvalue close to 0 is likely to be unstable. In order words, we need to optimize the product | 1 || 2 || 3 | or maximize the smallest eigenvalue. Using the three eigenvalues of ( ) as a function of as well as their product. As expected, when is smaller than 560nm or greater than 700nm, at least one eigenvalue tends to 0, meaning that at least two dyes are completely reflected or transmitted. We observe a maximum for the eigenvalues product at around 600nm, which corresponds to the theoretical optimal splitting wavelength. In the case of a real dichroic with non-idealized reflection and transmission characteristics, we have to rewrite and as We can then rank any dyes and dichroic combination and select the 1 ( ), 2 ( ), 3 ( ) and ( ) that produces the least singular matrix. In our case, the choice of ZT594RDC with a splitting at 594nm results in the eigenvalues: 1.05, 1.05 and 0.245 (product of 0.27), validating the choice of the dichroic for this specific set of dyes.
The generalization to more channels and dyes is straightforward and will just add additional vectorial components to this eigenvalue/vector analysis.
The unmixing matrix ( ) corresponding to the experiments and simulations in this work is diagonalized and Λ = λ , as the three eigenvectors form a basis. As ( ) needs to be invertible, we choose λ such that λ ≠ 0 (see discussion above). Thus, the dimension of the image of M is the same as the dimension of its domain and the rank of the unmixing matrix equals the number of color channels. This confirms that our cumulant analysis indeed provides an independent third channel allowing the unmixing of the three fluorophore species.

Co-registration based on calibration measurements with beads
An affine transformation and bilinear interpolation (Matlab) based on a calibration measurement with fluorescent beads that can be detected in both physical channels is applied to the transmission channel.
We typically use ∅ 0.2 μm TetraSpeck beads or ∅ 0.17 μm orange beads from the PSF calibration kit (Invitrogen) dried on glass and covered in the supplied immersion medium. For the measurements presented in the paper the co-registration precision (vector sum of coordinate displacement in original vs. co-registered channel) is ~ 10-30 nm. The virtual spectral cross-cumulant channel has by construction half the registration error with respect to the physical color channels. The unmixing step is a linear operation, thus the registration between the final unmixed color channels is comparable to the coregistration precision of the physical color channels -an order of magnitude better than the attainable resolution for second order analysis. It is noteworthy that careful co-registration is important for crosscumulant analysis, but sufficient accuracy is routinely achieved.

Co-registration based on image correlation
The temporal standard deviation of the transmission and reflection channels is computed to generate two background free images. An affine transformation and bilinear interpolation (Matlab) based on the normalized cross-correlation between the two images is then applied to all the frames of the transmission channel.

Supplementary Note 5: Simulations -multicolor SOFI with spectral unmixing
To this end, we investigated the influence of different photophysical properties of the fluorophores on the performance of our multicolor analysis in simulations (for details on the simulations, please consult the main text and the Methods section). We simulate thin, densely labelled filaments that are partially overlapping, mimicking the cytoskeleton of cells. We first confirm multicolour imaging with fluorophores that range from green to (far infra-) red emission and verify that different photobleaching and blinking kinetics do not impair multicolour imaging as long as cumulant analysis is appropriately performed. Last, we show that our concept is also able to discriminate between three fluorophores with largely overlapping emission spectra. We

Fluorophores with different photobleaching kinetics
To study the influence of differences in photobleaching, we considered Alexa Fluor 488, Atto 565 and Alexa Fluor 647 as above and only change the photobleaching time to 10, 40 and 80 s, respectively (see Supplementary Figure 6). This already covers almost one order of magnitude difference in photostability.
There is no noticeable change in the performance of our analysis, as the blinking kinetics and cumulant analysis are still appropriate for the fluorophore with the worst stability.

Fluorophores with different blinking kinetics
Since not all fluorophores show the same blinking performance under identical experimental conditions 6 , we tested the algorithm with different blinking off-times. As in the first simulations above, we considered

Fluorophores with different brightness
Similarly, it is difficult to achieve equal brightness for all fluorophores in experiments due to different spectral properties and blinking behaviour. As in the first simulations above, we considered Alexa Fluor

a) and b) Average of 4000 frames in the reflection and transmission channel. c) Composite image of the balanced second-order SOFI images with d) Alexa Fluor 488 (cyan hot), e) Atto 565 (yellow hot) and f) Alexa Fluor 647 (magenta hot).
Fluorophores with largely overlapping spectra Our multicolor approach is based on cumulant analysis and we used the same theoretical framework originally devised for spatially super-resolved SOFI 1 . The theoretically attainable resolution increase for second order analysis with respect to widefield imaging is two-fold after spectral cross-cumulation and post-processing (using deconvolution and linearization 2 (used here) or Fourier reweighing 3 ).
We estimated the resolution of our imaging using our implementation of image decorrelation analysis in ImageJ (default settings) 8 . This approach uses partial phase autocorrelation for a series of filtered images to determine the highest spatial frequency with sufficiently high signal in relation to noise. For the fixed cells in Figure 3