## Abstract

Super-resolution optical fluctuation imaging (SOFI) is a well-known super-resolution technique appreciated for its versatility and broad applicability. However, even though an extended theoretical description is available, it is still not fully understood how the interplay between different experimental parameters influences the quality of a SOFI image. We investigated the relationship between five experimental parameters (measurement time, on-time *t*_{on}, off-time *t*_{off}, probe brightness, and out of focus background) and the quality of the super-resolved images they yielded, expressed as Signal to Noise Ratio (SNR). Empirical relationships were modeled for second- and third-order SOFI using data simulated according to a D-Optimal design of experiments, which is an *ad-hoc* design built to reduce the experimental load when the total number of trials to be conducted becomes too high for practical applications. This approach proves to be more reliable and efficient for parameter optimization compared to the more classical parameter by parameter approach. Our results indicate that the best image quality is achieved for the fastest emitter blinking (lowest *t*_{on} and *t*_{off}), lowest background level, and the highest measurement duration, while the brightness variation does not affect the quality in a statistically significant way within the investigated range. However, when the ranges spanned by the parameters are constrained, a different set of optimal conditions may arise. For example, for second-order SOFI, we identified situations in which the increase of *t*_{off} can be beneficial to SNR, such as when the measurement duration is long enough. In general, optimal values of *t*_{on} and *t*_{off} have been found to be highly dependent from each other and from the measurement duration.

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

## 1. Introduction

Super-resolution optical fluctuation imaging (SOFI) is a data post-processing technique for fluorescence microscopy which takes advantage of fluctuations in the fluorophore emission to achieve a sub-diffraction resolution [1]. Key features of SOFI are its versatility, low computational demand and its fully quantitative model [2–5]. Its main requirement is that the fluorophores independently cycle between optically distinguishable states, which could be two molecular states that differ in fluorescence brightness, but the emission dynamics can also arise from other mechanisms such as changes in emission spectrum or polarization. Because each fluorescent molecule behaves independently, its contribution can be distinguished from that of other molecules by correlating the observed emission traces, which can be achieved by an analysis based on statistical cumulants. Overall, these relatively mild requirements render SOFI capable of coping with data recorded on any sufficiently sensitive microscope. The method can also be easily extended to 3D applications [6], a variety of fluorophores [7,8], and its effectiveness at short acquisition times makes it well suited for achieving the temporal resolution necessary for live cell imaging [4].

When compared with single molecule localization techniques like photo-activated localization microscopy (PALM [9]) and stochastic optical reconstruction microscopy (STORM [10]), SOFI attains a generally smaller spatial resolution increase [3]. However, the technique can work in situations where single-molecule localization algorithms collapse, such as when fluorophores display low photon counts or when the fluorophore distribution cannot be made sufficiently sparse, a situation that is difficult to handle for single-molecule localization but can also be mitigated using dedicated high-density techniques such as SPIDER [11,12] or Deep-STORM [13]. High fluorophore densities can easily arise when the blinking behaviour of the fluorophores cannot be finely tuned, or when there is little control over the density of the labeling. In light of all of this, SOFI can be considered a versatile tool that can be applied with minimal adaptations in a wide variety of conditions for super-resolution fluorescence imaging.

One issue that is still to be addressed in SOFI is how the interplay between different experimental parameters (such as the fluorescence image acquisition rate, the nature and kinetics of the blinking behaviour, etc.) influences the quality of the obtained super-resolution image. In fact, measurements under different conditions will yield images of widely different qualities. While the values of these parameters are often limited by the specifics of the samples and probes, in certain cases they can be adjusted by varying illumination conditions or changing the sample preparation. DNA-PAINT, for example, allows the blinking behavior to be tuned by the operator [14] by adjusting the composition of the DNA-dye conjugates that bind transiently to the sample structures. However, since in general the mathematical formulation of cumulants is highly non-linear, even when different experimental conditions can be tested, finding an analytical description relating acquisition parameters to the image quality is not straightforward. Therefore, practitioners usually resort to trial-and-error approaches in order to obtain a SOFI image of good resolution and contrast, separately adjusting individual experimental parameters until a satisfactory combination is found. In spite of its apparent simplicity, this procedure is not only tedious in practice, but also prone to error as it ignores possible interactions among parameters. This could lead to an erroneous individuation of their optimal conditions and also does not allow one to assess the causal relationship between the parameters and the image quality.

