Moments reconstruction and local dynamic range compression of high order Superresolution Optical Fluctuation Imaging

Super-resolution Optical Fluctuation Imaging (SOFI) offers a simple and affordable alternative to other super-resolution (SR) imaging techniques. The theoretical resolution enhancement of SOFI scales linearly with the cumulants’ order, while imaging conditions are less phototoxic to living samples as compared to other SR methods. High order SOFI could, therefore, be a method of choice for dynamic live cell imaging. However, due to cusp-artifacts and to dynamic range expansion of pixel intensities, this promise has not been materialized as of yet. Here we investigated and compared high order moments vs. high order cumulants SOFI reconstructions. We demonstrate that even-order moments reconstructions are intrinsically free of cusp artifacts, allowing for a subsequent deconvolution operation to be performed, hence enhancing the resolution even further. High order moments reconstructions performance was examined for various (simulated) conditions and applied to (experimental) imaging of QD labeled microtubules in fixed cells, and actin stress fiber dynamics in live cells.

compared the simulated data results with results from alternative methods, including (1) bSOFI, which utilizes the balanced cumulants to correct for the expanded dynamic range of pixel intensities of high-order SOFI cumulants, and (2) super-resolution radial fluctuation (SRRF) [42], which achieves super-resolution by calculating the degree of local gradient convergence (referred to as the radiality [42]) for each frame and characterize the corresponding temporal radial fluctuations using either temporal radiality average (TRA), or temporal radiality auto-cumulant (TRAC) of different orders. In section 5 we present 6th order moment reconstructions for experimental data together with deconvolution and ldrc. The data sets include quantum dot (QDs)-labeled microtubules in fixed cells and fluorescence protein-labeled β-actin in live cells. Our results are then compared to results obtained by operating the bSOFI [34] and SRRF [42] algorithms to the same data sets. A concluding discussion is given in section 6, summarizing our main findings: (I) even-order moments reconstruction is intrinsically free of cusp artifacts; (II) it can be independently combined with deconvolution without conflicting with the commonly used positivity constraint in image deconvolution; and (III) application of ldrc can correct for the expanded dynamic range of pixel intensities. These attributes allow for SR reconstruction of fast (~seconds) morphological changes in live cells.

