Time-resolved burst variance analysis

Quantifying biomolecular dynamics has become a major task of single-molecule fluorescence spectroscopy methods. In single-molecule Förster resonance energy transfer (smFRET), kinetic information is extracted from the stream of photons emitted by attached donor and acceptor fluorophores. Here, we describe a time-resolved version of burst variance analysis that can quantify kinetic rates at microsecond to millisecond timescales in smFRET experiments of diffusing molecules. Bursts are partitioned into segments with a fixed number of photons. The FRET variance is computed from these segments and compared with the variance expected from shot noise. By systematically varying the segment size, dynamics at different timescales can be captured. We provide a theoretical framework to extract kinetic rates from the decay of the FRET variance with increasing segment size. Compared to other methods such as filtered fluorescence correlation spectroscopy, recurrence analysis of single particles, and two-dimensional lifetime correlation spectroscopy, fewer photons are needed to obtain reliable timescale estimates, which reduces the required measurement time.


INTRODUCTION
The flexibility of proteins is key for their function. Resolving structural heterogeneity and quantifying the timescales at which proteins interconvert between different structural states has been a major goal in single-molecule fluorescence spectroscopy (1)(2)(3)(4). Single-molecule Förster resonance energy transfer (SmFRET) has particularly been used in the past 2 decades to study conformational changes in biomolecules (5,6). Most smFRET experiments use freely diffusing molecules. These experiments are easy to realize and avoid tethering of molecules to surfaces.
Naturally, a range of methods has been developed to extract dynamic information during the time molecules reside in the excitation volume of a confocal microscope ($1 ms). These methods range from dynamic photon distribution analysis (7), over maximum likelihood approaches (8)(9)(10)(11)(12) and equivalent Hidden-Markov model fitting such as H 2 MM (13,14) and multiparameter H 2 MM (15), fitting of FRET-histograms with different time binning (16), lifetime-filtered fluorescence correlation spectroscopy (fFCS) (17,18), two-dimensional lifetime correlation spectroscopy (19)(20)(21), recurrence analysis of single particles (RASP) (22,23), and lately a particularly promising approach using Bayesian nonparametrics (24)(25)(26). Each method has its merits and pitfalls. For instance, H 2 MM and maximum likelihood directly use the photon arrival times to optimize the parameters of a Submitted May 3, 2023, and accepted for publication July 3, 2023. WHY IT MATTERS Single-molecule fluorescence spectroscopy, particularly in combination with Förster resonance energy transfer, has been extremely successful in quantifying the dynamics of biomolecules. A toolbox of different methods is available to date that extracts dynamic information from the stream of photons emitted from donor and acceptor dyes. Yet, some of these methods require long integration times. In others, the presence or absence of dynamics is difficult to judge by eye and only fits with kinetic models provide this information. We therefore extended the popular method of burst variance analysis (BVA) to overcome some of these limitations. The new method termed timeresolved BVA quantifies dynamics from 5 ms to 5 ms at high accuracy with as little as 5000 bursts. Static and dynamic heterogeneity can be distinguished from each other, and even dynamics slower than the diffusion time can be quantified. Time-resolved BVA is a natural extension of classical BVA and therefore easy to implement by researchers in the field of single-molecule Förster resonance energy transfer. kinetic model and capture dynamics over a broad range of timescales. Dynamic photon distribution analysis computes FRET efficiency histograms by integrating the probability density that a molecule spends a certain time in each state of a kinetic model. The fit quality in these methods is often judged by generating FRET distributions from the model fit and comparing them to the experimental FRET histograms. Other methods such as lifetime-filtered fluorescence correlation spectroscopy, two-dimensional lifetime correlation spectroscopy, and RASP, first process the photon arrival times by computing correlation functions, frequency domain maps, or FRET histograms at different delay times. The preprocessed data are then used for model fitting. As an advantage, the presence of dynamics can already be inferred from the preprocessed data by eye, thus simplifying a model guess. On the other hand, these methods often require long measurements to obtain a high signal/ noise in the processed data.
Not standardly accounted for in these methods is static heterogeneity due to dye isomers or permutations of donor and acceptor positions. The latter is particularly prevalent in smFRET as donor and acceptor labeling is often done at cysteine residues, thus resulting in a mixture of labeling permutations. Burst variance analysis (BVA) (27) is a popular tool to identify both static and dynamic heterogeneity. Yet, BVA has mainly been used as a qualitative indicator for dynamics (3) as kinetic rates remain inaccessible. Here, we present an extension of BVA (27) termed time-resolved BVA (trBVA) that is also able to quantify kinetic rates from smFRET experiments of freely diffusing molecules between 200 ms -1 (5 ms) and 0.2 ms -1 (5 ms) with an error of a factor of 1.5. The method does not require long measurements and is easy to implement. To benchmark the robustness of trBVA, we performed smFRET simulations of dynamic particles and also applied the method to real single-molecule data of labeled DNA and protein.
We hope that trBVA will be a useful extension of the current smFRET analysis toolbox to identify biomolecular dynamics at timescales from micro-to milliseconds.