A more rigorous way to establish this causal relationship is to perform the parameter optimization using design of experiments (DOE). DOE is a well-established approach in many fields of applied sciences that offers striking advantages over trial-and-error strategies [15–21]. A DOE analysis aids in defining the amount of experiments required to estimate the underlying relationship between the experimental parameters, which are called *factors*, and the output of the measurement. This not only reduces the required effort, but can also take into account possible interactions among the factors [22], resulting in a functional model that describes how the variation of the different factors affects the outcome, such as the quality of a SOFI super-resolution image. A DOE model can describe linear relationships between factors and response, but can also include second- or higher-order terms to account for possible curvatures [22]. Developing a DOE model requires performing experimental trials for different values of the factors, where more trials are required when the number of factors increases. More advanced designs, such as D-Optimal designs, can be built *ad-hoc* to reduce the experimental load when the total number of trials to be conducted becomes too high for practical applications [23].

While DOE and related methodologies offer a clear potential to simplify the optimization of advanced imaging methodologies, so far this approach has seen little use in practice. In this work, we sought to explore the feasibility of using DOE to explore the relationship between the experimental input factors and the SOFI image quality. More precisely, we used a five factor D-Optimal experimental design [24] based on simulated SOFI experiments, since these allow exact control over all the input parameters. The resulting model can be exploited to forecast how SOFI imaging would perform for different combinations of experimental parameters within the boundaries imposed for their allowed ranges of values, and to investigate the interplay between those parameters. More generally, our work supports that DOE and related methologies can be readily applied to better understand the performance of advanced fluorescence imaging.

## 2. Theoretical background

#### 2.1 Super-resolution optical fluctuation imaging (SOFI)

In an ideal scenario, a time-lapse fluorescence microscopy signal recorded in the absence of noise can be mathematically expressed as:

The observed fluorescence $F$ depends on the position $\textbf {r}$ and time $t$. In this equation, $F$ is factored as the sum of the signal contributions of the $M$ emitters present in the two-dimensional field-of-view. The contribution of each emitter is weighted by the average emission intensity, $\epsilon$, and the diffraction limited point spread function, PSF, which can be approximated as a 2D Gaussian. The function $s_j(t)$ denotes the stochastic fluctuations in the emission of the $j$-th emitter. The second order SOFI signal $\kappa _{2}$ can then be written as:

#### 2.2 Signal-to-noise ratio (SNR)

Because SOFI is free of localization bias when its core assumptions are satisfied [5], the quality of a SOFI image can be readily assessed by quantifying the noise on the image. Here we achieved this using the signal-to-noise ratio (SNR) [25,26], which was calculated in the following way: for each set of parameters we simulated 20 different datasets, leading to 20 different SOFI images. The per-pixel SNRs were then calculated as the ratios of the average intensity and standard deviation seen for that particular pixel in the replicate images. The overall SNR was then calculated by averaging the per-pixel SNRs.

#### 2.3 Design of experiments (DOE) and optimal designs

We introduce the DOE methodology using a simple example. We assume that we are measuring the output of a particular experiment, $y$. This output is influenced by two different factors $x_1$ and $x_2$, where the interaction between these two factors also has an influence on $y$. We further assume that both $x_1$ and $x_2$ exhibit quadratic behaviour with respect to $y$. Based on such assumptions, $y$ can be formulated as:

The strategy is to perform experimental trials at specific combinations of $x_1$ and $x_2$. We encode the different factor values $x$ in the so-called design matrix $\mathbf {X}$, that has as many rows as there are trials and as many columns as there are terms in the model. We write the result of each corresponding trial in a vector $\mathbf {y}$. The model can then be written in matrix form as:

and the vector $\mathbf {b}$ containing the full set of model coefficients $\beta$ can then be computed as:DOE aims at efficiently answering the question of how many and which trial experiments should be performed to determine $\mathbf {b}$. First of all, we note that varying one factor at a time (keeping the others fixed) would not allow the assessment of factor interactions. Furthermore, a quadratic model like the one in Eq. (3) requires both $x_1$ and $x_2$ to be sampled at least at three distinct settings or *levels* (indexed as −1, 0 and 1 in Table 1). Therefore, to answer the previous question, experiments should ideally be performed at all possible combinations of three levels of $x_1$ and $x_2$ ($3^2=9$ combinations, as shown in $\mathbf {X_f}$ in Table 1), which defines a so-called *full factorial* DOE. In such a trivial situation, spanning all these combinations does not constitute an excessive burden. Nevertheless, if the number of factors increases, this might become unfeasible: e.g., if 4 factors are considered, a full factorial design would be constituted by $3^4=81$ experiments. For 5 factors, $3^5=243$ experiments would be required.

