Simulations of camera-based single-molecule fluorescence experiments

Single-molecule microscopy has become a widely used technique in (bio)physics and (bio)chemistry. A popular implementation is single-molecule Förster Resonance Energy Transfer (smFRET), for which total internal reflection fluorescence microscopy is frequently combined with camera-based detection of surface-immobilized molecules. Camera-based smFRET experiments generate large and complex datasets and several methods for video processing and analysis have been reported. As these algorithms often address similar aspects in video analysis, there is a growing need for standardized comparison. Here, we present a Matlab-based software (MASH-FRET) that allows for the simulation of camera-based smFRET videos, yielding standardized data sets suitable for benchmarking video processing algorithms. The software permits to vary parameters that are relevant in cameras-based smFRET, such as video quality, and the properties of the system under study. Experimental noise is modeled taking into account photon statistics and camera noise. Finally, we survey how video test sets should be designed to evaluate currently available data analysis strategies in camera-based sm fluorescence experiments. We complement our study by pre-optimizing and evaluating spot detection algorithms using our simulated video test sets.


Introduction
Since the initial proof of concept, single-molecule fluorescence techniques, in particular single-molecule Förster resonance energy transfer (smFRET), have proven powerful tools in probing biomolecular structures and dynamics [1][2][3]. Single fluorescent molecule sensitivity is predominantly achieved using two experimental configurations: (i) confocal microscopy in conjunction with single-photon detection (avalanche photodiodes or photomultiplier tubes) [4,5] and (ii) total internal reflection fluorescence [6] or wide-field microscopy with intensitybased detection using either an electron multiplying charge-coupled device (EM-CCD) [7,8] or a scientific complementary metal-oxide-semiconductor (sCMOS) camera [9]. Camera a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 detection is characterized by a time resolution in the lower millisecond range and a spatial resolution reaching the diffraction limit of visible light [7]. Much higher time resolution can be achieved when photon counting detectors are used in time-correlated single photon counting (TCSPC) experiments with pulsed excitation sources. Conversely, time-binned methods are rather straightforward to implement, comparably inexpensive and, importantly, allow for parallel recording of hundreds to thousands of single molecules at the same time [10,11].Consequently, camera-based fluorescence detection of single-molecules has become a widespread approach for high-throughput acquisition of single-molecule data, especially in the context of smFRET experiments.
Camera-based smFRET generates large and complex data sets that are henceforth referred to as single molecule videos (SMVs). SMVs are characterized by a low net signal and a low signal-to-noise ratio (SNR), a high molecule surface density (ρ), short intermolecular distances (IMD), and inhomogeneous background profiles. Much effort has been geared towards SMV data analysis in recent years [11][12][13], but their analysis is at present not standardized. Instead, individual researchers face a host of different data analysis strategies developed by different groups [12,[14][15][16][17]. Therein, SMV simulations in camera-based SM detection have been reported in the context of kinetic/dwell-time analysis, where method-specific kinetic rates are determined and evaluated against their simulated ground truth, usually as a function of the SNR [12,14,18,19]. Although these studies describe the simulation process used, recovering video simulation parameter (VSP) values from the manuscript or supplementary material to reproduce the simulated SMV is typically not straightforward. Independent and reliable assessment of such data analysis strategies requires a common standard to be defined. Thus, the systematic annotation of SMVs with metadata files, would ease and speed up the process of data reproduction. However, standardized, simulated smFRET data test sets annotated with metadata are, to the best of our knowledge, at present not available.
In the field of computer science, well-annotated and independently designed sets of test data allow the evaluation of individual algorithms [20,21]. Thereby, potential user-bias is minimized, and thus, the evaluation of algorithm via test data sets are accepted as common standard. We believe that standardized, simulated SMV are very suitable as test data sets as they inherently fulfill two conditions: First, they allow the reproducible evaluation of processing algorithms used in smFRET data analysis by providing a set of video simulation parameters (VSP) screening a large experimentally relevant parameter space [18]. Second, simulated SMVs do not contain experimental artifacts due to their unambiguous definition of molecular properties and instrument-specific parameters (Fig 1). A comparable approach has been presented by Sage et al. for the evaluation of single molecule localization microscopy methods [22] and by Preus et al. for the evaluation of background estimators in single molecule microscopy [23].
Here, we provide the theoretical groundwork to generate realistic SMVs covering a large, experimentally relevant VSP space. Specifically, the methodology presented herein addresses (i) the design of thermodynamic and kinetic models, (ii) the distribution of single-molecules in the cameras field of view (FOV) and point-spread function (PSF) modeling using a 2D-Gaussian approximation, (iii) the simulation of photobleaching, (iv) the introduction of signal contributions like background, spectral bleed-through and camera noise and (iv) the SMV export with respect to output file formats and SMV documentation. We present the concepts for performance evaluation of algorithms addressing multiple aspects of SMV data analysis, including SM localization, background correction, intensity time trace generation, trajectory discretization, as well as kinetic and thermodynamic modeling. We use our MATLAB-based SMV simulation tool within our Multifunctional Analysis Software for Handling single molecule FRET data (MASH-FRET) [13], which is easy-to-use and permits rapid See the README.md for the respective download links. All results for the method parameter optimization and the subsequent comparison of the evaluated SM localization algorithms are presented on an interactive web page. See the README.md for the respective web links. All other data are within the main manuscript and its supporting information files. and integrated analysis of camera-based smFRET experiments, to simulate well-annotated SMV test data sets. We apply them to optimize method parameters (MP) of different SM localization algorithms commonly used in the field of camera-based SM fluorescence and compare their performance as a function of the total emitted intensity of single molecules. Both MASH-FRET and the benchmarking data sets are freely available via https://github.com/ RNA-FRETools/MASH-FRET.

The principles of FRET and camera-based smFRET experiments
Förster resonance energy transfer (FRET) denotes the dipole-dipole coupling between a donor fluorophore and an acceptor that is usually, but not necessarily [24,25], also a fluorophore. According to Förster's theory developed in the 1940s, the efficiency of the process, i.e., the FRET efficiency or transfer efficiency E, depends on the inverse sixth power of the interdye distance r [26,27] where R 0 denotes the Förster radius, i.e., the distance specific to the FRET pair resulting in a FRET efficiency of 50%. Owing to the pronounced distance dependence of the FRET efficiency, FRET is frequently referred to as a spectroscopic ruler within a range of 3 to 10 nm [28]. For further information on FRET, please refer to dedicated reviews [3,6,15,29]. Camera-based smFRET typically involves total internal reflection (TIR) excitation of surface-tethered biomolecules that are fluorophore labeled [7]. Total internal reflection (TIR) of an incident laser beam can be achieved with a prism or directly through a suitable objective, Simulated SMVs can be applied to evaluate data analysis algorithms (left) and for comparison with experimental data (right). In the former application, MC simulations are performed to generate a set of SMVs characterized by defined VSPs, followed by analysis using the algorithm to be evaluated. The results (output parameters of the method) are then compared to the input parameters of the simulation to quantify the performance of the algorithm. Method specific parameters (MPs) are varied to maximize the agreement between input and output parameters in order to reach maximum accuracy and efficiency (see Section 4 for further details). Using preoptimized parameter sets to analyze experimental data (application) yields reliable results of the molecular system under study.
yielding an evanescent field with an exponentially decaying intensity profile for fluorophore excitation and low out-of-focus fluorescence [6]. Emitted photons are sorted according to their wavelength and detected using a camera. These raw videos are referred to as single-molecule videos (SMVs) throughout this article. In SMVs, fluorophore-labeled biomolecules act as point emitters of light and appear as PSFs in the FOV. In a perfect optical system without spherical aberration, their diffraction pattern can be satisfactorily described by a symmetric 2D Gaussian characterized by the full width at half maximum (FWHM) w det,D/A,0 [30].
Extracting donor (D) and acceptor (A) emission rates from SMVs permits to approximate the FRET efficiency by the apparent time-dependent transfer efficiency FRET(t) [7,31]  D; ex correspond to the photon emission of the acceptor and donor fluorophore at a time point t, respectively. In practice, the detected photon emission rates are affected by spectral bleed-through bt D/A , which results from the detection of donor photons in the acceptor channel, as well as direct acceptor excitation dE A . Both bt D/A and dE A can be determined from the overlap of the emission and absorption spectra of the fluorophores taking into account set of optical filters used, although they are generally determined experimentally [7]. Furthermore, the detected photon emission rates contain background signal with different physical and technical origin (see below for further details). Differences in donor and acceptor quantum yields QY D and QY A , as well as detection probabilities of both donor η D and acceptor η A are accounted for by the correction factor γ. A detailed overview on the most relevant corrections and their effects is given in the Supporting Information (Section C in S1 File). Applying these corrections eventually yields absolute FRET efficiencies E(r,t) FRET abs (t), which can be converted into mean inter-dye distances [31][32][33][34].

Simulation of SMVs
FRET efficiencies are calculated from fluorescence intensities. Consequently, we seek to simulate SMVs that contain time-binned fluorescence intensity trajectories, in which the stochastic nature of photon emission is accounted for by a Poisson probability distribution. Simulating SMVs requires the following VSPs to be defined (Fig 2): (i) The number of FRET states and the corresponding mean FRET efficiencies adopted by the molecular system under study, (ii) a kinetic model, which describes the transition rates between these FRET states as a Markovian process (Section 3.2), (iii) the photophysical parameters of the FRET pair, (iv) the number and location of single molecules within the FOV (Section 3.3.2), (v) the PSF model defined by the imaging system (Section 3.3.3), (vi) the background signal (Section 3.3.4), as well as (vii) a model to describe the noise of the camera (Section 3.3.6). The graphical user interface (GUI) of our MATLAB-based SMV simulation tool shown in Schematic A in S1 File allows the definition of all the above-mentioned VSPs in a straightforward manner. The simulation tool permits the user to export all video simulation parameters in Matlab-independent exchange file formats.

Workflow to simulate single molecule FRET data
In order to obtain SMVs and FRET efficiency time traces that closely match experimental data, kinetic Monte Carlo (MC) simulations were used to generate N fluorescence intensity trajectories (Fig 2) [35]. For this purpose, we defined the total number of conformational states (molecular states) J as well as the corresponding J(J-1) monoexponential rate constants k ij characterizing the interconversion between states i and j. Each state j was assigned to a discrete mean FRET efficiency FRET j . Simulation of each time trace n and frame l in a time series characterized by a total length L, followed by generating SMVs, was achieved in eight steps: 1. Based on the kinetic rates k ij , KMC simulations were performed to determine FRET(n,l) of each molecule n for each time bin l, i.e., the observation time a molecule spends in a certain molecular state. In the case of dwell times shorter than the time bin l, FRET(n,l) values of multiple states were averaged within the same bin l. Kinetic heterogeneity, i.e., state transitions characterized by multi-exponential kinetics [36], can be simulated by assigning different rate constants k ij with the same FRET j value but different states j. The kinetic model is further described in Section 3.2.
2. Donor and acceptor fluorescence intensities Iðn; lÞ D; em D; ex and Iðn; lÞ A; em D; ex were calculated according to Eq (2) based on the discrete FRET(n,l) values and assuming a total emitted donor fluorescence I tot,0 in the absence of FRET. Even though a constant I tot,0 for all SMs was used as default setting, it is possible to define molecule-to-molecule variations in I(n) tot,0 to model non-uniform illumination or to create test data sets featuring multiple SNR values (Section 3.3.3). Further, the probability of donor excitation in absence of FRET was assumed to be independent of time and the relative dipole orientation of both fluorophores. Additionally, the exponential decay in the z-direction of the evanescent wave generated in TIR was neglected, as immobilization typically confines the biomolecules in close proximity to the surface. Differences between quantum yields QY D and QY A , as well as between the detection efficiencies η D and η A were accounted for by the correction factor γ (Section C in S1 File). It is noteworthy, that the gamma correction introduces an indirect change of the total emitted fluorescence (Table C in S1 File) independent of the user-specified distribution I(n) tot,0 introduced above. A frequently observed phenomenon in smFRET are molecule-to-molecule variations. For example, variations can be observed with regard to the mean FRET value of a certain conformational state, the total emitted intensity and quantum yields [3,37]. We modeled cross-sample variability assuming a Gaussian distribution of the underlying VSP values characterized by a defined center and standard deviation σ, thus, FRET j and σ FRET,j , I tot,0 and σ Itot,0 , and γ and σ γ . This straightforward approach allows to simulate cross-sample variability in FRET states and trajectory SNR originating from instrumental artefacts, but also fluctuations of QY, interdye distances, local refractive indices, and fluorophore orientations that lead to κ 2 variations [18,38,39].
3. Instrumental imperfections resulting in direct acceptor excitation and spectral bleedthrough were also considered. The respective correction factors are listed in Section C in S1 File [40,41].
4. To simulate spontaneous dye photobleaching (see option in Schematic A in S1 File), trajectories were truncated using an exponentially distributed photobleaching time. 5. N molecules were distributed over a virtual FOV in the respective donor and acceptor channel (Section 3.3.1) using either pre-defined or random coordinates (Section 3.3.2). The position of each molecule (x n ,y n ) is time-invariant, i.e., focal drift was assumed to be absent. 6. Pixel values of the respective single molecule coordinates were set to the calculated donor and acceptor fluorescence intensities I(x n ,y n ,l) D/A in each frame and subsequently convoluted with the PSF of the virtual imaging system as discussed in Section 3.3.3. A spatially and/or temporally variable background was added to each frame to account for different sources of background as discussed in Section 3.3.4 and Section 3.3.5. Instrumental imperfections such as focal drift and chromatic aberrations can be simulated too.
7. Photon emission from fluorescent single molecules and background sources was implemented as all-or-none process with a constant probability to emit a photon according to Poisson statistics. In our MC simulations, the resulting photon shot-noise in donor and acceptor emission intensities was accounted for with a Poisson distribution of pixel intensities centered around the expected mean fluorescence intensity (Section 3.3.6).
8. The detected photon signal was finally convoluted with a realistic SNR function taking into account the sensitivity, linearity, and temporal noise of camera sensors according to "Emva 1288 Standard Release 3.1" [42] and Hirsch et al. [43]. (Section 3.3.6). All simulation parameters were set to default values inspired by values typically observed in an experimental setting ( Table 1).

Describing single molecule FRET trajectories as Markov chains
smFRET experiments rely on the assumption that a molecular system capable of adopting J conformations yields J unambiguously discernable FRET efficiencies. Therefore, we assigned each conformation j to a discrete FRET efficiency value FRET j in our simulations. Here, the thermodynamic equilibrium between these states depends on the kinetic rates k ij that characterize interconversion between state i and j for all states J. This interconversion was modeled as a Markov chain (Fig 2) [46]. In a Markov chain, state transitions are treated as homogenous processes, and the respective J × J matrix of transition probabilities p ij to transit from a state i to j is defined by the transition rates k i6 ¼j and the frame rate f according to [47]: where p i = j determines the probability to stay within the same state and p i6 ¼j to transit from state i to state j, respectively. For i 6 ¼ j, rates were assumed to adopt values ! 0, whereas rates were set to 0 for i = j. In the context of this study and for each simulated time trace n, we defined the state i of time bin l = 0 as the state with the highest probability p i = j according to the probability matrix in Eq (3). We assumed a total length L of the observation time and determined the conformational state at each time bin by means of a KMC simulation. Please note that Matlab's built-in function for HMM time series generation (hmmgenerate) parametrized using p jj' in Eq (3) to build a time-resolved kinetic model and to generate the corresponding FRET time trace does not allow for time averaged states of single bins or frames, respectively. The simulation of the respective states in the following time bins can be carried out in two different ways: (i) bin wise or (ii) dwell time wise. For (i) we generated a random value from a uniform probability distribution between 0 and 1 and compared it to the cumulative probability vector P i ðjÞ ¼ X j n¼1 p in of the respective state i in time bin l = 0. Thus, we defined the new state in time bin l = 1. The simulation was repeated for each time trace n until time bin L was reached. For (ii) we set p i = j = 0, generated a random value and compared it to the probability vector P i ðjÞ ¼ X j n¼1 p in . Thus, we obtained the next state j and generated the dwell time in state I from a monoexponential probability distribution, a hallmark of classical chemical kinetics, characterized by the respective rate coefficient k ij . All time bins within the generated dwell time were set to state i. In the case of dwell times exceeding the last time bin, the last dwell time was truncated to the total length L. Note that in particular transitions rates k ij ! f lead to the detection of artificial states due to time averaging [29]. Therefore, in the case of dwell times shorter than the binning time, we built a time-weighted average of consecutive states in time bin l and l+1 yielding an averaged FRET value (artificial states) typical for kinetics faster than the exposure time of the camera. Both approaches objectively define the conformational state j of the (bio)macromolecule under study n for each single time bin l. The number of simulation steps in the second approach is drastically reduced, which allows faster simulations. Therefore, it is recommended for the creation of long time traces, traces with small transition rates and traces with large differences between transition rates in order to perform a statistically meaningful kinetic analysis later on.
In intensity-based smFRET experiments, kinetic rates, and thus the transition probabilities, are experimentally accessible from dwell-time analysis [3] or by maximizing the likelihood of a matrix of rate constants in Eq (3) according to Gopich and Szabo [48]. Thus, experimentally obtained rate constants can be used to simulate the molecular system under study for comparison, an approach which was presented recently [24].

Instrument specific configuration
In smFRET, the detected signal originates from fluorescent molecules that are only a few nanometers in size. This results in a diffraction-limited spot, whose intensity profile is determined by the characteristics of the optical imaging system, including the microscope objective and beam expander that is typically placed in the detection pathway for magnification (please see Table C in S1 File for further examples). In the following, we describe how we accounted for these instrument specific parameters with an appropriate PSF model, by adjusting the number of single molecules in the FOV, the background signal and the noise model of the camera. All of these VSPs can be independently defined in our MATLAB GUI shown in Schematic A in S1 File.

Camera and video parameters.
SMVs directly dependent on camera-specific parameters, such as the spatial resolution, the pixel size, and the video length, all of which are independent from the actual molecular system. We defined the image size or FOV by the pixel size d pix and the dimensions (x res ,y res ) of the camera sensor. Here, the latter depends on the number of horizontal and vertical divisions on the sensor, x res and y res , i.e., the number of pixel within a Cartesian coordinate system (x, y). State-of-the-art EMCCD cameras such as an the Andor iXon3 897D (Oxford Instruments, UK) are characterized by a pixel edge width of 16 μm and an x res ×y res camera dimension of 512×512 pixel. Note that within an x res ×y res image matrix, x res defines the number of columns and y res the number of rows following Matlab (and other programming languages) notation where an image matrix is defined as lines by rows, i.e. y res ×x res . This corresponds to a total sensor area of 8190×8190 μm 2 . In a widefield microscope with 150-fold magnification (objective + beam expander, Table C in S1 File), the total sensor area corresponds to an area of 54×54 μm 2 in object space. It should be noted that p x p pixels can be combined to one superpixel (hardware binning), an approach that increases the SNR at the expense of decreased spatial resolution. In this study, we performed hardware binning of 2 x 2 pixel and split the detector into two spectroscopic channels characterized by a resolution of 256×128 pixel each. Thus, the simulated object area 54×27 μm 2 in each spectroscopic channel yields an object size of 0.21×0.21 μm 2 of the virtual camera pixel in image space perfectly matching the diffraction limit of visible light (Table C in S1 File).
The total number of frames L determines the length of the video. The frame rate f is a rather arbitrary value in simulated SMV, since greyscale videos are usually saved as x res × y res × L matrices. Here, each x res × y res image was associated with a binning interval, typically between 10 ms and 100 ms. The binning interval is important for visualizing the SMV with a video player and it affects the outcome of the kinetic analysis (Section 3.2).
Photoelectrons generated in response to incident photons are digitized, yielding an integer data type with a given bit rate (BR). The BR defines the range of values used for the assignment of pixel intensities in units of image counts, electron counts or photon counts. For example, at a BR of 14 bit per pixel, the intensity adopts values from 0 to 2 14 −1 = 16383 counts. It should be noted, however, that the BR defines a saturated maximum of pixel intensity, which is defined as the maximum of the measured relation between the variance of the gray value and the incident photon flux in units photons/pixel [42]. Thus, the maximum detectable photon flux is given by the variance of the respective image count value of the camera.
We simulated fluorescence intensity trajectories in units of photon counts (pc) per time bin (pc/time bin) or per frame (pc/frame). This enables the evaluation of methods with SMV test data sets independent from the camera specific image count signal. The MASH-FRET analyzing software does not require pc as unit for the intensity; units are arbitrary and can be pc or any other measure for the intensity. However, in the case where pc conversion is not provided by the acquisition software of the camera, i.e., the detected fluorescence intensity is given in units of image counts (ic) or electron counts (ec) per bin time (ec/time bin = cps) or per frame (ec/frame = cpf), individual measures of the camera should be converted into photon [49]. Equations for converting photon counts into image counts and vice versa are provided in Section 3.3.6 or [42].
3.3.2 Distributing single molecules within the FOV. We used two approaches to distribute N simulated molecules within one half of the FOV: (i) We positioned molecules at predefined pixel locations (Fig 3A), yielding regular patterns similar to the ones obtained when nanostructured materials like zero-mode waveguides are used for surface-immobilization [50,51]. Here, overlap between individual PSFs is absent. (ii) We randomly distributed molecules (Fig 3B), typically leading to overlapping PSFs as observed when surface immobilization is achieved via biotinylated PEG or BSA [36,52].
To obtain a single molecule density ρ, which is defined as the number of SMs per area, and to avoid overlap of different spots, the number of simulated molecules N was adjusted as a function of the image dimensions, i.e., the FOV. To simulate a typical single molecule experiment, we used a low surface molecule density (0.035 μm -2 ρ 0.063 μm -2 ). On average, this corresponds to 82 N 146 molecules distributed over the simulated FOV of 256 x 128 pixel (Fig 3B). Please note that we used subpixel accuracy to distribute molecules within the FOV, resulting in different appearance of single molecules in SMVs depending on their localization relative to the grid of pixels (Fig 3A). Subpixel accuracy in simulated SMVs is particularly important for the evaluation of spot detection algorithms in super resolution microscopy techniques [16,53].

PSF model.
The intensity profile of a diffraction-limited spot within the conjugated image plane using an ideal imaging system with a high numerical aperture (NA > 1) can be approximated with a normalized, symmetric 3D Gaussian function [30]. Here, we described the detection probability per pixel of a molecule n by the detection point spread function PSF det in the focal plane: where the Gaussian width w det,D/A,0 of the PSF det, n, D/A depends on the detection system of the microscope and the emission wavelength of the fluorophores (Table C in S1 File). We defined the position of the individual SM, i.e., the center position (x n ,y n ), with subpixel resolution, and we implement longitudinal chromatic aberration z 0,D/A and linear focal drift along the optical axis with z 0 (t) = m Á t + z 0,D/A where m is the focal drift rate. A single fluorophore imaged with a magnifying microscope objective yields, according to Abbe's Law, a diffraction-limited spot of size w det,0 % λ em /(2 NA). Imaging single Cy3.5 molecules (λ em = 600 nm) with a water-immersion objective (60-fold magnification, NA = 1.2), results in a diffraction-limited spot size of 15 μm, a value that matches the camera-specific pixel size (Fig 3A top and Table C in S1 File). Upon increasing the magnification to a factor of 150, a typical value in the context of TIRFM, a diffraction limited spot spreads over a detection area of 3×3 pixel (Fig 3A bottom), or in the case of 2x2 pixel hardware binning over 1.5×1.5 pixel. The latter is used as default in our simulation software. It is important to mention that for a microscope-independent evaluation of SM methods, the Gaussian width of the PSF is given in pixel rather than micrometers.

Background simulation.
Background signal arises from various environmental and molecular sources, as well as photophysical processes [23]: (i) Environmental light (stray and ambient light) is usually suppressed via blackout-material (blinds, coating, masking tape) along the beam path and especially around the sample and/or the objective. (ii) The incident laser excitation light leads to Rayleigh or Raman scattering from the sample, the optics or even the detectors surface of the CCD camera [7]. The wavelength of scattered excitation light does usually not overlap with the spectral range of interest, i.e., the emission wavelength of the fluorophores, and can thus be suppressed by a suitable set of optical filters. (iii) Objective and sample chamber autofluorescence depends, for example, on the type and quality of glass. In Simulations of camera-based single-molecule fluorescence experiments addition, the sample itself can be autofluorescent when single-molecules are imaged in living cells or cell extracts [54] albeit the contribution of cellular autofluorescence can be reduced by choosing a dye pair that is excited at λ ex ! 520 nm [55]. Moreover, spurious adsorption of fluorescent impurities to the surface is observed on a regular basis. (iv) Studying intermolecular reactions by smFRET may require the presence of a fluorophore-labeled species at super-picomolar concentrations in solution [36,56]. Thus, freely diffusing and transiently adsorbing fluorescent particles contribute to the fluorescence that adds to the overall background signal. (v) The signal detected by a camera with closed shutter, so called dark images, numerically contributes to the background with a bias-offsetμ d . These dark images are not necessarily homogenous and the particular pattern is camera dependent [43]. Relevant contributions in a particular experiment, however, need to be identified on a case-by-case basis.
In our simulations, we modeled the background as spatially variable bg(x,y) D/A , which accounts for background sources (ii-iv). Here, the background is a function of the intensity of the light used for sample excitation. As laser beams with Gaussian cross-section profile are typically used for excitation [10], the excitation profile was described by a normalized, asymmetric 2D Gaussian function PSF ex , centered in the middle of the respective FOV [10,15]: where (x c ,y c ) is the center position of the FOV and w 0,ex, x/y is the width of the Gaussian background profile. Again, the defocusing function z 0 (t) accounts for lateral chromatic aberration and linear focal drift. It should be noted that this model neglects possible pattern of the dark images (v), as these systematic contribution to the measured image count signal are easily corrected for in an experiment [42,43]. To account for individual experimental conditions (optics, magnification, diameter, etc.) of the incident laser beam, w 0,ex,x/y can be adapted accordingly depending on the actual TIR profile of the microscope used (Fig 3B, related parameters in Schematic A in S1 File). To obtain a spatially invariant excitation profile and constant background, PSF ex must be set to 1.

Photon flux incident on the camera detector. The total detected intensity I(x,y,l) D/
A of all incident photons on the camera detector including the contribution of a background signal was defined as Iðx; y; lÞ D=A ¼ Iðx n ; y n ; lÞ D=A PSF det; D=A ðx; y; z ¼ 0Þ þ bg D=A Á PSF ex ðx; y; z ¼ 0Þ: The intensity per pixel in the respective donor and acceptor channel was calculated by the 2D integral over the corresponding pixel dimensions: 3.3.6 Camera noise model. Camera noise is composed of photon or photoelectron shot noise, amplification noise, spurious charge and camera-specific read-out noise, which contribute in an additive fashion [42,57,58]. Herein, we modeled EMCCD camera noise according to Hirsch et al. [43], taking only noise components into account, which are experimentally relevant: (i) We considered shot noise from all light sources contributing to photons reaching the detector pixel I (x,y,l) D/A . The probability of n ph incident photons to be detected on the detector pixel (x,y) during a frame l is given by a Poisson distribution P(n ph ;μ ph = I(x,y,l) D/A characterized by the mean number of incident photons μ ph . The simulation of sole photon noise is commonly used for the simulation of single photon detection applications [19,45,[59][60][61][62]. (ii) The second noise contribution is the camera shot noise σ pe common to all (EM)CCDs or sCMOS [42,43]. The probability of detecting n ph photons, i.e., producing n pe photoelectron, is the quantum efficiency η of the detector already included via the gamma correction, and follows a binomial distribution. The joint probability of the Poisson distribution and the binomial distribution is again a Poisson distribution P(n pe ;μ pe = ηI(x,y,l) D/A characterized by the mean number of accumulated photoelectrons. (iii) During the read-out process, electrons are shifted through the pixel of the detector area, the readout-register and the EM register by means of changing electrode voltages. This shift creates unwanted electrons, a so-called clock induced charge (CIC) as one component of spurious charge. The probability of observing CIC is very low, though, CIC nonetheless adds non-negligible intensity spikes to the detected signal. These may introduce artefacts in SMV analysis. (iv) Thermal noise, as another component of the spurious charge, is usually suppressed by cooling the detector to very low temperatures (-80˚C). Therefore, we did not consider it further. According to [43], the number of input electrons n ie in the EM register is a convolution of the two Poisson distributions of photoelectrons and spurious charge P(n ie ;μ ie = ηI(x,y,l) D/A + CIC) (v) Amplification noise of an EM register with gain g was modeled as a gamma distribution G(n oe ;μ oe = gμ ie ) where the mean number of output electrons μ oe = gμ ie consist of the Poisson distributed sum of all input electrons arriving in the EM register. (vi) Readout-noise σ d results from the conversion of the analog electron signal, i.e., the number of output electrons n oe , into a discrete image value μ ic = μ oe /s. We described readout noise with a Gaussian distribution N(n ic ;μ ic = μ ic,dark + μ oe /s, σ d ) characterized by a standard deviation σ d and the amplifier sensitivity or analog-to-digital factor s. We added a constant bias offset μ ic,dark which is usually added electronically to avoid negative image counts. (vii) Quantization noise σ q of the analog-to-digital converter was neglected, because moderate amplifier sensitivities are typically used. Thus, we yield the Poisson-Gamma-Normal (PGN) noise model of the EMCCD camera described earlier [43]  In difference to reference [43], we found that the CIC noise probability is well described by an exponential distribution Exp(CIC) = A CIC Á exp(−I/τ CIC ) that features the terms A CIC and τ CIC (Fig 4). Here, A CIC depends on the EM gain and has usually very small values whereas τ CIC varies with the hardware binning and the vertical clock speed to shift the electrons to the read-out register. Further, we described the number of output electrons n oe of the EM register with a normal distribution N(μ oe = I(x,y,l) D/A η D/A g, σ oe ) with mean image counts μ oe and variance σ oe 2 = μ oe. We approximated the PGN model by a In order to estimate the model parameters from experimental values, we recorded the signal of single surface-immobilized, Cy3 labelled RNA molecules with varying fluorescence intensity as described elsewhere [3,36]. The mean signal intensities μ ic and standard deviations σ ic are shown in Fig 4A. Following the nomenclature introduced in reference [43] and the standards Simulations of camera-based single-molecule fluorescence experiments for the characterization of image sensors and cameras [42], we defined the temporal variance σ ic 2 of the digital signal as the sum of the aforementioned noise sources: Here, σ d combines the noise related to the sensor read-out, the amplifier circuits and the spurious charge that is considered to be independent from the incoming photon signal. K = g/s is the overall system gain, i.e., a linear factor converting the charge into the digital output signal μ ic as follows: where μ d denotes the mean number of electrons and μ ic, dark = Kμ d the mean grey or offsetvalue of a camera accessed experimentally in the complete absence of light (Fig 4A and 4C).
The combination of Eqs (11) to (13) yields the variance of the image count signal of the camera assuming μ pe = σ 2 pe = ημ ph : The determination of K from regression analysis as depicted in Fig 4A is known as photon transfer method [57,58]. K was further used to calculate the mean number of incident photons μ ph as a function of the measured or simulated intensity signal of the camera μ ic and vice versa in Fig 4B. With Eqs (12) and (13) one can write the SNR ratio of the camera signal as follows: Fig 4B shows the SNR comparison of a real and an ideal sensor. An ideal sensor is characterized by a detection efficiency η = 1, negligible dark σ d = 0 and discretization noise σ q /K = 0. Here, the SNR of the investigated EMCCD as a function of photon counts is close to the square root function of an ideal sensor. We validated the introduced noise models Eqs (8) to (10) by comparing them to experimental noise distributions in absence and presence of light (Fig 4 and 4D). Except for the simple N noise model in Eq (10) in case of dark images, we found all noise models to be in good agreement with the experimentally derived signal intensity distribution.
Since recently, low-noise sCMOS cameras featuring megapixel resolution achieve similar or even higher frame rates as their low-resolution CCD counterparts [11]. This allows for a greater SM sensitivity and larger FOVs given the same magnification of the optical system. Even though low-noise sCMOS cameras are devoid of amplification noise, the SNR is comparable to EMCCDs. All relevant VSPs related to the camera noise model are implemented in the GUI shown in Schematic A in S1 File.

A generalized approach for comparing algorithms using simulated SMVs
Simulated SMVs may serve for comparing, evaluating and/or optimizing algorithms in smFRET data analysis, such as spot detection, background correction, time trace generation, and model selection. Herein, we used simulated SMV as ground truth providing the required information as video simulation parameters to calculate evaluation measures. These measures allow the direct comparison of algorithm performance. Thus, a systematic variation of relevant VSPs enables to find the confidence range of algorithms and their optimal parameterization. We optimized these parameters for a number of commonly used SM localization algorithms (Fig 5).

Evaluation measures and classification
To compare different algorithms used in the context of SMV data analysis, we first classified the true values known from the input VSPs as ground truth (GT) (Fig 5A, panel I). The output of the algorithm was then classified as true positive (TP, Fig 5A, panel II), i.e., correctly detected values, false positive (FP, Fig 5A, panel III), i.e., erroneously detected non-existing values, false negative (FN, Fig 5A, panel IV), i.e., undetected existing values which are incorrectly classified as non-value, or true negative (TN, omitted in Fig 5A as the herein tested algorithms do not classify TNs), i.e., values correctly classified as non-value. We then quantified the performance of different algorithms by computing the precision (also: positive predictive value) and the recall (also called sensitivity) [63,64]: In addition, the detection efficiency was quantified using the accuracy: Recall, precision and accuracy adopt values between 0 and 1 and reach 1 when TP 6 ¼ 0 and FP = FN = 0. Since most smFRET-specific algorithms do not produce a TN class output, TN = 0 in most cases. In the context of method evaluation, the accuracy can be used to optimize the input MP in order to assess the method-specific limits of its particular output parameters, to quantify maximum accuracies and to test computation times. Since it is generalizable and easily applicable to other algorithms, we believe that it may contribute to the standardization, and hence, the reproducibility of smFRET data analysis.

SMV test sets with different VSPs for method evaluations
We generated a number of SMVs which served as GT for the comparison of different algorithms. We classified the simulated SMVs depending on which VSP was varied and which SMV property was changed, respectively (Table D in S1 File). In particular, we varied the following VSPs: (i) The single molecule fluorescence intensity I tot,0 , which affects the SNR, (ii) the PSF width w det,0 , which mimics different imaging quantities, (iii) the molecule-to-molecule distances IMD and (iv) the molecule surface density ρ, both of which affect SM spot overlap, (v) the spatially variable background that affects the signal-to-background ratio, (vi) the trace length L and (vii) the ratio of forward/backward transition rate constants k ij / k ji , both of which affect the kinetic model determination and (vii) the heterogeneity associated with state-specific mean FRET values and the total emitted fluorescence intensity, both of which may influence the inferred model. All simulated SMVs are freely available for download at https://github.com/ RNA-FRETools/MASH-FRET for the individual assessment of smFRET-specific algorithms.

Application of simulated SMVs to optimize single molecule detection parameters
To achieve their optimal performance, SM localization algorithms require optimized input MPs (Fig 1). We optimized the input MPs of four spot detection algorithms, Houghpeaks and II) false negatives (FN, green) based on closest distance to a GT coordinate (greedy approach) and being located within the tolerance radius. Example III) and IV) compare a greedy and Gale-Shapely based approach where both show the same detection of SM but different classifications. For details, see text. (B) A subimage of 256×30 pixels from a 512×512 pixels video with a total of 24 × 12 single molecules. The total emission intensity (given in photon counts per frame) decreases from the left to the right. White circles mark the GT coordinates, blue, red and green mark TP, FP and FN, respectively. (C) Optimization of ISS algorithm parameterization for SM detection using SMV category (i) (Table D in S1 File) with w det,0 = 1 pixel. (C, left) Variation of model parameters (MP): 35 combinations of ISS input parameters I thresh (intensity threshold) and NhoodSize (spot size) and their obtained color-coded recall, precision, and accuracy are shown for molecules of 40 pc total emitted intensity I tot,0 . The respective range for a maximum (= 1) in recall, precision and accuracy are indicated by red bounding boxes. (C, middle and right) Variation of input parameters: Heat maps represent optimal parameterization of I thresh,opt (middle) and NhoodSize opt (right) to achieve maximum recall (1 st row), precision (2 nd row), and accuracy (3 rd row) as function of PSF size w det,0 and SM total emitted intensities I tot,0 . Note that I tot,0 is varied from 1-300 pc within a SMV, while w det,0 varies from video to video (category (i), Table D in S1 File). The corresponding example parameter optimization results for I tot,0 = 40 pc and w det,0 = 1 pixel on the left are highlighted by grey squares in all heatmaps. https://doi.org/10.1371/journal.pone.0195277.g005 Simulations of camera-based single-molecule fluorescence experiments (HP), in-series screening (ISS), Schmied (SCH) and TwoTone (2tone) ( Table 2), and used simulated SMVs of the category (i-ii) suitable for testing the ability of an algorithm to detect SMs given a heterogeneous distribution of the total emitted intensities I tot,0 . All four algorithms apply an intensity cut-off I thresh to the SM image prior to SM localization. This ensures that the algorithm can be used at different signal-to-background ratios. All tested detection algorithms process SMVs either sequentially (frame-by-frame) or upon averaging over all frames (frame/time-averaged SMV). The algorithms differ, however, with regard to the number of input parameter, the size parameter related to the PSF width (HP and ISS) or make use of a band pass kernel width (2tone). Different classification algorithms for TP, FP, and FN have been presented: (i) the Greedy approach which classifies TPs using the nearest neighbor criterion (type III in Fig 5A), and (ii) the Gale-Shapley approach, also known as stable marriage problem, which aims at maximizing TPs (type IV in Fig 5A) [63,65]. We used the Greedy approach, which is easy to implement and commonly applied. The maximum number of detected spots per frame was fixed to an upper bound larger than the GT to allow for the detection of false positives The resulting values for TP, FP, and FN were used to determine the algorithm-specific combination of input MPs yielding maximum precision, recall, and accuracy (Eqs 15 and 16). reach the maximum achievable value of 1, whereas the MP range for the maximum accuracy is more limited than for recall and precision. From this range, maximum accuracy is achieved exemplarily for I thresh = 13 pc and NhoodSize = 9 pixels (Fig 5C, solid grey square). I thresh is the more critical value for high recall values within the examined parameter space in category (ii) simulations co-varying I tot,0 and w det,0 . The NhoodSize parameter of ISS has no impact on the mean recall, precision and accuracy for w det,0 < 3 pixel and I tot,0 > 80pcpf. However, for conditions of decreasing SM signal (per pixel), i.e. towards w det,0 = 3 pixel and decreasing I tot,0 below 80 pcpf, maximal precision and accuracy can only be achieved by decreasing I thresh and increasing NhoodSize (Fig 5C, middle and right). Please note that there is a strong drop of the maximum accuracy upon decreasing I thresh (Fig 5C, left).
In order to optimize the input parameters of the other spot detection algorithms assessed, both I tot,0 and w det,0 , as well as the horizontal intermolecular distance IMD and w det,0 were covaried. The results are shown in Fig B in S1 File. For a constant fluorescence intensity, the parameter I thresh has a greater influence on accuracy than NhoodSize, in particular if IMD >> w det,0. . However, optimum values for I thresh and NhoodSize depend on the signal intensity I(x,y, l), the camera offset μ ic, dark , and the camera noise σ ic , and are thus specific to the present set of simulation parameters assessed. The parameter optimization for 2tone is special, because we found two parameter regimes yielding optimal results. The band pass kernel size BpSize of 1 and 3 requires I thresh parameters of 5.5 and 1.5, respectively. Such an anti-correlated effect between two parameters yielding maximum accuracy was demonstrated for ISS at increasing w det,0 and decreasing I tot,0 values too. Without going into methodical details, a modest increase of I thresh usually reduces FPs, which can be a consequence of image artifacts brighter than the local background caused by image processing steps used in a SM detection method. 2tone shows a clear preference of BpSize = 1 for the I tot,0 dependence which is not the case for the IMD dependence.
In summary, reducing only the parameter I thresh in the four investigated SM localization methods primarily results in more TPs and FPs, and thus, higher recall but lower precision. In general, a PSF theoretically leaks out of the ROI defined by the neighborhood size parameter NhoodSize, making e.g. ISS and HP prone for the detection of FPs at the boundary of the Nhoodsize-ROI. SM localization methods may contain image processing steps (usually implemented as black-box) that entail similar PSF smear effect (e.g. kernel size of 2tone). Thus, highest accuracy is often a compromise between I thesh and other parameters like Nhoodsize. In smFRET, for low molecule densities where PSF overlap is statistically rare, one might consider prioritizing precision over recall, a strategy to avoid duplicates (FPs located close to GT). min. distance to image edge 3 TwoTone (2tone) designed for analyzing smFRET SMVs Intensity threshold I thresh 1-6 a.u. [15] bandpass filter kernel BpSize (4-6) a.u. 1, 3, 5, 7, 9 Input parameters in brackets are specific to the analysis of time-averaged SMVs. a "The suppression neighborhood is the neighborhood around each peak that is set to zero after the peak is identified." Compare Matlab documentation.
https://doi.org/10.1371/journal.pone.0195277.t002 After parameter optimization, SM localization methods are compared in terms of recall, precision and accuracy using their individual set of optimized MP and varying SM total emitted intensity I tot,0 (Fig 6A). To enhance performance contrast, the individual result for recall, precision and accuracy were ranked between the four assessed methods with 1 to 4 points. In case of equal performance, ranking points are equally shared between the respective methods ( Fig  6B). Both methods 2tone and ISS have a low false detection rate and thus a high precision. This is is beneficial for the detection of SM spots and a subsequent SM trace generation to avoid "empty" time traces which contain only background signal. ISS achieved best values for accuracy, however, comparable with 2tone. In detail, we find ISS to perform better for higher and 2tone for lower intensity values I tot,0 (Fig 6B). Therefore, ISS and 2tone are both recommended for SM localization in SMVs akin to category (i). Note, that for other categories evaluation results may differ. To provide a deeper insight into the parameter optimization and method comparison presented in this section, we refer to https://github.com/RNA-FRETools/ MASH-FRET.

Conclusion
The first part of this article describes theoretical and practical aspects of simulating camerabased smFRET videos. Simulated SMVs presented herein feature true experimental conditions like realistic camera noise based on Emva Standard 1288 Release 3.1, the simulation of complex kinetic models of multi-state systems potentially showing kinetic heterogeneity, inhomogeneous Gaussian background, and variable SM surface densities. In addition, SM cross sample variability is accounted for in terms of Gaussian distributions of the total emission intensities, FRET states and gamma correction factors. The SMV simulation tool is integrated in our Matlab-based Multifunctional Analysis Software for Handling single molecule FRET data (MASH-FRET), which is freely available for download https://github.com/RNA-FRETools/MASH-FRET.
In the second part, we used the simulation tool to generate well-annotated test data sets that are independent of operating system and software. These test data sets are available online under https://github.com/RNA-FRETools/MASH-FRET and can be used as GT to evaluate computational methods in smFRET data analysis. We illustrated how simulated SMVs can be used to optimize the performance of four spot detection algorithms. For this purpose, we adapted the Greedy approach to categorize detected coordinates as TPs, FPs, and FNs, a classification that was used to define recall, precision and accuracy. This approach provides a standardized way of scoring the performance of spot detection algorithms. We provided a quantitative summary of optimized parameters of the methods assessed as a function of different PSF widths, SM emission intensities and intermolecular distances. We observed that method accuracies and parameterizations of all spot detection algorithms assessed are distinct functions of w det,0 , I tot,0 and IMD. Based on our results we want to emphasize the importance of a careful method parameterization prior to SM detection. Supporting information S1 File. Cumulative supporting information. Included: Abbreviations, Variables and Data Types/Section A; Simulation Software GUI/Section B and Schematic A; FRET theory and correction factors for smFRET type measurements [67,68]/ Section C and Tables A and B; Imaging properties for different optical imaging systems/Section D and Table C; Fitting algorithm and concepts/Section E; Camera noise model-analytical solution/Section F; SMV categories/ Section G and Table D; Export including metadata/Section H; Spectral representation of optical elements/Section I and Fig A; Method specific optimization of input parameters/Section I and Fig B. (PDF)