MATERIALS AND METHODS Theory
A photon burst i from a biomolecule labeled with donor (D) and acceptor (A) that diffuses through the confocal volume of a microscope contains d i donor and a i acceptor photons. The total number of detected photons in the burst is n i ¼ a i þ d i (including background photons), and the total number of bursts is N. We denote the uncorrected FRET efficiency as e and the corrected FRET efficiency as E (corrected for the differences in quantum yield of the dyes, cross talk between channels, background, and acceptor direct excitation; see section burst identification and data preprocessing). The idea of classical BVA is to partition photons of a burst into segments of m (typically m ¼ 5Þ consecutive photons. For each of these M i ¼ bn i =mc photon segments, the uncorrected FRET efficiency e ij (segment index j) is computed. Finally, we then calculate the variance of e using all segments of the N bursts The expected FRET variance of these segments in the absence of both dynamic and static heterogeneity (notably, Eq. 2 is also correct in the limit at which multiple states interconvert at timescales faster than the interphoton time), i.e., assuming the presence of only a single state, is due only to shot noise, and is given by The excess variance due to conformational heterogeneity is then given by the difference between Eqs. 1 and 2: Importantly, the analysis can also be performed with a subset of the N bursts. For instance, in a FRET-resolved trBVA version, the excess variance (Eq. 3) is computed for a set of bursts that lie within a chosen FRET efficiency range. If S 2 > 0, the FRET variance exceeds the shot noise expectation, thus indicating static or dynamic heterogeneity. The basic idea of trBVA is to vary the length m of the photon segments ( Fig. 1 A and B). Clearly, both variances s 2 and s 2 will change with m, but these changes will not be identical such that S 2 is itself a function of m. This function therefore contains information about the heterogeneity among and within bursts, which either is static or dynamic, i.e., time dependent. To extract this information, we derived an analytical expression for the excess variance of the subset of m-photon segments with specific time duration t, which we call the "t-specific excess variance" (Appendix I). Here, t is defined as the length of the time interval between the first and the last photon of a segment. Writing the FRET autocorrelation function as gðtÞ ¼ Cdeð0ÞdeðtÞD with deðtÞ ¼ eðtÞ À CeD, we obtain Importantly, for the ensemble of all m-photon segments, the time window t is a random variable with a conditional probability density function PðtjmÞ. Once PðtjmÞ is known, the excess variance due to conformational dynamics as function of m can be calculated from the following: The change of S 2 with increasing m can therefore be computed by knowing the autocorrelation function Cdeð0ÞdeðtÞD and the distribution PðtjmÞ. The autocorrelation function can be easily computed for any kinetic model. If K is the rate matrix of the model, p eq is the population vector of conformational states at equilibrium (Kp eq ¼ 0), and e is a diagonal matrix with the same dimensions as K whose diagonal elements are the FRET efficiencies of each conformational state, then the FRET autocorrelation function can be expressed as (8) where 1 is a vector of ones. For instance, a model, in which two states with FRET efficiencies e 1 and e 2 interconvert with rates k 12 and k 21 , has the correlation function A fit of S 2 with Eqs. 4, 5, and 7 would provide the two unknown quantities Cde 2 D and k obs ¼ k 12 þ k 21 if PðtjmÞ was known. In fact, this distribution can be extracted from the experimental data FIGURE 1 Scheme of the trBVA procedure and simulated FRET efficiency histograms (corrected) for a freely diffusing dynamic particle. (A) In trBVA, the photons from acceptor (red) and donor (green) in a burst i are partitioned into segments of length m. For each segment, apparent (uncorrected) FRET values e ij are computed, and the variance of these FRET values is studied as function of m. (B) Illustration of the trBVA variance analysis. The FRET efficiencies of individual segments within a burst i are depicted as function of time for two scenarios: a hypothetical burst without dynamics, i.e., only including shot noise (left column), and a "measured" burst with a conformational transition (right column). The FRET efficiencies are shown for three values of segment lengths m (indicated). The variances of FRET efficiencies are depicted as gray shaded areas. The trBVA excess variance is the difference between the measured variance (right) and the shot noise variance (left). Importantly, trBVA excess variance for a given segment length m is computed from the segments of all bursts, not only for a single burst as shown in B. (C) Brownian dynamics simulation of FRET efficiency histograms (corrected) for a particle diffusing freely through a confocal spot including bleaching of donor and acceptor. The particles switched between two states with corrected FRET efficiencies E 1 ¼ 0.9 and E 2 ¼ 0.1. The kinetic forward (k 12 ) and backward (k 21 ) rates were assumed to be identical. The FRET efficiency histograms are shown for different values of k 12 and k 21 (bottom). directly. We first determine the time duration of all photon segments of length m for all bursts or a subset of bursts within a chosen FRET window E l % E < E l þ DE in the FRET-resolved version. A histogram of these times Hðt i jmÞ for equally spaced time bins t i with i ¼ f1; 2; 3; .; Kg then provides a reasonable estimate for PðtjmÞ. For data fitting, we therefore use Eq. 5 in discrete form: Hðt i jmÞ: (8) For completeness, we also provide the explicit forms of Ds 2 ðm; tÞ for a two-state and a three-state system in Appendix II. For comparison, we also computed the donor-acceptor cross correlation function G DA ðtÞ ¼ Cn D ðt 0 Þn A ðt 0 þtÞD =Cn D DCn A D for the selected bursts. Here, n D ðt 0 Þ and n A ðt 0 Þ are the photon counts at time t 0 . To extract the relaxation time, G DA ðtÞ was fitted with the empirical function Here, k obs ¼ k 12 þ k 21 is the observed rate of conformational changes, t D is an empirical timescale to describe the decay of G DA ðtÞ due to diffusion, and b is a stretching exponent.

Data simulation
To test the accuracy of trBVA in extracting kinetic rates from smFRET experiments, we simulated photon time traces of diffusing particles that switch between two conformational states (1 and 2) described by kinetic rate coefficients k 12 and k 21 . The FRET efficiencies of the two states were E 1 ¼ 0:1 and E 2 ¼ 0:9, respectively. The diffusion of the particle through the confocal volume was modeled via Brownian dynamics simulations with the software package Fretica (https://schuler.bioc.uzh.ch/programs/), developed by Daniel Nettels and Benjamin Schuler (University of Zurich). The Stokes radius of the particles was set to 4.3 nm, which corresponds to a medium-sized protein, and the particles diffused in a solvent with the viscosity of water at 25 C, i.e., 1 mPas, resulting in a diffusion coefficient of 5 x 10 À 5 mm 2 =ms. The simulation was initialized by randomly placing particles in a simulation sphere with a radius of R ¼ 3 mm. The number of initial particles was drawn from a Poisson distribution with a mean n 0 ¼ 4 3 pRc 0 with a bulk particle concentration of c 0 ¼ 50 pM. The simulation was performed in spherical coordinates assuming for simplicity radial symmetry of the confocal volume, which is located at the origin. Brownian motion is simulated using the following: Here, rðtÞ is the radial distance at simulation steps t ¼ 1.T, where T is the length of the simulation in steps of Dt ¼ 1 ms, i.e., the time between two simulation steps, D is the diffusion coef-ficient, and Dr is a random distance drawn from a normal distribution with zero mean and a variance s 2 Dr ¼ 2DDt. Each particle is simulated until it leaves the simulation sphere. To ensure a constant mean concentration of particles near the center of the sphere, the particle loss at the sphere's surface is compensated by periodically (periodicity T new ) placing new particles inside the sphere near the boundary. The distribution of new particles c new ðrÞ that entered the sphere after time T new is obtained by solving the radial diffusion equation with the initial condition cðr < R; t ¼ 0Þ ¼ 0 and the boundary conditions cðr ¼ R; tÞ ¼ c 0 and cðr /0; tÞ ¼ 0. The solution is known (28) and given by cðr; tÞ The mean number of new particles entering the sphere is then computed by integrating over the volume of the sphere: After each time interval T new , a random number of new particles was drawn from the Poisson distribution with mean n new . The particles were placed at radial distances randomly chosen from the distribution with the density function P new ðrÞ ¼ 4pr 2 c new ðrÞ=n new for r < R. In total, we simulated particle trajectories for 1800 s. Once the particle trajectories were simulated, we added conformational dynamics simulated according to the rate equation dp dt ¼ Kp; where p is the population vector of four states: low FRET (DA 1 ) with FRET efficiency E 1 , high FRET (DA 2 ) with FRET efficiency E 2 , donoronly (D), and acceptor-only (A) in the basis fD; DA 1 ;DA 2 ;Ag. The rate matrix K is a combination of the rate matrix K 0 for conformational transitions between DA 1 and DA 2 and the rate matrix K bl describing photophysical effects, photobleaching in our case, with the bleaching rates k a and k d for acceptor and donor fluorophores located at the origin ðr ¼ 0), respectively. We assumed a bleaching timescale of k a ¼ k d ¼ 5 Â 10 À 4 ms À 1 for the simulations. The position-dependent profile IðrÞ that accounts for the illumination intensity at different positions in the confocal volume is given by For each particle with the diffusion trajectory rðtÞ and starting time t 0 , a random state trajectory s t is simulated according to Eqs. 14-18 with the program Fretica. The initial state sðt 0 Þ was chosen randomly according to the initial probabilities for the four states given by the vector p 0 with the same basis as p. We chose an equal distribution of high-and low-FRET species and the same number of donor-only and acceptor-only molecules with p 0 ¼ 0:1; 0:8 k21 k12þk21 ; 0:8 k12 k12þk21 ; 0:1 T . In addition, we set the total photon rate at the center of the excitation volume to l tot ¼ 0:4 ms À 1 and introduced realistic background photon rates of l d ¼ 5:6 10 À 3 ms À 1 for the donor channel and l a ¼ 3 10 À 3 ms À 1 for the acceptor channel. To model the experimental situation in a realistic fashion, we also introduced different detection efficiencies for the where Q a;d and h a;d are the quantum yields and detection efficiencies for acceptor and donor dye, cross talk (leakage) between donor photons in the acceptor channel (b ¼ 0:054), and the probability to directly excite the acceptor with the donor excitation laser (a ¼ 0:048). As we introduced donor-only and acceptor-only molecules together with the possibility of photobleaching, we also simulated pulsed-interleaved excitation (PIE) of both dyes with g PIE ¼ 2 (6,29). To this end, experimental instrumental response functions were used to generate the photon distributions after donor and acceptor excitation within one PIE period. Finally, a time-tagged time-resolved file containing the simulated photons was generated. Simulations of a three-state model were performed in the same manner.

Burst identification and data preprocessing
After simulating photon traces based on the kinetic model described above, the time-tagged time-resolved file was processed with standard single-molecule analysis tools (6) for generating corrected FRET efficiency histograms. Importantly, for the calculation of variances for BVA, raw photon counts, without correction, were used to calculate apparent FRET efficiencies, also known as proximity ratios. Unless stated otherwise, the photon trajectory was binned into time windows of 100 ms. A burst is defined as a collection of consecutive bins with more than two photons per bin and a total photon number of at least 100 photons after donor excitation. The corrections included background, differences in the brightness of donor and acceptor, channel cross talk, and acceptor direct excitation. The procedure is described in detail elsewhere (6,30). The corrected photon numbers of donor (n DD ) and acceptor (n DA ) after donor excitation were used to compute the FRET efficiency of the burst via To exclusively identify molecules that contain both dyes, we computed the stoichiometry ratio for each burst via Only bursts with S PIE < 0:65 were retained for further analysis. Since bursts were identified based on photon counts after donor excitation, molecules without donor were automatically excluded from the analysis. To also exclude bursts in which the acceptor bleached during the transit of the particle through the confocal volume, we further selected bursts in which the mean detection time of photons was similar after donor and acceptor excitation. We define where Ct Dex D and Ct Aex D are the mean detection times (in ms) after donor and acceptor excitation, respectively. Including shot noise, the asymmetry value a PIE has a standard deviation given by where the prime indicates the uncorrected photon counts. We chose a restrictive threshold of s PIE < 0:15 to exclude bursts with bleached acceptors.

Global trBVA
To test the ability of trBVA ( Fig. 1 A) in quantifying timescales of conformational dynamics, we simulated the photon emission process for freely diffusing molecules in a photon-by-photon manner. We modeled molecules that switch between two conformational states 1 and 2 with "forward" rate k 12 and "backward" rate k 21 .
The corrected FRET efficiencies (E) of the two states were E 1 ¼ 0:1 and E 2 ¼ 0:9, which corresponds to the uncorrected values e 1 z0:2 and e 2 z0:9. For simplicity, we assumed identical rates in both directions. At a slow exchange rate of k 12 ¼ k 21 ¼ 0:1 ms À 1 , i.e., one transition per ten milliseconds on average, the FRET efficiency histogram shows two well-separated peaks with shot noise limited width at the expected FRET efficiencies (Fig. 1 C). Intermediate values between the dominant peaks become prominent with increasing exchange rates, as more molecules change their conformation while diffusing through the confocal volume. At higher rates, the FRET peaks start to coalesce, and at the highest exchange rate of k 12 ¼ k 21 ¼ 50 ms À 1 , the FRET peaks merged completely, thus giving the impression of a single conformational state. To analyze these data with trBVA, we computed the variance of FRET fluctuations by partitioning bursts into consecutive segments with m photons (Fig. 1 A). As outlined in the theory section, computing the variance of these segments and subtracting the shot noise contribution one would have if there was a single state with a FRET efficiency equal to the population weighted mean of the states, provides the excess variance S 2 (Eqs. 1-3). Fig. 2 A demonstrates that S 2 first increases and then decreases with increasing size of the photon segments m. The trBVA traces obtained from the data (Fig. 2 A) can now be used to determine the apparent relaxation time t ¼ ðk 12 þ k 21 Þ À 1 of the conformational fluctuations using the experimentally determined distribution Hðt i jmÞ of the time duration of m-photon segments. Examples are shown in Fig. 2 B. For m ¼ 2, the distribution is a decaying function as expected based on photon counting theory (8).
For higher values of m, Hðt i jmÞ shows a clear maximum due to the fact that a successive emission of several photons causes a delay between the first and the m th photon that leads to the rise at short times. To fit the trBVA traces, we use Eqs. 4, 7, and 8 to compute S 2 for each value of m and minimize the least The fit contains two parameters, the amplitude of the FRET correlation function Cde 2 D and the kinetic rate k obs ¼ k 12 þ k 21 (Eq. 7), i.e., the eigenvalue of the rate matrix. The fits provide an excellent description of the simulated data over a broad range of exchange rates (Fig. 2 A). An alternative method to determine kinetic rates would be to compute the FRET autocorrelation function directly or analogously, the donor-acceptor cross correlation G DA ðtÞ ¼ Cn A ðtÞn D ðt þtÞD=Cn A DCn D D for the data (Fig. 2  C). Distance dynamics lead to a rise of the cross correlation amplitude since donor and acceptor signal are anticorrelated. Yet, the finite burst duration causes an additional decay in G DA ðtÞ at the timescale at which molecules diffuse through the confocal spot. This diffusion amplitude dominates G DA ðtÞ, and slow dynamics at timescales close to the diffusion are difficult to identify (Fig. 2 C). This problem is circumvented with trBVA. A comparison of the apparent relaxation times t ¼ ðk 12 þ k 21 Þ À 1 from the trBVA analysis with the true values used in the simulation demonstrates an excellent agreement (Fig. 2 D). For dynamics across three orders of magnitude (5 ms to 5 ms), trBVA provides estimates of t with less than twofold deviation from the ground truth (Fig. 2 E). Even dynamics slower than the diffusion of molecules through the confocal spot can be obtained. The reason for this surprising result is that S 2 is bounded by two limits. For dynamics much faster than the experimental interphoton time, the lower boundary is given by S 2 ðmÞ ¼ 0 (Appendix III). Yet, for extremely slow dynamics, the FRET autocorrelation function is approximately constant (gðtÞzCde 2 D) but different from zero. Under this condition, the excess variance is given by S 2 ¼ Cde 2 Dð1 À m À 1 Þ, which is an increasing function of m and represents the upper boundary (Fig. 2 A, top). Notably, this increase is not in conflict with the central limit theorem. The total variance indeed decreases with increasing m (Appendix III). Instead, the increase of S 2 results from an inaccurate estimate of shot noise (Eq. 2) in the presence of static heterogeneity (Appendix III). Importantly, even slight deviations from the ð1 À m À 1 Þ dependence require a finite decay time in gðtÞ, which explains the success of trBVA at slow timescales. Notably, this is a helpful feature to identify static heterogeneity. For instance, labeling proteins with donor and acceptor is often done via two cysteine residues, which unavoidably results in two labeling permutations. If the molecular brightness of the dyes differs in the two variants, they will exhibit different FRET efficiencies, and S 2 will follow the ð1 À m À 1 Þ dependence. In comparison to trBVA, the relaxation times from the donor-acceptor cross correlation function G DA ðtÞ are highly inaccurate at the diffusion timescales (Fig. 2 D). Compared to the 1.5fold error in trBVA, the cross correlation analysis deviates from the ground truth sevenfold at a relaxation time of 1 ms.

FRET-resolved trBVA
Similar to regular BVA, also the time-resolved version can be used to investigate dynamics in different regions of the FRET efficiency histogram. In FRETresolved trBVA, the segments of bursts within a partic-ular FRET range are analyzed. Importantly, selecting bursts within a FRET range means selecting trajectories according to their mean FRET efficiency. In a two-state system, S 2 for bursts with FRET values different from the ensemble average will therefore be biased. Bursts with FRET values substantially lower than the ensemble average will contain trajectories with longer dwell times in the low-FRET state and shorter dwell times in the high-FRET state (Fig. 3 A).
The opposite happens when bursts with substantially higher FRET than the ensemble average are being selected. As the observed rate is a sum of forward and backward rate, the faster rate, i.e., the shorter dwell time, dominates. Hence, at the flanks of the FRET efficiency distribution, the observed exchange rates will in general be higher than the correct value ( Fig. 3 B, top). with an arbitrarily chosen uncorrected FRET value e are given by p 0 1 ¼ ðe 2 À eÞ =ðe 2 À e 1 Þ and p 0 2 ¼ ðe À e 1 Þ =ðe 2 À e 1 Þ, respectively. The primes indicate that these occupancies differ from those of the whole ensemble of molecules. The amplitude of the FRET autocorrelation at this FRET value is Cde 2 D ¼ p 0 1 p 0 2 ðe 2 À e 1 Þ 2 (see also Eq. 7), which can be re-written as Hence, the amplitude follows a second-order polynomial in e where the roots identify the position of the states (Fig. 3 B). Notably, this relationship is independent of the true relative populations of the two states (p 1 and p 2 ). The amplitude analysis is therefore suited to identify the (uncorrected) FRET values of the interconverting states e 1 and e 2 . In general, FRETresolved trBVA experiments can be used to identify the positions of FRET states. However, kinetic rates should always be inferred from S 2 using all bursts and not from FRET-resolved trBVA! This is important as the FRET-dependent rates will always exhibit a minimum at a FRET value centered between e 2 and e 1 , i.e., the point at which p 0 1 ¼ p 0 2 , irrespective of the abundance of both conformers in the whole ensemble. Moreover, the observed rate at the minimum is higher than the eigenvalue of the system (Fig. 3 B, top) because trajectories without transitions (bursts with e 1 and e 2 ) are underrepresented in this FRET range. To exemplify this deviation, we simulated a more complicated system in which three states with different FRET efficiencies (E 1 ¼ 0:1, E 2 ¼ 0:5, E 3 ¼ 0:9) interconvert at different timescales (Fig. 3  C). We assume that state 1 and 2 exchange at a slow timescale with the rates k 12 ¼ k 21 ¼ 1 ms À 1 , whereas state 2 and 3 exchange an order of magnitude faster with k 23 ¼ k 32 ¼ 10 ms À 1 . A comparison with the case in which exchange is hundredfold slower than the diffusion time through the detection volume shows how drastically dynamics can alter the appearance of FRET efficiency distributions (Fig. 3 C). In the presence of fast exchange at two different timescales, the FRET efficiency histogram shows a major peak at an apparent FRET efficiency value of 0.7, a minor peak at 0.1, and a floor of events in between the peaks. In a quantitative global analysis, we first computed S 2 for all bursts. As expected, the trBVA trace increases and decreases with m (Fig. 3 D). A fit with a singleexponential FRET correlation function (Eqs. 4, 7, and 8) already provides a reasonable fit (Fig. 3 D, top). Yet, the residuals clearly show discrepancies between data and fit. Indeed, a fit with a double-exponential correlation function, which corresponds to the correct three-state model (Appendix II), provides an excellent description of the data (Fig. 3 D, bottom) and gives the correct eigenvalues (Fig. 3 E). To exemplify how static heterogeneity would manifest in trBVA, we set the fitted rates in the correlation functions to zero (Fig. 3 D). The comparison shows that dynamics lower the amplitude of the trBVA trace and introduces the decay at large m. In the more qualitative FRETresolved rate analysis, we calculated trBVA traces for bursts with different FRET efficiency values. An empirical fit with a single-exponential FRET correlation function provides apparent exchange rates for the individual FRET efficiency values. These rates exhibit a nontrivial FRET dependence (Fig. 3 E). A minimum is observed at FRET values between state 1 and 2. Starting from the minimum, the exchange rates increase toward lower FRET values as expected (compare to Fig. 3 B, top). However, although the rates also increase toward higher FRET values, a flattening of this dependence between state 2 and 3 is found. The position coincides with the position of the major peak at high FRET, which can be taken as indication that molecules in this peak dynamically switch at a fast timescale. Yet, the analysis is qualitative as the rates at both minimum and flattening point are substantially higher than the eigenvalues (Fig. 3 E).
As a rule of thumb, steep changes in exchange rates along the FRET coordinate indicate regions with biased trajectories and therefore regions close to the positions of the FRET states. FRET-independent exchange rates (minima or flat regions in the rate profile) indicate trajectories with strong exchange between states. Yet, care has to be taken as 1) flattening of the rate profile might not always be clearly visible, and 2) states in exchange rarely have identical populations such that exchange rates should not be inferred from the rate-FRET profile but always from the trBVA decay of the whole ensemble.