D-Optimal designs provide a way to tackle an unfeasible number of experiments required by a full design. In Optimal designs, the number of experiments to perform is a custom parameter that can be set by the users according to their resources and time. A D-optimal design identifies the best set of experiments to be conducted under this constraint, so as to minimize the uncertainty associated with the model coefficients by maximizing the determinant $\mathbf {D}=|\mathbf {X^{T}}\mathbf {X}|$. The method operates using an iterative search algorithm, by incrementally changing an initial design matrix **X** and exchanging its rows with rows from the design matrix of the corresponding full factorial design. In doing so, at each step the determinant **D** is increased. This procedure was used in Table 1 to obtain a reduced D-Optimal design matrix $\mathbf {X_o}$ from the full factorial design matrix $\mathbf {X_f}$. Here, the number of trial experiments has been set to six (the minimum needed to estimate a six term model) and therefore the resulting matrix is made by the six experiments from the full-factorial design that allow an estimation of $\mathbf {b}$ with comparable reliability [23].

Investigating the resulting polynomial model allows one to identify if a factor should be increased or decreased in order to maximize the measurement result. If a factor is involved only in first order terms, it will be always linearly related to the output. By contrast, if a factor is involved in terms of higher order, then the polynomial partial derivative with respect to it can indicate an inversion point, such that this factor may need to be increased or decreased depending on the experimental conditions. Finally, interaction terms represent synergies or antagonism between two different factors, depending on the sign of the corresponding regression coefficient.

## 3. Materials and methods

#### 3.1 Simulation software

The simulation software used here has been extensively used in previous investigations [5,27–29]. Briefly, the software simulates movies of fluorescence images from a sample labeled with blinking fluorophores. More in detail:

- 1. A predetermined number of emitters are positioned over the field-of-view according to a user-defined distribution.
- 2. Emitters are allowed to blink stochastically, where the survival times of both the fluorescent and non-fluorescent states are assumed to follow an exponential distribution with time constants $t_{\textrm {on}}$ and $t_{\textrm {off}}$.
- 3. The number of photons emitted by each fluorophore in a particular fluorescence image is determined by summing all of the durations during which it was in the fluorescent state in that image (exposure time 0.02 seconds), and multiplying this by the emitter brightness $\epsilon$.
- 4. The distribution of these photons on the detector is modeled as a 2D Gaussian function with a standard deviation of 80.1 nm based on an assumed numerical aperture of 1.4 and emission wavelength 520 nm. The contribution of this molecule to each detector pixel (optical pixel size of 100 nm) is then given by integrating this distribution over this pixel.
- 5. These signal contributions are calculated for all fluorophores and added to the detector pixels which results in an overall fluorescence image.
- 6. Poisson noise is added to each detector pixel.
- 7. Steps 2 to 6 are repeated for each fluorescence image that is to be simulated.

The reference structure for the simulations carried out in this article is displayed in Fig. 1: a square image of 50 pixels by side, with each pixel being 100 nm wide. It was conceived to be reminiscent of the endoplasmatic reticulum.

#### 3.2 Experimental design

The experimental response that was considered in this work is the SNR of the SOFI images obtained at different settings of five factors, measurement duration $T$, $t_{\textrm {on}}$, $t_{\textrm {off}}$, brightness $\epsilon$ and background $b$. The experimental domain for these parameters was defined based on a preliminary screening design, and is reported (in terms of lower and upper limits of the parameter range) in Table 2.

The conditions explored in the design were conceived to be comparable to DNA-PAINT. As previously mentioned, DNA-PAINT represents a particular case in which the blinking behaviour of the probe can be varied during sample preparation and the effect of different measurement conditions can be clearly seen on the SOFI images. For the sake of illustration, a concrete example of this behaviour is given and discussed in the Supplement 1.

We performed the DOE analysis using a D-optimal design encompassing 45 trial experiments. The results were log-transformed, after adding a constant of one:

Log-transformation is commonly used in regression when non-linear relationships exist between the independent and dependent variables. In such cases, the logarithm allows a linear model to be used while keeping the modeled relationship non-linear. Moreover, the log-transformation compensates for the fact that the SNR distribution was found skewed towards low values, which could result in a biased model if left uncorrected.

We developed two full models for second- and third-order SOFI respectively, as they were expected to behave differently with respect to the kinetic behaviour of the probe [5]. We then eliminated irrelevant terms in the full model by performing an F-test on the results of an initial fit. An F-test is a common statistical test used, in this case, to establish whether taking into account the term under consideration leads to significantly reduced residuals compared to those obtained with a so-called *naive* model that comprises only the intercept term. If it does not then the term does not explain a statistically relevant contribution and can be discarded. The performance of the final model was then estimated using the coefficient of determination $R^2$. This is a common metric used in regression problems to assess how effectively the model replicates the observations. It is calculated according to the following equation.

To avoid that the disparity between the magnitudes of the different factors could result in bias in the model, all factors were scaled by converting them to coded units such that the values ranged from −1 to +1 over their domain [22]. Equation 8 shows the relationship between the coded units and the original units values.

Here $x_i$ is the level of the factor expressed in real units, $z_i$ is the same level expressed in coded unit, $x_1^0$ is the factor level $0$ expressed in real units and $\Delta x$ is half the difference between the highest level and the lowest level expressed in real units.

A set of Matlab routines to produce a D-Optimal design and to reproduce the analysis presented in this work are available at the GitHub repository: https://github.com/dcevoli/doe4sofi

## 4. Results and discussion

The influence of the chosen parameters can be seen directly in experiments using DNA-PAINT [14]. In DNA-PAINT the $t_{\textrm {on}}$ is defined as the average binding time between the two oligonucleotide strands, while $t_{\textrm {off}}$ is the average time passing between two successful bindings. While the $t_{\textrm {on}}$ depends uniquely on the affinity between two strands, the $t_{\textrm {off}}$ can be easily manipulated by the operator by modifying the concentration of dye in solution (and therefore of free strands available for the binding). Hence, DNA-PAINT could constitute a good candidate to test for the effect of modifying $t_{\textrm {off}}$ on the SNR in experimental conditions. An experimental DNA-PAINT dataset is presented in the Supplement 1 for the sake of illustration. What can be seen when such experiments are performed is that changing the $t_{\textrm {off}}$ produces a direct effect on the SNR, which seems to increase when the $t_{\textrm {off}}$ decreases. However, increasing the concentration of the dye in solution not only influences the $t_{\textrm {off}}$, but also increases the amount of stray light. Also, since a quantitative relationship between the $t_{\textrm {off}}$ and dye concentration is not know, it would be difficult to estimate actual values of $t_{\textrm {off}}$ from the experimental data alone. Therefore, relying on experimental DNA-PAINT data for a more rigorous optimization would be extremely difficult: a complete experimental model based on simulations is the most straightforward way to assess the influence of each parameter without confounding effects.

For the simulated datasets, we applied our methodology to both second and third-order SOFI analysis, which resulted in two different polynomial expressions describing the relation between the measurement parameters and the SNR of the resulting SOFI images, within the investigated range for each parameter. For the sake of simplicity, we will concentrate first on the model for third-order SOFI, which is simpler and more readily interpretable.

#### 4.1 Third-order SOFI

For third-order SOFI, the final model contains significant terms for three factors: measurement duration $T$, $t_{\textrm {on}}$ and $t_{\textrm {off}}$. The resulting model has $R^2=0.96$. This model was obtained by fitting a full model of 21 terms and then skimming terms found non significant by the F-Test. The full model is reported in the Supplement 1.

We can now investigate this model, focusing especially on the dependence of the SNR on each parameter. To explain further, a factor can be linearly proportional to the SNR, and in this case the direction of improvement needs to be determined, or the relationship may not be linear, and then more than one optimal condition may be found. To this end partial derivatives can be used to isolate the relationship of each factor to the $y$ output, as in Eqs. (10)–12.

Within the domain under study the partial derivative with respect to $T$ (Eq. (10)) is always positive. This means that the SNR always benefits from an increase in measurement duration, independently from the other experimental parameters. This is reasonable since longer measurements on a stationary sample will allow for better estimations of the cumulant values.