Review of SOFI, correlations, cumulants, and moments
We briefly repeat here the SOFI theory [15] but re-cast it in a form that affords the virtual emitter interpretation of SOFI at high orders that we proposed in our accompanying manuscript [33]. This re-casting provides insight into high order SOFI cumulants and the proposed moments reconstruction.
In the practice of SOFI, the sample is labeled with stochastically blinking emitters. This labeled sample is then imaged and consecutive frames are recorded. The data set is then SOFI processed to yield the SOFI image. Given a sample with N emitters that independently blink, the fluorescence signal captured at location r  and time t is given by: where k is the index of the emitter, ϵ k is the 'on'-state brightness of the k th emitter, b k (t) is the stochastic time dependent blinking profile of k th emitter where: 1 when the emitter is in the ''on'' state ( ) , 0 when the emitter is in the ''off'' state  is the point-spread-function (PSF) of the imaging system, and k r  is the location of the k th emitter. In SOFI calculations, we calculate the correlation functions along the time axis with time lags ( 1 2 , ,..., n τ τ τ ) and pixel locations 1 2 ( , ,.... ) n r r r    : 1 (Fig. 1 (a)), from which all the possible partitions are identified as shown in Fig. 1(b). Partitions can possess different numbers of parts where each part can possess different numbers of elements (1st and 2nd columns in Fig. 1 (3 rd column in Fig. 1(b)). Each contributes one term to a summation series to construct the joint-cumulant, where each term can be expressed as the product of two factors. This is shown in the 4 th and 5 th column in Fig. 1(b). The first factor f 1 depends on the size of this partition (denoted as q in 1 st column in Fig. 1(b)) and is defined as: Fig. 1(b)). The second factor f 2 is the product of all the joint-moments of each part within this partition, as illustrated in the 5 th column in Fig. 1 Note here that in Eq. (2.5) above, the joint-moments G(I p ) are essentially the lower order correlation functions discussed in the original SOFI paper [15]. If a partition contains a part that has only one element, we have the corresponding G(I p ) as ( ) 0 As a result, the corresponding f 2 factor will be 0, and this partition will not contribute to the joint-cumulant. The calculation of C 5 (I) is shown in Fig. S1 Detailed derivation of Eq. (2.6) can be found in Appendix 1 [45]. Once the distance factor is solved and divided from both sides of Eq. (2.6), the cumulant value at location gc r  is obtained.
The SOFI pixel location vector is equivalent to the vector average of the selected pixels' locations (in case of pixel repetitions, repeat the corresponding location vectors as well). The choice of pixel combination imposes a trade-off between noise contribution and the attenuation imposed by the distance factor 1 ( ,..., ) n n W r r =   (defined in (2.6)). On one hand, noise could potentially contribute to the resultant cumulant value if there is pixel repetition in the selection. On the other hand, when the selected pixels are distributed too far away from each other, the distance factor becomes small and attenuates the correlation signal. Existing approaches have been focused on avoiding the noise contribution from duplicated pixels [46], but here we explore and present the opposite of this trade off, where we want to diminish the effect of the distance factor at the cost of potential noise contribution. A detailed explanation for our choice of pixel combinations for high order SOFI is given in Appendix 2 and Fig. S2 [45]. Under the framework of virtual emitter interpretation [33], the physical meaning of the joint-cumulant calculated for a set of pixels (either with or without pixel repetition) is taken to mean as the image formed by virtual emitters at the locations of the original emitters, but having virtual brightnesses. These virtual brightnesses are the products of ϵ n (meaning the n th power of the original 'on-state' brightness of the emitter) and w n (δb k (t)) (meaning the n th order cumulant of the blinking profile of the corresponding emitter with the time lags defined for the overall joint-cumulant function). Because the blinking statistics of emitters across the image are not necessarily spatially uniform, the 'on-time ratio', defined as the percentage of time the emitter spent at 'on' state, can vary, causing cumulant values to have different signs at different parts of the image (Fig. 2). Since images are usually presented with positive pixel values, the absolute value operator could yield an image with cusp-artifacts, degrading the image quality of high-order SOFI cumulants [33].