Probing the dynamics of double-stranded DNA
As an application of trBVA, we probed the dynamics of double-stranded DNA (dsDNA) breathing. Structural fluctuations in dsDNA have previously been measured using fluorescence quenching (31). A relaxation time of $50 ms was found for these local opening-closing motions, a timescale well within the regime that can be probed with trBVA. We performed smFRET experiments on dsDNA at neutral and acidic pH. At acidic pH, dsDNA is known to be destabilized (32) due to the protonation of DNA bases, and we expect a significant difference in the amplitude and/or timescales of these motions between pH 7 and pH 4. We generated 12 dsDNA samples of 84 bp length each that were derived from a naturally occurring promoter sequence in Bacillus subtilis (33). The samples were site-specifically labeled with AlexaFluor488 as donor and AlexaFluor594 as acceptor at varying positions, thus spanning the full FRET efficiency range from low to high values. We performed short 5-to 10-minlong experiments using PIE (29) and identified bursts as described in the materials and methods section. As expected, the FRET efficiency histograms of these samples span the full FRET range (Fig. 4 A). Notably, the widths of the FRET efficiency histograms are significantly increased at pH 4 compared to pH 7, suggesting that the drop in pH either alters the timescales of distance dynamics or the amplitude or both (Fig. 4  A). We then used trBVA to analyze FRET fluctuations in these samples. A comparison of S 2 at m ¼ 5 shows a substantially increased fluctuation amplitude at pH 4 compared to pH 7 (Fig. 4 A). This variance is reduced at m ¼ 46, suggesting a pronounced microsecond decay. An overview of the decays indeed demonstrates the presence of structural fluctuations that are intensified at low pH (Fig. 4 B). A fit with a twostate model provides an empirical description of these data with relaxation times that are in rough accord with the previous estimate of 50 ms at dsDNA samples with intermediate FRET efficiencies, whereas substantially larger relaxation times were found for samples with extremely low and high FRET values. However, the fits do not properly capture the trBVA decays. To obtain a better description of the traces, we also fitted with double-exponential FRET autocorrelation functions, which is equivalent to a model with three states. This model describes all trBVA traces well and results in two relaxation times (Fig. 4 C). The fast relaxation time is closest to the previous estimate of 50 ms at samples with low FRET efficiencies (Fig. 4 C, bottom). Yet, for samples with high FRET efficiency, the fast relaxation time drops to values in the order of 2-10 ms. This very fast timescale could be caused by transitions of the dyes into photophysical triplet states or by direct contacts between donor and acceptor that lead to quenched dye complexes (Fig. 4 C). However, we also identify a slow timescale in the order of 500-2000 ms, which apparently represents slower motional modes of the structure of the DNA. In fact, previous results demonstrated that dsDNA breathing motions exhibit nonexponential dynamics (31) such that our three-state model only provides a simplified description of the true dynamics. As a second example, we determined the folding-unfolding dynamics of the B-domain of protein A (BDPA) from Staphylococcus aureus, a protein that had previously been used to benchmark RASP (22). The particular variant used here (F13W/Y14C/G29A/P57C) has a folding relaxation time of 0.93 ms -1 at 2.5 M of the denaturant guanidinium chloride (GdmCl) at 37 C. The protein was labeled at positions 14 and 57 using AlexaFluor488 and AlexaFluor594, and other details of the experiments (buffer, laser intensity, etc.) can be found in Hoffmann et al. (22). Since the experiment was not performed with PIE, we selected the bursts for trBVA based on their FRET value to exclude molecules with inactive acceptor (Fig. 5 A, inset). The trBVA trace cannot be described with a single-exponential FRET correlation function (Fig. 5 A), and a double-exponential function was required. Whereas the fast rate (l 1 ¼ 411 ms -1 or 2.4 ms) is associated with the smaller amplitude (36%) and is well in the regime of dye triplet blinking, the slower rate (l 2 ¼ 0.9 ms -1 or 1.1 ms) dominates the amplitude and indeed corresponds to the timescale observed with RASP (1.4 ms -1 ) and temperature jump experiments (0.93 ms -1 ).
In summary, the relaxation times of DNA breathing and of the folding and unfolding of BDPA obtained with trBVA agree well with previous measurements. Compared to our simulations, a very fast relaxation component at timescales of a few microseconds is found in both data sets and might reflect the triplet blinking of our dyes.