The remaining two partial derivatives (Eqs. (11) and 12) have zeroes within the experimental domain: these identify thresholds for the inversion of the dependency between the factors ($t_{\textrm {on}}$ and $t_{\textrm {off}}$) and SNR. These thresholds are represented as colored lines in Fig. 2(a), which schematically displays the domain subspace relative to $t_{\textrm {on}}$ and $t_{\textrm {off}}$.

The partial derivative with respect to $t_{\textrm {off}}$ (Eq. (11)) is positive for values of $t_{\textrm {on}}$ over 0.30 seconds. This threshold is marked as a vertical pink line in Fig. 2(a). When $t_{\textrm {on}}$ is shorter than 0.30 seconds, in order to increase the SNR it is necessary to decrease the $t_{\textrm {off}}$. On the contrary, when $t_{\textrm {on}}$ is higher than this threshold, then the SNR increases proportionally with $t_{\textrm {off}}$.

The partial derivative with respect to $t_{\textrm {on}}$ (Eq. (12)) carries the dependency on both $T$ and $t_{\textrm {off}}$, resulting in a linear equation with negative slope and whose intercept is dependent on $T$, marked as a green line in Fig. 2(a). This line identifies another threshold between two zones: for values above this threshold (upper-right corner of Figs. 2(a)) the SNR is directly proportional to $t_{\textrm {on}}$, while for values below this threshold the SNR is inversely proportional to $t_{\textrm {on}}$.

To summarize, within the investigated range, the two thresholds identify three possible ranges of factor conditions within which different optimization criteria hold (zones (1) to (3) in Fig. 2). For $t_{\textrm {on}}$ below 0.30, regardless of the $t_{\textrm {off}}$, *i*. *e*. zone (1), the highest SNR is obtained with the shortest values of $t_{\textrm {on}}$ and $t_{\textrm {off}}$. When both values of $t_{\textrm {on}}$ and $t_{\textrm {off}}$ are within zone (2), the best results are obtained with the shortest $t_{\textrm {on}}$ and longest $t_{\textrm {off}}$. By contrast, for $t_{\textrm {on}}$ and $t_{\textrm {off}}$ within zone (3), the best results are obtained with the highest value for both factors. During experimental observations using DNA-PAINT (see Supplement 1) results comparable with the optimal conditions in zone (1) of Fig. 2(a) have been found.

The SNR values predicted by the model for each $t_{\textrm {on}}$ and $t_{\textrm {off}}$ combinations are visualised in Figs. 2(b) and 2(c). Among all the experiments, SNR values for our structure have been reported to vary between 0 and over 10. Results below 0.5 were generally not exploitable. The structure started to become recognizable for results over 0.5, and acceptable images could be found only for SNR over 1. The best looking images were found at SNR values larger than 4. Each different line in Figs. 2(b) and 2(c) represents the trend of the SNR values relative to the range of one of the two interacting factors at different levels of the second (see caption). It can be easily seen that the best overall results are obtained for the shortest values of both $t_{\textrm {on}}$ and $t_{\textrm {off}}$, corresponding to an SNR of 5.25. This representation also highlights how the thresholds marked in Fig. 2(a) relate to the SNR trends. In Fig. 2(b), it can be seen as the threshold between zone (2) and zone (3), marked as a green line, represents a point of minimum for each of the colored trends corresponding to different $t_{\textrm {on}}$ values. In Fig. 2(c) instead the trends taken at levels of $t_{\textrm {on}}$ lower than 0.30 seconds, *i*. *e*. above the pink line, are inversely proportional to $t_{\textrm {off}}$, whereas the trends below the pink line are directly proportional to $t_{\textrm {off}}$.