High-order moments reconstruction: theory and Interpretation
Inspired by the interchangeable relation between cumulant and moment [35], we investigated the statistical behavior of high-order moments of emitter blinking trajectories expressed as a function of the 'on time ratio' ρ in a similar way to cumulant analysis [33]. Considering only the blinking profile (with unit brightness) as shown in Eq. (2.2), the 'on' state signal is 1, the 'off' state signal is 0, the time average of the blinking trajectory is ρ, therefore, after subtraction of the time average ρ, the blinking trajectory exhibits (1-ρ) and (-ρ) for the 'on' and 'off' states respectively. We call this signal after subtraction of the time average as the 'center-shifted' signals. Additionally, the percentage of 'on' and 'off 'states in the overall trajectory is ρ and (1-ρ) respectively, providing the weighting factors for both states when calculating the moments. The n th order moment can be readily calculated as the weighted summation of the n th power of the 'center-shifted' signal for both states: Figure 2 shows moments of different orders as a function of ρ ( Fig. 2(a)) in comparison to cumulants of different orders as a function of ρ ( Fig. 2(b)). While cumulants exhibit oscillation between positive and negative values, even-order moments have pure positive values (and odd-order moments are bi-modal with a single node). In practice sampled, ther [33], leading therefore elim reconstructing reconstruction signals from eliminating c (originated fr pure sign (pu deconvolution resolution. In  where W is the 'emitter distance factor', whose analytical form is the same with that of the distance factor [15], and is dependent on the mutual distances between different pixels as shown in (2.6). Detailed derivations of Eq. (3.3) is given in Appendix 3 [45]. We also define m r   as:  From Eq. (3.6) we deduce that this portion of the signal (M1 n ) is equivalent to an image formed by virtual emitters that are located at the same locations as the original emitters with brightnesses ϵ k n M n (δb k (t)) (for the k th virtual emitter). These brightnesses differ from the ones derived for cumulants [33]: ϵ k n C n (δb k (t)). For the k th virtual emitter, its virtual brightness is the product between the n th power of its on-state brightness ϵ k n multiplied by the n th order moment (instead of cumulant) of its blinking fluctuation δb k (t). Because ϵ k n is always positive, even order moments are always positive, therefore the virtual brightness for this portion of the moments signal are always positive.
The second, non-physical part of the summation series in Eq. (3.3) is the case where the partitions contains non-identical emitter location vectors. The corresponding virtual emitters are located at locations where there are no real emitters (unless by coincidence the mass center overlaps with the location of a real emitter with rare chances). This part represents additional virtual (artificial) emitters at locations vectors that are not identical. It originates from cross-ter emitters at ne 'ghost'-emitte distance facto distance facto

High-orde
To take a clo emitters, we s of the simulat tabulated in F the 6 th order m Eq. (3.3) from Fig. 3(d) with emitters. We camera's inte introduced to the scope of th Besides, E brightnesses o in between of [45], where t emitters' inten [45]), indicati critical distan placed in para total number between the a order momen attenuated be not much less Nonetheless, S7 [45]. The t pixel sizes an rms of signals ew locations (d ers). The brigh or, ranging from or [15].

er moments
ose-look of the simulated 3 ne ted movie ( Fig  Fig. 3(a) Fig. 3(b). In which is calcula ment also conf f the PSFs of ning introduce inning effect [35] [45]. All reconstructions were compared to the ground truth image. As shown in Fig. 4, bSOFI reconstructions exhibits discontinuities in the simulated filaments while SRRF artificially narrows them down. moments reconstructions yield a more faithful representation of the simulated data as compared to the ground truth. Further reconstructions results for a variety of simulated challenging image conditions are summarized as supplementary figures in the appendix [45], including for different labeling density (Fig. S8 [45]), increased filaments thickness (or equivalently labeling uncertainty) (Fig. S9 [45]), increased nonspecific binding emitter density (Fig. S10 [45]), and various signal levels (Fig. S11 [45]).
Details of the simulations are given in Appendix 5 [45]. We further tested the 3D sectioning capability on an additional set of simulations where acquisitions of the same simulated sample at 100 different focal planes were generated [37] and processed independently and subsequently combined for 3D reconstruction. ldrc together with moments reconstruction have yielded better sectioning performance than SRRF when compared to the ground truth of the simulation (Fig. 5, Visualization 1, and Visualization 2).

High-order moments reconstruction of experimental data
High-order moments reconstruction (6 th order) in combination with ldrc and deconvolution were applied to experimental data of quantum dots-labeled α-tubulin filaments in fixed Hela cells. The results are compared to bSOFI and SRRF results (Fig. 6). As shown already in the previous section, SRRF exhibit the highest visual resolution enhancement, but at the expense of introduction of distortions, while ldrc-M6 exhibits more faithful results (to the average image). The ldrc-M was performe (Fig. S16  application of deconvolution to the reconstruction, independent of the dynamic range compression process. The theoretical resolution enhancement factor for even order moments is at least 2 , and once combined with subsequent deconvolution algorithm, such factor can be improved to 2n . Lastly, we have demonstrated a super-resolved M6-reconstructed live cell movie with a temporal resolution of 6 seconds per frame (requiring only 200 frames of the original movie for each frame of the reconstructed movie) using a conventional wide field fluorescence microscope.