CONCLUSION
We presented a time-resolved version of BVA termed trBVA and developed a theoretical framework to apply trBVA in a quantitative manner to smFRET experiments of diffusing molecules. TrBVA is capable of identifying dynamics in biomolecules at timescales from 5 ms up to 5 ms with remarkable accuracy. Using simulated data, we also showed that trBVA can be used in a FRET-resolved manner to identify the FRET values of states that are in exchange. In more complicated cases in which more than two states exchange, FRET-resolved trBVA merely provides qualitative information about the FRET efficiency values of the states. In general, FRET-resolved trBVA is a qualitative tool to understand the complexity of the dynamics at hand.
Finally, we demonstrated the ability of trBVA to identify dynamics in real experiments using the examples of the breathing motions in dsDNA and of fast folding/ unfolding kinetics of a protein. We are therefore convinced that trBVA is an excellent addition to the existing toolset of smFRET.

APPENDIX I
To arrive at Eq. 4, we start from three common assumptions: 1) The total photon rate does not fluctuate in time.
2) The probability, eðxðtÞÞ, of observing an acceptor photon is determined by the spatial distance xðtÞ between donor and acceptor dyes. 3) Dye excitation-emission cycles are fast compared to the interphoton time.

A B
For an m-photon segment with a given set of photon arrival times ft i g i ¼ 1.m resulting from a single trajectory xðtÞ, the probabilities for the individual photons to be detected in the acceptor channel are given by fe i g i ¼ 1.m ; where e i ¼ eðxðt i ÞÞ. Our goal is to first compute the first and second moment of the distribution of e for a single trajectory and then to average them over all trajectories. The probability to observe a acceptor photons is therefore given by the Poisson binomial distribution (34) P pb ða fe i g m i ¼ 1 Þ, which generalizes the ordinary binomial distribution in that the probabilities for individual trials do not need to be equal. The mean of the distribution is known to be CaD ¼ P m i ¼ 1 e i ; and the variance is For the given set of arrival times ft i g; we then get the mean and variance of e as For the second moment, we can then write the following: We now need to average over all sets of arrival times. To this end, consider photon segments that have a fixed duration t between the first and last photon, whereas all m À 2 photons in between the first and last photon have random arrival times. For the moment, we only consider a specific trajectory xðt 0 Þ of a molecule, but we will average over all trajectories at a later stage. To compute the total probability of obtaining a acceptor photons in a photon segment, we need to average over all possible arrival times ft i g i ¼ 1.m . Note that the arrival times are not an ordered set. Whereas the arrival times are fixed for the first and last photon (t 1 ¼ 0 and t m ¼ t), the arrival times of the remaining m À 2 photons are independently and uniformly distributed random values between 0 and t. The probability density function of such sets ft i g i ¼ 1.m is given by As mentioned above, the number of acceptor photons for a given set ft i g obeys Poisson binomial statistics. The mean FRET efficiency for the trajectory xðt 0 Þ is then given by CeDðm; t; xðt 0 ÞÞ ¼ Similarly, with the use of Eq. I.3, the second moment is dt 0 dt 00 ½eðxðt 0 ÞÞeðxðt 00 ÞÞ ' (I. 6) At this stage, we average both results I.5 and I.6 over the ensemble of distance trajectories and denote this average as C.D x . In this notation, the FRET efficiency averaged over all trajectories CeðxðtÞÞD x and deðxðtÞÞ are then written as CeD x hCeðxðtÞÞD x ; deðxðtÞÞ After averaging, the result for the mean FRET efficiency is Similarly, averaging Eq. I. 6 gives The variance of segments with length m and time duration t is then (I. 10) where the first term is the variance of a binomial distribution. The second term is the t-specific excess variance that contains information about the conformational fluctuations and with gðtÞ ¼ Cdeð0ÞdeðtÞD x is given by ðt À t 0 Þgðt 0 Þdt 0 ! ; (I.11) which is identical to Eq. 4 in the main text. It is noteworthy that the above result is similar to the result obtained by Gopich and Szabo for the case of fixed time bins instead of variable segment lengths (16). The main difference is that all photons in a fixed time bin are randomly distributed, whereas in BVA, the arrival times of the first and the last photon in a segment are fixed. When relaxing this constraint, the probability density function in Eq. I.4 becomes and solving the integrals in I.5 and I.6, the t-specific excess variance becomes