Figure 3 reports results obtained in different simulated conditions, showcasing how the SOFI images change following the direction of improvement suggested by the model. Figure 3(a) represents the same subspace as in Fig. 2, but as a contour plot. The same information on the SNR present in Figs. 2(b) and 2(c) is now summarized in one graph, where the $x$ and $y$ axes represent the two factors $t_{\textrm {on}}$ and $t_{\textrm {off}}$, while the SNR values are represented with a color code ranging from blue (low SNR values) to yellow (high SNR values). As already noticed, the best overall results are obtained for the shortest $t_{\textrm {on}}$ and $t_{\textrm {off}}$ values (corresponding to the red square and Fig. 3(f)), and therefore shorter $t_{\textrm {on}}$ and $t_{\textrm {off}}$ are always preferable within the investigated range. However, when measuring with a label characterised by a $t_{\textrm {on}}$ value greater than 0.30 seconds, the best results can be obtained with longer $t_{\textrm {off}}$, corresponding to the green square. This can be seen by comparing how in Fig. 3 the quality of the obtained images improve when the $t_{\textrm {on}}$ is kept fixed at 0.39 seconds while the $t_{\textrm {off}}$ increases from 25.30 seconds (Fig. 3(g)) to 480.70 seconds (Fig. 3(c)). On the contrary, if measuring with $t_{\textrm {on}}$ lower than 0.30 seconds, the best images are obtained when measuring with the shortest $t_{\textrm {on}}$ and $t_{\textrm {off}}$, as it can be noted comparing Figs. 3(b)–3(d)–3(f), obtained using decreasing values of $t_{\textrm {on}}$ and $t_{\textrm {off}}$. It can be also noted as Fig. 3(g) is affected by artefacts that may look like cusp-artifacts. These artifacts are due to the positive and negative oscillations in cumulants of third or higher-order SOFI, which result in the separation of the signal from a single emitter in two lobes delimited by a black boundary [30]. In the image obtained in the optimal condition the same artefacts appear to be absent.

#### 4.2 Second-order SOFI

The final model for second-order SOFI (Eq. (13)) contains significant terms for four factors: measurement duration $T$, $t_{\textrm {on}}$, $t_{\textrm {off}}$ and background $b$. It has $R^2$=0.98. The partial derivatives of the model are reported in Eqs. (14)–17. As previously explained, the final model was obtained by fitting a full model containing 21 terms and then skimming the terms found non significant by the F-Test. The full model is reported in the Supplement 1.

In Fig. 4, an overview of the interpretation of the model is presented. In the first column (Figs. 4(a), 4(d), 4(g) and 4(j)), four different sub-spaces are schematically represented, relative to the four significant interactions present in the model. Each corresponding sub-figure in the second and third column represents the trend of the SNR values relative to that interaction, *i*. *e*. the trend of the SNR within the range of one of the two factors, color-coded at different levels of the other. For example, Fig. 4(a) schematizes the $T$-$b$ interaction. Figure 4(b) shows the trend of SNR as a function of $T$ and different values of $b$ (color-coded, see caption). Conversely, in Fig. 4(c), SNR is plotted as a function of $b$ while $T$ is represented by the color-code. For each interaction, some thresholds have been identified by solving Eqs. (14)–17. For example, the threshold for the $T$-$b$ interaction is visible as a pink line in 4(a) and 4(c). As it can be seen in 4(c), the SNR is inversely correlated with background when the measurement time is shorter than 339 frames, and directely correlated to the background when the measurement time is longer than 339 frames. However, these conditions correspond to very low SNR values, up to the point that the image structure is barely recognizable; therefore they are scarcely relevant for optimization purposes.

To summarize, the model obtained for second-order SOFI offers the following conclusions: as in third-order SOFI, the SNR in second-order SOFI is always directly proportional to the measurement duration (Figs. 4(b) and 4(e)). Background is generally inversely proportional to the SNR, except for high $t_{\textrm {on}}$ and low measurement duration, conditions in which the SNR is anyway extremely low. Investigating the interplay between $t_{\textrm {on}}$ and $t_{\textrm {off}}$ reveals a more complex situation, but the best overall SNR values are again obtained when both are at their respective lowest values. More specifically, the optimal SNR value (over 14.5) was reached with the longest measurement duration, the lowest $t_{\textrm {on}}$ and $t_{\textrm {off}}$, and the lowest background. These conditions yielded the super-resolved SOFI image shown in Fig. 5.

## 5. Conclusions

In this work we exploited the principles of DOE to construct an empirical model linking five fluorescence image acquisition parameters to the quality of the SOFI images they yield, as evaluated through the use of the SNR. The model was obtained using a D-Optimal design of experiments.