APPENDIX II
Here we provide the formulas for the t-specific excess variance Ds 2 ðm; tÞ used for data fitting with the two-state and three-state model. As outlined in Eq. 7, the correlation function of a two-state model is given by g 2 À State ðtÞ ¼ a e À lt ; where a ¼ k12k21 k12þk21 ðe 1 À e 2 Þ 2 and l ¼ k 12 þ k 21 : Similarly, the correlation function for a three-state model can be obtained from Eq. 6 using the appropriate rate matrix K and the diagonal matrix containing the FRET efficiencies of the three states e. For instance, for the model shown in Fig. 3 C (top), we have The rate matrix K has two nonzero eigenvalues l 1 and l 2 , and the correlation function can be written as g 3 À State ðtÞ ¼ a 1 e À l1t þ a 2 e À l2t . Inserting these expressions in Eq. 4 and solving the integrals gives

APPENDIX III
For dynamics that are slow compared to the duration of a segment t, we can approximate the correlation function by a constant number Cdeð0Þdeðt 0 ÞD x zCde 2 D x for t 0˛ð 0; tÞ. In this case, the t-specific excess variance (I.11) simplifies to x is the variance of the FRET efficiency over all conformational states. Note that Eq. III.1 is an increasing function with m. Following Eq. I.10, the total t-specific variance in the static case is then Ceð1 À eÞD x m : From the last expression, it is clear that the total t-specific variance is a decaying function of m (third term in III.2), as expected based on the central limit theorem. Using the expression CeD x ð1 À CeD x Þ=m to estimate shot noise using binomial statistics, which is only correct in the absence of heterogeneity, causes the factor ð1 À m À 1 Þ in the t-specific excess variance (III.1). A better estimate for shot noise in the static case is therefore Ceð1 À eÞD x = m.
The other limit is given when the timescale of conformational dynamics (t) is much faster than the interphoton time. Hence, the time duration of any segment (t) is much longer than t or, equivalently, PðtjmÞ ( 1 for t $ t. When evaluating the integral in Eq. 5, the distribution PðtjmÞ has nonzero weights only for those values of Ds 2 ðm; tÞ for which gðtÞ $ 0; and therefore Ds 2 ðm; tÞ $ 0. Correspondingly, Eq. 5 evaluates to S 2 $ 0, which is the lower boundary of S 2 .

APPENDIX IV
To judge the accuracy of the parameters that can be obtained from trBVA traces, we provide a lower limit for the amplitude of the autocorrelation function. In a two-state system, the amplitude of the FRET autocorrelation function will increase with De 2 ¼ ðe 2 À e 1 Þ 2 (see Eq. 7, Fig. S1 A). Clearly, the higher the FRET-separation of the two states is, the higher is the amplitude of the autocorrelation function. A lower limit of De should be given by shot noise. For simplicity, we assume that the two states exchange with slow dynamics compared to the diffusion time through the confocal volume. In addition, we assume that their shot noise variance is identical and given by s. To distinguish the two states in a FRET efficiency histogram, the separation between the states should exceed the combined shot noise variance of both states. We therefore require De > 2s: (IV.1) From this expression, we can get a lower limit for the amplitude of the FRET autocorrelation function. The amplitude is given a ¼ k 12 k 21 De 2 =l 2 ; where l ¼ k 12 þ k 21 . Given a value for the difference De, the best possible case (the highest amplitude) would be at k 12 ¼ k 21 . Any other combination of k 12 and k 21 would result in even smaller amplitudes. With k 12 ¼ k 21 , the amplitude reduces to a ¼ De 2 =4; and with Eq. IV.1, we would require a > s 2 : (IV.2) Importantly, even if s of a single state cannot be reliably determined because of a high overlap between the states or because of fast exchange dynamics, the above inequality would still provide a lower limit for the amplitude below which a determination of the FRET autocorrelation function becomes unreliable. For instance, for a measured FRET distribution centered at e ¼ 0:5 and a threshold of 100 photons, one would ideally like to have an amplitude of a > 0:5 2 =100 ¼ 2:5 10 À 3 and De > 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 0:5 2 =100 p ¼ 0:1. We tested this estimate with simulations assuming identical kinetic rates (k 12 ¼ k 21 ¼ 1 ms À 1 ) but different values of De (Fig. S1). The apparent kinetic rate l is obtained with good accuracy down to a value of De ¼ 0:2. However at Dez0:1 and at an amplitude of az3:5 10 À 3 , i.e., close to our accuracy estimates of De > 0:1 and a > 2:5 10 À 3 , the fitted rate exceeds the ground truth threefold ( 2 Â a 1 e À l 1 t þ a 2 e À l 2 t Ã þ 4ðm À 2Þ " a 1 ð1 À e À l 1 t Þ l 1 t þ a 2 À 1 À e À l 2 t Á l 2 t # þ 2ðm À 2Þðm À 3Þ " a 1 ðe À l 1 t À 1 þ l 1 tÞ l 2 1 t 2 þ a 1 ðe À l 2 t À 1 þ l 2 tÞ l 2 2 t 2

#)
: (II.4) ( Fig. S1 B), indicating that the parameters of the FRET autocorrelation function cannot be reliably determined. In addition, we tested the sensitivity of trBVA to the photon threshold used to identify bursts. Simulations show that the photon threshold has no impact on the determined kinetic rates from trBVA as long as the FRET efficiency separation between the states fulfills Eq. IV.1 (Fig. S1 B, inset). Yet, at an extremely low separation Dez 0:1; a doubling of the photon threshold from T ¼ 100 to T ¼ 200 indeed lowers the discrepancy between fitted rate and ground truth.