In our findings, second- and third-order SOFI images showed a different behavior, though both exhibited increased SNR with increasing measurement durations, which is logical since the estimates of the cumulants will be more accurate when stationary samples are measured over a longer duration. Moreover, the results obtained from both orders of SOFI show that the best SNR values were found when $t_{\textrm {on}}$ and $t_{\textrm {off}}$ were minimized. By contrast, we did not find an effect of the emitter brightness within the range of conditions examined here. Our models also revealed nonlinear behavior with respect to the factors $t_{\textrm {on}}$ and $t_{\textrm {off}}$. In fact, both orders exhibit regions within the experimental domain in which a relatively high SNR is found at high values of $t_{\textrm {on}}$ and $t_{\textrm {off}}$. The maximum SNR value reachable at high values of $t_{\textrm {on}}$ and $t_{\textrm {off}}$ is nevertheless not as good as the global optimum.

This work demonstrated the usefulness of a DOE-based framework for the investigation of SOFI, and super-resolution algorithms in general, given its ability to deal with non-linear relationships and interactions among measurements parameters. The methodology showcased can be replicated by operators (tailoring the factors considered and the experimental domain explored) to obtain useful optimization indications or to alleviate the experimental time investment. It could also be equally valuable when trying to gain insights in the behaviour of SOFI from experimental data, or (when coupled with simulated data) to probe its effectiveness in extreme or ideal conditions.

## Funding

International Associate Laboratory for High Performance Fluorescence Microscopy UnivLille-KU Leuven; European Research Council (714688 NanoCellActivity); Fonds Wetenschappelijk Onderzoek (G090819N).

## Acknowledgments

The authors thank Peter Goos (KU Leuven) for fruitful discussions. We thank Ralf Jungmann, Thomas Schlichthärle, and Alexander Auer (Max Planck Institute for Biochemistry) for providing the DNA-PAINT sample. This work was supported by the FLAG-ERA grant SENSEI via FWO grant G0G2819N.

## Disclosures

The authors declare no conflict of interest.

## Supplemental document

See Supplement 1 for supporting content.

## References

**1. **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]

**2. **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]

**3. **S. Geissbuehler, C. Dellagiacoma, and T. Lasser, “Comparison between SOFI and STORM,” Biomed. Opt. Express **2**(3), 408–420 (2011). [CrossRef]

**4. **H. Deschout, T. Lukes, A. Sharipov, D. Szlag, L. Feletti, W. Vandenberg, P. Dedecker, J. Hofkens, M. Leutenegger, T. Lasser, and A. Radenovic, “Complementarity of PALM and SOFI for super-resolution live-cell imaging of focal adhesions,” Nat. Commun. **7**(1), 13693 (2016). [CrossRef]

**5. **W. Vandenberg, M. Leutenegger, S. Duwé, P. Dedecker, and P. Dedecker, “An extended quantitative model for super-resolution optical fluctuation imaging (SOFI),” Opt. Express **27**(18), 25749–25766 (2019). [CrossRef]

**6. **S. Geissbuehler, A. Sharipov, A. Godinat, N. L. Bocchio, P. A. Sandoz, A. Huss, N. A. Jensen, S. Jakobs, J. Enderlein, F. Gisou van der Goot, E. A. Dubikovskaya, T. Lasser, and M. Leutenegger, “Live-cell multiplane three-dimensional super-resolution optical fluctuation imaging,” Nat. Commun. **5**(1), 5830 (2014). [CrossRef]

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

**8. **K. Grußmayer, T. Lukes, T. Lasser, and A. Radenovic, “Self-blinking dyes unlock high-order and multiplane super-resolution optical fluctuation imaging,” ACS Nano **14**(7), 9156–9165 (2020). [CrossRef]

**9. **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]

**10. **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]

**11. **S. Hugelier, J. J. de Rooi, R. Bernex, S. Duwé, O. Devos, M. Sliwa, P. Dedecker, P. H. C. Eilers, and C. Ruckebusch, “Sparse deconvolution of high-density super-resolution images,” Sci. Rep. **6**(1), 21413 (2016). [CrossRef]

**12. **S. Hugelier, P. H. C. Eilers, O. Devos, and C. Ruckebusch, “Improved superresolution microscopy imaging by sparse deconvolution with an interframe penalty,” J. Chemom. **31**(4), e2847 (2017). [CrossRef]

**13. **E. Nehme, L. E. Weiss, T. Michaeli, and Y. Shechtman, “Deep-storm: super-resolution single-molecule microscopy by deep learning,” Optica **5**(4), 458–464 (2018). [CrossRef]

**14. **J. Schnitzbauer, M. T. Strauss, T. Schlichthaerle, F. Schueder, and R. Jungmann, “Super-resolution microscopy with DNA-PAINT,” Nat. Protoc. **12**(6), 1198–1228 (2017). [CrossRef]

**15. **M. Kuterbekov, P. Machillot, F. Baillet, A. M. Jonas, K. Glinel, and C. Picart, “Design of experiments to assess the effect of culture parameters on the osteogenic differentiation of human adipose stromal cells,” Stem Cell Res. Ther. **10**(1), 256 (2019). [CrossRef]

**16. **S. A. Weissman and N. G. Anderson, “Design of experiments (DoE) and process optimization. a review of recent publications,” Org. Process Res. Dev. **19**(11), 1605–1633 (2015). [CrossRef]

**17. **R. Leardi, “Experimental design in chemistry: A tutorial,” Anal. Chim. Acta **652**(1-2), 161–172 (2009). [CrossRef]

**18. **M. A. Bezerra, R. E. Santelli, E. P. Oliveira, L. S. Villar, and L. A. Escaleira, “Response surface methodology (RSM) as a tool for optimization in analytical chemistry,” Talanta **76**(5), 965–977 (2008). [CrossRef]

**19. **E. S. Hecht, A. L. Oberg, and D. C. Muddiman, “Optimizing mass spectrometry analyses: a tailored review on the utility of design of experiments,” J. Am. Soc. Mass Spectrom. **27**(5), 767–785 (2016). [CrossRef]

**20. **M. Mäkelä, “Experimental design and response surface methodology in energy applications: A tutorial review,” Energy Convers. Manage. **151**, 630–640 (2017). [CrossRef]

**21. **P. Yu, M. Y. Low, and W. Zhou, “Design of experiments and regression modelling in food flavour and sensory analysis: A review,” Trends Food Sci. Technol. **71**, 202–215 (2018). [CrossRef]

**22. **A. Dean, D. Voss, and D. Draguljić, * Design and Analysis of Experiments* (Springer Texts in Statistics, 2017).

**23. **P. Goos and B. Jones, * Optimal Design of Experiments: A case Study Approach* (Chichester, 2011). OCLC: 846512756.

**24. **P. F. de Aguiar, B. Bourguignon, M. S. Khots, D. L. Massart, and R. Phan-Than-Luu, “D-optimal designs,” Chemom. Intell. Lab. Syst. **30**(2), 199–210 (1995). [CrossRef]

**25. **B. Moeyaert, W. Vandenberg, P. Dedecker, and P. Dedecker, “SOFIevaluator: a strategy for the quantitative quality assessment of SOFI data,” Biomed. Opt. Express **11**(2), 636–648 (2020). [CrossRef]

**26. **W. Vandenberg, S. Duwé, M. Leutenegger, B. Moeyaert, B. Krajnik, T. Lasser, and P. Dedecker, “Model-free uncertainty estimation in stochastical optical fluctuation imaging (SOFI) leads to a doubled temporal resolution,” Biomed. Opt. Express **7**(2), 467–480 (2016). [CrossRef]

**27. **W. Vandenberg and P. Dedecker, “Effect of probe diffusion on the SOFI imaging accuracy,” Sci. Rep. **7**(1), 44665 (2017). [CrossRef]

**28. **Y. Peeters, W. Vandenberg, S. Duwé, A. Bouwens, T. Lukeš, C. Ruckebusch, T. Lasser, and P. Dedecker, “Correcting for photodestruction in super-resolution optical fluctuation imaging,” Sci. Rep. **7**(1), 10470 (2017). [CrossRef]

**29. **G. C. H. Mo, B. Ross, F. Hertel, P. Manna, X. Yang, E. Greenwald, C. Booth, A. M. Plummer, B. Tenner, Z. Chen, Y. Wang, E. J. Kennedy, P. A. Cole, K. G. Fleming, A. Palmer, R. Jimenez, J. Xiao, P. Dedecker, and J. Zhang, “Genetically encoded biosensors for visualizing live-cell biochemical activity at super-resolution,” Nat. Methods **14**(4), 427–434 (2017). [CrossRef]

**30. **X. Yi and S. Weiss, “Cusp-artifacts in high order superresolution optical fluctuation imaging,” Biomed. Opt. Express **11**(2), 554–570 (2020). [CrossRef]