Unraveling multi-state molecular dynamics in single-molecule FRET experiments -- Part II: Quantitative analysis of multi-state kinetic networks

Single-molecule F\"orster Resonance Energy Transfer (smFRET) is ideally suited to resolve the dynamics of biomolecules. A significant challenge to date is capturing and quantifying the exchange between multiple conformational states, mainly when these dynamics occur on the sub-millisecond timescale. Many methods for the quantitative analysis are challenged if more than two states are involved, and the appropriate choice of the number of states in the kinetic network is difficult. An additional complication arises if dynamically active molecules coexist with solely static molecules in similar conformational states. To address these problems, we developed an integrative analysis framework that combines the information from FRET-lines, time-correlated single photon counting, and fluorescence correlation spectroscopy. While individually, these methodologies provide ambiguous results for quantitatively characterize the dynamics in complex kinetic networks, the global analysis enables accurate determination of the number of states, their kinetic connectivity, and the transition rate constants. To challenge the potential of smFRET to study multi-state kinetic networks, we benchmark the framework using synthetic data of three-state systems with different kinetic connectivity and interconversion rates.


Analysis of multi-state kinetic networks
Biomolecular dynamics are often complex, involving multiple conformational states and sub-states that interconvert over a wide range of timescales from nanoseconds to minutes and hours. Single-molecule FRET (smFRET) experiments provide a wealth of information about the molecular system and are ideal for resolving these dynamics. In multiparameter fluorescence detection (MFD) experiments on freely diffusing molecules, this information is encoded in the time-ordered sequence of the detected photons recorded with picosecond resolution. Various analysis methods have been developed over the years to obtain quantitative information on the structural dynamics of the biomolecular systems. The most widely used methods are the statistical analysis of FRET-efficiency histograms (Photon Distribution Analysis, PDA) 1-4 , intensity-based fluorescence correlation spectroscopy (FCS) [5][6][7][8] , and time-resolved fluorescence decay analysis (TCSPC) 9 , but other approaches have been applied as well [10][11][12] . Each representation of the experimental data and the corresponding analysis method has its strengths and weaknesses to determine the fluorescence properties of the species, detect their kinetic connectivity, and quantify the rate constants (Table 1). Established methods such as FCS and TCSPC are computationally fast. They rely on established algorithms to find the optimal parameters of a physical or empirical model that describe the experimental data. While TCSPC is ideally suited to resolve the FRET efficiencies of the contributing states, conformational dynamics from nano-to milliseconds can be resolved by fluorescence correlation spectroscopy (FCS), fluorescence lifetime correlation spectroscopy (FLCS), and filtered-FCS, which utilize statistical weighting to recover species-specific correlation curves 6,13,14 . These correlation approaches work well for homogenous samples of dynamic molecules interconverting between two conformational states. However, they are challenged by the increased complexity of many biological systems that involve three or more states interconverting on different kinetic timescales or contain heterogeneous mixtures of static and dynamic molecules. Therefore, to unravel the complex dynamics of such systems, a holistic approach combining multiple methods is required. Existing methods for the quantitative analysis of dynamics are applied to a reduced representation of the single-photon-counting data. At the same time, the full potential of the multidimensional dataset is not utilized. This multidimensional information is revealed in the pairwise histograms of averaged fluorescence observables -although the informational content is likewise reduced due to the averaging performed over each single molecule event 15 . One example of how the multidimensional information can be utilized is the pairwise plot of the intensity-based FRET efficiency E and the intensity-weighted average donor fluorescence lifetime ⟨ D(A) ⟩ F . The correlation between these two FRET indicators enables the detection of conformational dynamics. Further, using FRET-lines to describe the exchange between different states, the connectivity within the kinetic network is revealed as described in Part I 3 . Ideally, a multi-state kinetic model would be directly fit this multidimensional data set. However, to our knowledge, a quantitative description of the complete multidimensional histogram is currently limited to computationally expensive stochastic simulations. The stochastic noise of Monte-Carlo simulations also makes this approach difficult to apply in optimization routines, which converge more rapidly if analytical expressions are employed. While such expressions are known for simple cases 16 , they are currently unavailable for the multi-state networks discussed in this work. Here, we take a step towards a holistic analysis framework by combining FRET-lines, TCSPC, and FCS in a global approach to quantify the exchange in multi-state kinetic networks. In a first step, the correct kinetic model is identified by a graphical analysis using FRET-lines, defining the number of states and their linkage. The exchange rates are then quantified using a global analysis of the donor fluorescence decay and the color correlation functions. The framework is applied to simulated data sets of multi-state systems with a binary exchange between two states in the presence of a background of static molecules. When only the TCSPC and FCS information is used, ambiguous solutions are obtained that differ in the kinetic connectivity of the species and the fraction of molecules participating in the dynamic exchange. To resolve this ambiguity, FRET-lines provide a graphical analysis of the kinetic connectivity of states and permits the estimation of the equilibrium constant from the peak of the dynamic population in binary systems. For systems involving a fast dynamic exchange between more than two states, additional information is required. Using simulations of three-state systems, we illustrate the potential of filtered-FCS to detect the direct exchange between different species in complex networks and deduce the kinetic linkage, even in this challenging case. Lastly, we derive relations between the correlation amplitudes and the single-molecule FRET indicators E and ⟨ D(A) ⟩ F , highlighting the connections between the different representations of the data and the future possibility to extend this holistic approach to data analysis.

Nomenclature for multi-state systems in smFRET
Purely static or dynamic biomolecules are rarely found in nature. It is often observed that biomolecules can be activated through allosteric effects such as binding of cofactors or regulators, post-translational modifications, or conformational changes in associated domains, switching the molecule from a static into a dynamic state [17][18][19] . In such situations, molecules with the same FRET efficiency may either be static or participate in the conformational dynamics, introducing a degeneracy into the analysis in which the same observed FRET species may belong to different states that are either static or dynamic.
In FRET experiments on freely diffusing single molecules, the accessible timescale of dynamics is limited by the diffusion time to >0.1 ms -1 causing additional complications because transitions between conformational states on slower timescales could appear as pseudo-static populations in the analysis. To avoid confusion about the physical description of the biomolecular system as static or dynamic and to classify observed populations in the experiment, we propose a concise nomenclature for smFRET experiments in Figure 1. A conformational state C (i) is defined as a distinct structural state of the biomolecule that can be classified as static or be in dynamic exchange with other conformational states. The alternation between dynamic and static states of the biomolecule may be subject to allosteric regulation, biomolecular interactions, or covalent modifications. In the smFRET experiment, conformational states are observed indirectly through the fluorescence properties of the covalently linked dyes, such as the FRET efficiency, fluorescence lifetime, or fluorescence anisotropy of the donor and acceptor fluorophores. An observed fluorescence species o (i) is generally assigned to a single conformational state. However, due to quenching or sticking of the fluorophores, different fluorescence species may belong to the same conformational state. On the other hand, multiple conformational states may belong to the same fluorescence species if the fluorescence properties do not change significantly (see Supplementary Note 1 for an overview of potential ambiguities). In smFRET experiments, fluorescence species are observed as populations p (i) in the one-or two-dimensional histograms. Due to dynamic averaging during the diffusion time, a population may originate from a mixture of different fluorescence species. Dynamic and static populations may be distinguished in a plot of the FRET efficiency E against the donor fluorescence lifetime ⟨ D(A) ⟩ F (Figure 1, bottom). The static populations, originating from the fluorescence species o (1) and o (2) , lie on the static FRET-line. In contrast, the dynamic population p (1,2) shows the characteristic dynamic shift (ds) from the static FRET-line, as introduced in Part I. The heterogeneity within the dynamic populations can be resolved by a subensemble analysis of the fluorescence decays. The assignment of static and dynamic populations is complicated when the sample contains a mixture of static and dynamic conformational states of identical FRET efficiencies (Figure 1, right). In the limit of fast dynamic exchange, a dynamic population is shifted from the static FRET-line and separated from static populations. As the conformational exchange becomes slower and approaches the diffusion time of the molecule, there is a probability that dynamic molecules do not undergo a conformational change during the observation time. While originating from dynamic conformational states, these singlemolecule events will show as a pseudo-static population on the static FRET-line and are difficult to separate from actual static populations. As shown below, the fraction of dynamic molecules is a central parameter in the quantitative analysis of such heterogenous multi-state systems by correlation methods. In the next section, we will first show how a graphical analysis can be employed to estimate the equilibrium constant of the dynamic exchange in the background of static states. Table 1: Comparison of methods to analyze conformational dynamics of multi-state kinetic networks in smFRET. The methods are compared to their ability to identify the different states (fluorescence properties and fractions), kinetic connectivity, dynamic exchange, and the accessible timescales. MLE: Gopich-Szabo photon trajectory analysis using maximum likelihood estimation 20 , BVA: burst variance analysis 21 , FRET-2CDE: FRET twochannel kernel density estimator 22 . * Timescale of dynamics 500 µs. ** The lower limit depends on the average inter-photon time. Faster timescales are accessible for higher signal count rates.

Method
Identification of states

Graphical analysis of dynamic populations
Quantitative information on conformational dynamics is encoded in the shape of the FRET efficiency histogram, as in dynamic photon distribution analysis (PDA) 3,23,24 . These analyses, however, are challenged if the experiment contains a mixture of static and dynamic molecules due to the difficulty of distinguishing actual static and pseudo-static molecules. Pseudo-static molecules are dynamic molecules that, by chance, remained in one conformational state during the transit through the observation volume. This section will describe how the separation of static and dynamic molecules in the E-〈 〉 histogram can provide quantitative information on the equilibrium constant by a graphical analysis of the peak of the dynamic population. In the description of FRET-lines, the timescales of the dynamics are not considered explicitly. For a dynamic system, the distribution of the state occupancies depends on the microscopic interconversion rates and the observation time 3,16,25 . For the calculation of dynamic FRET lines in Part I, we have instead considered all possible values for the state occupancy, ∈ 0,1 . In other words, we have replaced the true distribution of the state occupancies by a uniform distribution with equal probability for all values of ( const.). FRET lines may, however, still be used to address the timescale of dynamics qualitatively. In the absence of dynamics, the two-dimensional histogram will reveal distinct static populations as limiting states, which fall onto the static FRET line. In the case of fast exchange between distinct FRET states, the conformational dynamics are averaged for every singlemolecule event, resulting in a single population representing the equilibrium. "Fast" interconversion relates to the timescale of diffusion (~1-5 ms), and generally classifies processes on a timescale of 100 µs and below. This single peak will be shifted from the static-FRET line in the two-dimensional histogram, as described before. The slow transition between limiting states, such as the dynamics on the timescale of diffusion or slower, leads to a broadening of the observed distributions; thus, the shape of the distribution depends on the timescale of the dynamics. (E-H) for different timescales of conformational dynamics from slow (k 12 = k 21 = 0.01 ms -1 ) to fast exchange (k 12 = k 21 = 10 ms -1 ). The population's shape and location indicate the dynamic timescale, while the dynamic FRET-line is independent of the rate constants and describes all timescales. The diffusion time in all simulations was 1.5 ms.

8
To illustrate this effect, we performed simulations of a two-state system with FRET efficiencies of 0.2 and 0.8 ( Figure 2). We set the backward and forward rates equal ( ) and varied them from 0.01 ms -1 to 10 ms -1 , at a constant diffusion time = 1.5 ms. When the rate constants are significantly slower than the inverse diffusion time, 1/ = 0.67 ms -1 , the two subpopulations are separated because molecules do rarely interconvert during the observation time (Figure 2 A,E). With increasing rate constants, the molecules are more likely to change their state during the observation time, resulting in single-molecule events with intermediate FRET efficiencies (Figure 2 B-C, F-G), while complete averaging is observed at fast interconversion rates (Figure 2 D,H). However, the dynamic FRET line describes all possible mixing ratios between the involved states, regardless of the magnitude of the rate constants, and applies to all cases. In the first part of the paper, we defined the observed deviation of a population perpendicular to the static FRET-line as the dynamic shift ( Figure 3A). Interestingly, the dynamic shift of the population for the simulated system with does not reach its maximum possible value even for fast dynamics ( = 10 ms -1 , Figure 3A and Figure 2D). In the limiting case of fast dynamics, the observed dynamic shift of the dynamically-averaged population depends on the FRET efficiencies ( and ) and the species fractions ( and ) of the two states (see Supplementary Note 2): where is the species fractions of state 1. The dependence of the observed dynamic shift on the average FRET efficiency of the dynamic population, , is illustrated in Figure 3B. The maximum dynamic shift, ds , is given by (see Part I): which is obtained at a species fraction of: The dynamic FRET-line for a two-state system with FRET efficiencies = 0.8 (red) and = 0.2 (blue) reaches a maximum dynamic shift, ds max (cyan), as indicated by the arrow, at a FRET efficiency of = 0.6, corresponding to a species fraction of = 2/3. A smaller dynamic shift, ds obs (magenta), is observed for the simulation with rate constants of k 12 = k 21 = 10 ms -1 = 0.5, compare Figure 2D). The static FRET-line is given in black, and the dynamic FRET-line is color-coded according to the species fraction . B) The dependence of the dynamic shift on the average FRET efficiency for the two-state system shown in panel A. 9 corresponding to a value of 2/3 in the given example. The maximum dynamic shift defined by eq. (2) thus serves as a figure of merit to judge whether dynamics could be detected in the experiment in the ideal case, while the observed value of the dynamic shift also depends on the equilibrium constant according to eq. (1). Beyond the qualitative assessment of dynamics by visual inspection of the two-dimensional histogram, it is possible to extract quantitative information about the dynamic equilibrium from the twodimensional histograms. The equilibrium species fractions and (i.e., the fraction of molecules in state 1 or 2) can be determined from the average FRET efficiency for all single-molecule events, 〈 〉 , by where and are the FRET efficiencies of the limiting states. From the species fractions, the equilibrium constant can be calculated, which relates to the kinetic rates by: .
While for purely dynamic systems 〈 〉 is readily calculated from the experimental dataset, the presence of additional static states will result in incorrect values for the species fractions. In the case of a mixture of dynamic and static molecules, it would be advantageous if the equilibrium fraction could be obtained from the position of the dynamic population alone, which is most easily defined by its peak value or mode. For fast interconversion, the mode of the population in the two-dimensional histogram directly corresponds to the average FRET efficiency 〈 〉 . If the timescale of kinetics becomes comparable to the diffusion time, the mode of the dynamic distribution deviates from the actual value of the species fraction . To study this effect, we consider the distribution of FRET efficiencies for a dynamic system explicitly. The average FRET efficiency within a single-molecule event, depends on the fraction of time spent in the different states, equivalent to the species fractions : 1 Here, is the species fraction of state 1 in a single-molecule event, which is different from the equilibrium species fraction discussed before. If we know the state occupancy distribution of , , we can calculate the distribution of FRET efficiencies . In general, takes a complex mathematical form (a complete derivation is given in Supplementary Note 3), but can be simplified to the sum of three terms: Here, the first and the second terms describe the probability that the molecule was in state 1 or 2 at the beginning of the observation time and did not interconvert (pseudo-static states). Also, the term describes the dynamic part of the distribution of , such as in the case that molecules switched between the states at least once during its observation time. This term is given by 3, 25-27 : where 1 and and are the modified Bessel functions of the first kind of order zero and one, respectively. For fast interconversion (or long observation times), the pseudo-static terms vanish while only the dynamic term remains. On the other hand, for slow dynamics, the static terms tend to the equilibrium fractions of the two states, while the dynamic term vanishes. Theoretical FRET efficiency distributions for a two-state dynamic system are given in Figure 4A for rate constants of = 1 ms -1 and = 3 ms -1 at different observation times . For short observation times ( = 0.25 ms), only the pseudo-static peaks remain, while complete averaging is observed for long observation times ( = 100 ms). The average FRET efficiency 〈 〉 of 0.65 relates to the equilibrium species fractions by eq. (4). Considering the case that additional static species contribute to the average FRET efficiency, we would like to infer the equilibrium species fraction only from the dynamic part of the distribution ( Figure 4B). The property of that is most easily inferred from the two-dimensional histograms is its maximum or peak (compare Figure 2). For an observation time = 100 ms, the modal value corresponds to the equilibrium fraction of state 1, 0.75 ( Figure 4B, top). However, as the observation time decreases, deviates from the equilibrium fraction to the point where the modal value coincides with the pseudo-static population ( = 1, Figure 4B bottom).
Fortunately, the relationship between the actual equilibrium fraction of state 1, , and the modal value of the distribution, , is approximately linear, enabling a simple conversion between the two quantities by ( Figure 4C and Supplementary Note 3): where , is the limiting equilibrium fraction at 0 (i.e., the ordinate intercept in Figure 4C) that depends only on the average number of transitions during the observation time given by : where .
The model value of the species fraction, , can be obtained from the modal value of the FRET efficiency distribution, , obtained by graphical analysis if the FRET efficiencies of the limiting states are known: (11) If the sum of rates is known from FCS analysis or other methods, the equilibrium fraction and the equilibrium constant can be determined from equations (9) and (10), which allows determining the microscopic rate constants quantitatively using the equilibrium information obtained from the graphical analysis of the two-dimensional histograms. While the approach in principle requires observation time windows of identical length , it may, as an approximation, be set to the diffusion time for datasets of single-molecule events of freely diffusing molecules.

Analytical description of FCS curves
Fluorescence Correlations Spectroscopy (FCS) relies on the fluctuations of recorded signals to characterize molecular interactions such as binding and unbinding, chemical kinetics, and diffusion of fluorescent molecules 5,8,28 . Also, when combined with FRET, FCS, the conformational dynamics can be quantitatively determined [29][30][31][32][33] . Typically, the fluorescence signals are collected over specific spectral detection windows. Here, we refer to the correlation analysis of the fluorescence intensities of a donor and an acceptor fluorophore (monitored in two detection channels commonly named "green" and "red") as color-FCS. We avoid the conventional abbreviation FRET-FCS to differentiate it from the related method of filtered-FCS, which relies on FRET to distinguish different species but does not explicitly use color channels. Analytical models for color-FCS are usually limited to kinetic networks involving two states due to the increased number of parameters of multi-state systems and the limited experimental information available 32 . Advanced correlation methods take advantage of the lifetime information available with pulsed laser excitation, allowing one to interrogate biomolecular dynamics by twodimensional maps of fluorescence decays 34,35 or filtered correlation algorithms (fFCS) 6,13,36,37 . By using an additional dimension of the collected data, these methods offer the potential to interrogate more complex kinetic networks. Quantitative analysis of FCS experiments requires a set of model-specific analytic functions that describe the time-evolution of the correlations. Often, smFRET and FCS experiments are used to study unimolecular reactions, wherein a biomolecule switches between different conformational states during the observation time. These dynamic molecules may be found together with molecules that are stable on the timescale of seconds to minutes and are thus considered as static in the single-molecule experiment.
In our model, we thus consider the coexistence of static and dynamic molecules with identical properties of the respective conformational states (Figure 1). Dynamic molecules may change between different conformational states during the observation time, resulting in variations of the fluorescence properties such as the FRET efficiency and donor fluorescence lifetime, while these properties remain constant for static molecules. For unimolecular reactions, the time evolution of the different states is described by a system of linear differential equations, which for a system with three dynamic states is given by: The rate constants describe the rates of transition from state j to state . In general, this is expressed in matrix notation as: , (12) where is the transition rate matrix and is the vector of the total fractions of the species. In the following, we denote the fraction of static molecules by , and the fraction of dynamic molecules by . Both fractions are normalized to one, i.e. ∑ 1 and ∑ 1. Hence, the total fraction of a species is given by the sum of the fractions and , weighted by the total fraction of dynamic molecules as: where and , are the vectors of the fractions of the dynamic and the static states, and describes the fraction of molecules that participate in dynamic exchange.
The correlation function, , is modeled based on the set of reaction rate constants, fluorescence properties, and the population of the static states. The general definition of the correlation function between two time-dependent signals and is given by: where 〈… 〉 denotes the time average over a long measurement. In the following, we assume that the signal fluctuations due to the diffusion of molecules and the conformational dynamics arise from independent sources (i.e., that they are statistically independent), which is generally fulfilled for experimental systems. The contributions of diffusion and dynamics can thus be treated separately, and the correlation function is given by the product: where is the correlation time, is the average number of molecules in the observation volume, the factor , describes the kinetic exchange and the diffusion of the molecules. For a 3D Gaussian detection profile, the factor is given by: where is the diffusion time. The parameters and are the width of the focal and the axial plane of the detection volume, respectively, where the intensity decays to 1/ of the maximum value.
The kinetic factor , depends on the correlation matrix that describes the time evolution of the fluctuations of the species populations, and the signals and of the observed signal of the different species (given as column vectors). Then, the kinetic part of correlation can be expressed in matrix notation as: , (17) where ⋅ ∑ is the average of signal over the species fractions , which corresponds to the time average of the signals under the assumption that the system is ergodic. In color-FCS, the vectors and correspond to the green and red signal intensities and . In filtered-FCS, and are the fractional fluorescence intensities of the species, obtained by weighting the signal based on the fluorescence decay using the filter functions for each species a and b, respectively. The general solution for the correlation matrix of the kinetic network in the presence of static states is given by: can be solved using the eigenvalue decomposition (EVD) of the transition rate matrix: where are the eigen-matrices and λ the eigenvalues of which relate to the measured FCS relaxation times by: Equation (18) can then be expanded into a sum of exponential terms and substituted into eq. (17) to obtain the following general expression for the kinetic correlation function in the presence of static states: , 1 ; ; , 1 … 1 ; where the matrix elements are given by: (22) and the factors describe the contrast between the species and and are given by:

14
For details on the derivations, see Supplementary Note 4.
To obtain the correlation functions, we need to define the correlated signal vectors S a and S b . In color-FCS, these signals are the detected "green" (donor) and "red" (acceptor) signal intensities S G and S R : where is a vector whose elements correspond to the FRET efficiencies of the fluorescence species, is the molecular brightness of the donor in the absence of FRET, is the crosstalk from the donor fluorophore into the red detection channel. is a combined correction parameter relating the donor and acceptor fluorescence quantum yield and the detection efficiencies of the green and red channels [38][39][40] .
For simplicity, we assume that the crosstalk of the donor fluorescence into the red detection channel of the acceptor is zero, that the -factor is one, and consider no background noise. In the expressions of the normalized correlation function, the scaling factor cancels out. For the simplest case of a twostate dynamic system in the presence of static states, we then obtain the general expression for the kinetic correlation function: where the pre-factor ∂ depends on the channels that are correlated: where the average FRET efficiency is given by the average over the total species fractions . The complete derivation of the analytical form of the correlation function for two-and three-state systems is outlined in Supplementary Note 4. Corresponding expressions had previously been obtained for twostate dynamic systems in the presence of a third static state 31 .

Ambiguities in color-FCS
Before applying the formalism derived in the previous section for the quantitative analysis of the simulated datasets, we emphasize why the combination of FCS and TCSPC is needed. Color-FCS (or FRET-FCS) is generally underdetermined as there are more model parameters than experimentally accessible parameters. Thus, it is required that the FRET efficiency of at least one of the two states is known, but better results are obtained if both FRET efficiencies are restraint 6 . The origin of this ambiguity is outlined in the following. For a purely dynamic two-state system, i.e., in the absence of static molecules, the expression for the correlation function given in equation (25) simplifies to: In an experiment, we measure two autocorrelation functions (GG, RR) and two cross-correlation functions (GR, RG). The cross-correlation functions contain identical information, and the time constant of the exponential term is shared between all correlation functions. Thus, we determine three correlation amplitudes, ∂ 1 , and the decay rate of the exponential term, . As is evident from eq. (26), however, the cross-correlation amplitude relates to the auto-correlation amplitudes by: and contains no independent information. The system is thus underdetermined, as we have access to only three experimental observables compared to the four parameters of the model ( , , and

15
. The ambiguity between the model parameters takes a complex form and is illustrated in Supplementary Note 5. This ambiguity is resolved if the FRET efficiencies of the states are known from single-molecule FRET efficiency histograms or fluorescence decay analysis. In the following, we explore the combination of FCS with TCSPC to restrain the FRET efficiencies of the states using the information provided by the fluorescence decays, which enables quantitative analysis of the kinetics by FCS.

Joint analysis of fluorescence decays and FCS
To unambiguously resolve all contributing states and the interconversion, we combine the information provided by FCS and TCSPC and optimize all model parameters globally. While FCS is sensitive to the relaxation rate constants, TCSPC informs the FRET efficiencies and the total species fractions. Thus, the two methods provide orthogonal information that defines the FRET efficiencies of the states and the transition rate matrix. The global analysis is also expected to stabilize the optimization algorithm and reduce the uncertainty of the model parameters.
The donor fluorescence decay of the FRET sample | depends on the FRET efficiencies of the species and the associated total species fractions . We assume that the time scale of fluorescence and the time scale of dynamics are decoupled. In other words, the fluorescence lifetime is much shorter than the relaxation time of the kinetic processes. Therefore, the fluorescence decay of the ensemble of molecules can be described by the total fractions of the FRET species and their FRET efficiencies : To optimize the model parameters, we define a global goodness-of-fit function, , as the sum of the squared weighted deviations for TCSPC, , and FCS, : Here, the weighting factors, and , account for the non-uniform noise in the data. For TCSPC, the weighting factor is estimated based on the experimental counts under the assumption of Poissonian counting statistics as

Analysis of three-state kinetic networks
To benchmark the global analysis framework for the analysis of multi-state systems, we simulated a series of experiments. We consider a heterogeneous mixture of various static and dynamic FRET species ( Figure 5, Tables S1-S2). No linker dynamics were included in the simulations. First, we consider four distinct static species with low, medium, and high FRET efficiency (LF, MF, MF' and HF species, respectively, Figure 5A and S9). As expected, the four populations in the two-dimensional histograms lie on the static FRET-line, and no indication for conformational dynamics is seen. We then simulated heterogeneous mixtures of three FRET states, two of which are in dynamic exchange ( Figure 5 B-C). The kinetic rates were chosen such that the resulting dynamic population has an average FRET efficiency that is identical to the MF' population in the static mixture, resulting in almost indistinguishable one-dimensional FRET efficiency histograms for the three simulations. By overlaying the dynamic FRET-lines connecting the static species, the interconverting states of the dynamic population can be assigned (solid lines). Dynamic FRET-lines of species that are not in dynamic exchange do not intersect with the dynamic population (dashed lines). Moreover, as the dynamic populations are positioned directly on the limiting binary dynamic FRET-lines, we can exclude the possibility of ternary exchange between all three states, which would instead result in a population positioned within the area defined by the three limiting lines (compare section 3.5 of Part I). The graphical analysis by FRET lines thus provides a simple approach to determine the kinetic connectivity of the network. For the moment representation ( Figure 5, bottom), dynamic FRET-lines can be drawn as simple lines. While the equilibrium constant may be extracted from the plots given the fast dynamics in these examples (as described in section 3), additional information is required to quantify the interconversion rates of the kinetic network. In the following, we apply the global analysis of FCS and TCSPC to the simulated datasets shown in Figure 5. For the simulation of four static species ( Figure 5A), the absence of conformational dynamics  Table S1 and S2. is confirmed by the absence of a kinetic contribution in the FCS curves (SI Figure G2). For the simulations with dynamics between two species (Figure 5 B-C), static and dynamic species are indistinguishable in analyzing the fluorescence decays, which are well described by a three-state model ( Figure 6 A,C). The green/red cross-correlation curves show a pronounced anti-correlation, which reveals conformational dynamics (Figure 6 B,D). In the analysis of the correlation curves, a single relaxation time of ~100 µs is sufficient to describe the data in both cases. This implies that only two of the three species are in dynamic exchange, consistent with the graphical analysis performed by FRET-lines. We thus consider three possible kinetic schemes that differ in the assignment of the purely static state: In this notation, MF, HF ⇌ LF refers to the scheme where the HF and the LF species are exchanging, and the MF species is purely static. To distinguish amongst the possible kinetic schemes in eq. (31) and to quantify the microscopic parameters, we use the global analysis of FCS and TCSPC. In total, there are nine model parameters: three FRET efficiencies for the FRET species HF, MF, and LF; two exchange rate constants to describe the exchange among the two dynamic states, two independent fractions of the static states (the third is determined by the other two), and the probability that a molecule is in the dynamic state . These microscopic parameters define the experimental observables such as the slope of the fluorescence decays or the amplitudes and relaxation timescales of the FCS curves. The fluorescence decays are fully described by the FRET efficiencies and total species fractions (see eq. (29)), while the relation between the microscopic parameters and the amplitudes of the FCS curves is more complex. In the description of the FCS model function, we had split the contributions of the FRET efficiencies of the different states from the quantities that depend only on the parameters of the kinetic network (see eq. (25)). The relaxation time of the kinetic amplitude is given by the inverse of the sum of the interconversion rates, . The amplitudes of the auto-and cross-correlation curves depend mainly on the total species fractions and FRET efficiencies, which are determined from the information provided by TCSPC. As described in detail in Supplementary Note 6, the only new information obtained from the amplitudes of the FCS curves is the relative amplitude of the kinetic term (the pre-exponential factor in eq. (94)), given by . Thus, only seven parameters are available from the experiment: the FRET efficiencies and static fractions obtained from TCSPC, and the sum of the rates and the product from the FCS curves. The system is hence inherently underdetermined, and ambiguity is expected between the fraction of dynamic molecules and the interconversion rates that define the dynamic fractions :

Resolving complex kinetic networks using the global analysis framework
To test this prediction, we sampled the probability distribution of the parameters for the possible kinetic networks given in eq. (31) using a Markov chain Monte Carlo (MCMC) approach. The resulting distributions are shown as two-dimensional contour plots in Figure 7A. Indeed, the experimentally accessible parameters show defined, narrow distributions due to the high signal-to-noise ratio of the simulated data. However, the distributions of the microscopic parameters exhibit the expected ambiguity for the total fraction of dynamic molecules . Due to the ambiguity of the experimental parameters, two different kinetic schemes (blue: LF⇌MF, red: LF⇌HF) are consistent with the data as indicated by the identical (right). Each of these two kinetic schemes is additionally consistent with a range of model parameters. If the correct kinetic scheme (blue: LF⇌MF) and the total fraction of dynamic molecules are known, the unique solution can be identified (gray dashed line). C) Ambiguities arising in the global analysis. The first ambiguity arises from the fact that the system is inherently underdetermined, preventing both and to be resolved at the same time (dashed lines). Additionally, due to the ambiguity of the FCS amplitudes, different kinetic schemes can be compatible with the data, resulting in two branches of the solution space in the given case (red: false solution, blue: true solution). D) A graphical analysis by FRET-lines resolves the connectivity within the kinetic network. The position of the dynamic population on the dynamic FRET-line between the states LF and MF shows that the kinetics occur between the LF and MF states, while the HF state is static (blue line in C). The dynamic FRET-line for the competing solution (red line C) is shown as a red dashed line. To resolve the ambiguity between the parameters and in the general case of , it is additionally required to determine the dynamic species fractions . In the case of fast dynamics as given here, this parameter can be estimated from the position of the dynamic population along the dynamic FRET-lines. The quantities a, b and c represent the length of the vectors connecting the LF state and the dynamic population (orange), the MF state and the dynamic population (green), and the LF and MF states (blue), respectively. For slower dynamics, the mode of the population deviates from the true mean and the procedure described in section 3 should be applied.
( Figure 7B) which arises because only the product can be quantified. Additionally, a second ambiguity arises between the different realizations of the kinetic network in eq. (31). By considering the permutations of the kinetic scheme, different assignments of the FRET efficiencies to the static and dynamic states can result in identical FCS amplitudes, splitting the solution space into two branches (see Supplementary Note 6). For the given example of simulation 2, two competing solutions exist between the schemes LF⇌MF, HF, and LF⇌HF, MF (colored blue and red in Figure 7B). In contrast, the third permutation results in non-physical solutions for the microscopic parameters. These two solutions are indistinguishable in the analysis framework, as is evident from the identical reduced (Figure 7 B). The observed ambiguities can be described analytically based on the analytical model functions if the actual parameters are known. The resulting relation between the parameters and is shown in Figure 7C and described in detail in Supplementary Note 6.
The question remains how these ambiguities can be resolved. To decide between the two branches corresponding to the different state assignments, FRET-lines provide the required information by identifying the kinetic connectivity of the network. From the dynamic FRET-lines, we had identified LF⇌ MF, HF as the true solution ( Figure 5B and Figure 7D), allowing us to eliminate the competing solution LF⇌ HF, MF. To resolve the ambiguity arising from the presence of purely static states (eq. (32)), it is required to determine either the total fraction of dynamic molecules or the rate constant (see dashed lines in Figure 7C). For fast dynamics, is directly accessible from the twodimensional histograms as the fraction of molecules in the dynamic population that deviates from the static FRET-line. However, this approach does not apply to slower dynamics due to pseudo-static states on the static FRET-line. The fraction of dynamic molecules may also be obtained from a detailed analysis of the FRET efficiency histogram by PDA using a combination of static and dynamic populations 3 . In the given example, the knowledge of resolves the ambiguity because the interconversion rates were chosen equal ( 5 ms -1 ). However, an ambiguity remains for the general case of regarding the assignment of the interconversion rates to the dynamic states (see Supplementary Note 6 for details). Knowledge of the rate (or ) resolves the ambiguity in the analysis in all cases. To define the rate , it is sufficient to know the FCS relaxation time and the dynamic fraction . As we have shown in section 3, this information may be obtained from the mode of the dynamic distribution for intermediate to fast dynamics. For the fast dynamics of the system discussed here, one may also estimate the dynamic fraction directly from a graphical analysis in the moment representation from the position of the dynamic population along the dynamic FRET-line connection the LF and MF states ( Figure 7D). The dynamic population is positioned at the center of the line, and the resulting dynamic fractions are 0.5.
In summary, we demonstrated that, even for simple two-state kinetic networks in the presence of a background of static molecules, a global analysis of TCSPC and FCS provides ambiguous solutions. These ambiguities can partially be resolved using FRET-lines to eliminate models that are incompatible with the data. Also, it is required to know either the total fraction of dynamic molecules or the equilibrium constant of the dynamic process to fully determine the microscopic parameters.

Resolving complex kinetic networks by filtered FCS
In the previous section, we showed how binary kinetic exchange in the presence of a background of static molecules could be resolved by integrative analysis of fluorescence decays, FCS curves, and FRET-lines. To provide more challenging test cases, we performed a series of simulations with an exchange between three species in a linear reaction scheme. The possible kinetic schemes are MF⇌HF⇌LF, HF⇌MF⇌LF, and MF⇌LF⇌HF (Figure 8 B-D). Due to the fast kinetic exchange, complete averaging is observed for the dynamic population (dynF). The interconversion rates were chosen such that the FRET efficiency and fluorescence-weighted average lifetime 〈 〉 are identical for the dynamic population. Consequentially, the four scenarios are in principle indistinguishable for the two-dimensional histograms of 〈 〉 vs. (top row) or the moment representation (bottom row). While the corresponding dynamic FRET-lines indicate the correct kinetic pathways in Figure 8 B-D, the fast exchange renders it impossible to resolve the kinetic network in this case. However, the position of the dynamic population between the limiting FRET-lines of the binary exchanges is a clear indication for a three-state exchange. As described in section 3 the equilibrium fractions of the contributing species can be determined by a graphical analysis in this case.
To analyze the kinetics in the complex scenario of fast three-state dynamics, we apply filtered fluorescence correlation spectroscopy (fFCS) 6,13 . Filtered FCS exploits the information contained in the fluorescence decays to increase the selectivity and contrast of the correlation functions. By characterizing the different species in the mixture by their fluorescence decay patterns, filters are constructed that allow one to separate the contributions of the different species to the correlation function using statistical weights. Due to the orthogonality of the filters, the resulting correlation functions are species-specific. Thus, it is possible to resolve the binary exchange between different species even in complex mixtures and obtain information on the respective relaxation times . A distinct advantage of fFCS is that the number of correlation curves increases to the square of the number of contributing species. In contrast, for cFCS the number of correlation functions is limited by the number of color detection channels (four for two-color detection). However, fFCS requires prior knowledge of the number of species and depends crucially on the quality of the filters, which require precise knowledge of the fluorescent properties of each species. Given these prerequisites, fFCS can reveal the kinetic connectivity and quantify the exchange rate constant of the kinetic network.
To show the potential of fFCS, we return to the previous test cases of binary exchange between two species in the kinetic scheme MF, LF⇌ HF ( Figure 8A). For this simulation, the two-dimensional histogram reveals four peaks. Three of the four peaks (HF, MF, and LF) are located on the static FRET line, corresponding to molecules with constant fluorescence properties during the observation time. The dynamic population (dynF) is positioned on the dynamic FRET-line describing the exchange between the LF and HF populations and reveals the dynamic exchange between these species. Only the crosscorrelation function between the LF and HF species shows a positive signal ( Figure 8E). Also, the corresponding species autocorrelation functions reveal a positive correlation term that matches the timescale of the rise of the cross-correlation function (SI Figure G3). From a global analysis of the species auto-and cross-correlation functions, we obtain a single relaxation time, , which relates to the interconversion rates by 1/ → → . The cross-correlation functions that interrogate the MF⇌HF or LF⇌MF transitions, on the other hand, show no amplitude, proving that there is no exchange between these species.
For fast exchanging processes, it is not possible to resolve the kinetic network from visual inspection of the two-dimensional histograms (Figure 8 B-C). To address this problem, we computed all possible cross-correlation functions using specific filters for the three species HF, MF, and LF ( Figure 8 F-G). There is no direct connection between the LF and MF species for the linear kinetic scheme MF⇌HF⇌LF (Figure 8 F). Correspondingly, the exchange between these species is delayed compared to the direct transitions between the MF/HF and HF/LF species, as is evident from the delayed rise of the LF-MF species cross-correlation function (sCCF) compared to the HF-LF and the HF-MF sCCF. Identical observations are made for the kinetic networks LF⇌MF⇌HF and MF⇌LF⇌HF (Figure 8 G-H), showing a delayed rise of the cross-correlation for the indirect pathway. We can thus obtain qualitative information about the connectivity in the network through the relaxation times of the species crosscorrelation functions and identify dominant pathways in the network. However, it is impossible to exclude the possibility of exchange between species from the delayed interconversion alone. Similar results would, for example, be obtained for fully connected kinetic networks with a slow exchange between two species. The number of relaxation times required to describe the data provides information about the minimum number of states of the network. Since the relaxation times correspond to the inverse of the non-zero eigenvalues of the transition rate matrix, a kinetic network involving states shows 1 relaxation times in the correlation functions.
In summary, fFCS enables the direct interrogation of transitions between distinct species and allows to recover the relaxation times and the connectivity within the kinetic network. Indirect transitions show a delayed rise of the SCCF compared to direct transitions. Non-unique solutions of the joint analysis presented in section 4.1.2 are thus resolved, enabling the analysis of complex kinetic networks by smFRET. Practically, the analysis can be limited by the noise of the experimental data, which affects the quality of the resulting filters. Thus, the separation of the species and the signal-to-noise ratio of the resulting correlation functions. Another practical challenge is the identification of the contributing species for the design of the filters. Here, FRET-lines are essential to identify and assign static and dynamic species from the two-dimensional histograms. To verify that the proper solution is attained, control simulations of ambiguous solutions could be performed, enabling a direct comparison of the experimental two-dimensional histograms to that of the obtained solutions 42 .

Connecting FCS amplitudes and FRET indicators
We have shown how the combined information from FCS, TCSPC, and FRET-lines can resolve ambiguities in the analysis. The single-molecule information encoded in the two-dimensional histograms has been used to estimate the kinetic connectivity graphically, but its full potential has not been exploited. We derive the relationship between the two-dimensional histograms and the correlation amplitudes, which provides an additional restraint to the analysis. The kinetic correlation functions represent the time-dependent (co)variance of the signals. The initial amplitude at zero lag time may thus be expressed as: where Cov , is the covariance between the signals and . In the ideal case, the signals in the donor and acceptor channel are defined by the FRET efficiency: resulting in the following expressions for the amplitudes of the correlation functions: where we have used the relationship Var Var 1 Cov , 1 . Thus, the correlation amplitudes represent the normalized variance of the FRET efficiency. In the first part of the paper, we have shown that the variance of the lifetime or FRET efficiency distribution can be determined from the FRET observables and 〈 〉 by: where 1 〈 〉 is the FRET efficiency calculated from the intensity-weighted average lifetime.
This implies that we can calculate the amplitudes of the correlation function directly from the singlemolecule FRET indicators, under the condition that the dynamics are fast compared to the diffusion time. For slower dynamics, the variance is underestimated due to the limited sampling within a singlemolecule event (see section 3).
The connection between single-molecule FRET indicators and the correlation amplitudes is illustrated in Figure 9 for a three-state system with fast dynamics in the absence of static states, showing a single population that falls between the binary FRET-lines in the ( , 〈 〉 ), moment or variance representations (Figure 9 A-C). Using equations (35), a molecule-wise estimate of the correlation amplitudes is obtained that can be compared to the actual correlation amplitudes obtained from FCS analysis (Figure 9 D-G). The static and dynamic FRET-lines can likewise be converted into the equivalent of FCS amplitudes (Figure 9 D-F). The information encoded in the single-molecule FRET indicators could be used as an additional restraint in the analysis.

Conclusions
Using synthetic datasets, we challenged the capabilities of conventional analysis methods used in smFRET experiments that rely on one-dimensional data representation. In our integrative approach, FRET-lines serve as visual guidelines for interpreting experiments and the classification of populations of single molecules as static or dynamic. For slow exchange kinetics, FRET-lines directly resolve the connectivity of states. For kinetics that is significantly faster than the integration time, molecule-wise histograms, combined with FRET-lines, help distinguish dynamic averages from static populations. We developed a global analysis framework of FCS and TCSPC, which was not sufficient to identify unique solutions for two-state kinetic networks in the presence of static states. While it was possible to detect the presence of dynamics and quantify their timescale, the network connectivity and the corresponding static and dynamic fractions were not unambiguously recovered. Here, FRET-lines provided the required information to identify the limiting states of the dynamic exchange and their connectivity within the kinetic network. As a next step, the global analysis framework could be extended to utilize species-correlation function from filtered-FCS and include photon distribution analysis. We also showed that the equilibrium constant of dynamic processes could be estimated from a graphical analysis of the -〈 〉 plot even in the presence of a background of static molecules. Together with the relaxation time obtained by FCS, the microscopic rate constants can thus be quantified without requiring a precise determination of the FCS amplitudes. More complex kinetic networks consisting of three fastexchanging states could be resolved by filtered-FCS, providing species-specific auto-and crosscorrelation curves that reveal the connectivity from the patterns of the relaxation times and amplitudes. Unlike color-FCS, the number of species correlation functions in fFCS increases with the number of states, allowing robust analyses of multi-state networks.
Our global analysis framework is a step towards a self-consistent, holistic model of the two-dimensional histogram of the observables and 〈 〉 in smFRET experiments. The prediction of the moleculewise distribution of these parameters currently relies on Monte Carlo simulations, which introduce stochastic noise into the analysis that poses a problem for most optimization algorithms. Although specialized algorithms from the field of machine learning, such as evolutionary algorithms or simulated annealing, may be used to overcome the convergence problem, these algorithms require many iterations for convergence. Deterministic and efficient algorithms are therefore needed. We envision that the future of smFRET will rely on a holistic analysis of the complete experimental information, wherein the kinetic information encoded in the multidimensional histograms of molecule-wise parameters will be an essential first step for proposing candidate models that are subjected to further analysis. Using this approach, it will be possible to quantify the kinetics in complex networks, paving the way towards understanding the intrinsic dynamics of biomolecules and addressing fundamental questions relating to their function.

Supplementary Note 2 -Dependence of the dynamic shift on the mixing fractions
In this section, we derive the explicit expression for the dynamic shift as a function of the equilibrium fractions of the dynamic system. In the first part of the paper, we have defined the dynamic shift of a two-state system as the maximum distance of the dynamic FRET-line from the static FRET-line. Specifically, the dynamic shift is measured perpendicular to the static FRET-line in the parameter space defined by the FRET efficiency, , and the normalized donor fluorescence lifetime, 〈 〉 ⁄ , where 〈 〉 is the intensity-weighted average donor fluorescence lifetime in the presence of the acceptor and is the donor-only lifetime. For the difference between dynamic and static FRET-lines along the FRET efficiency axis, Δ , we obtained (compare eq. 19 in the SI of Part I): By replacing 〈 〉 using the expression for the dynamic FRET-line, (compare eq. 17 in the SI of Part I): and using the parametric expression for the average FRET efficiency as a function of the equilibrium fractions and 1 of the two states: 1 , we obtain: Δ 1 1 1 1 Then the "dynamic shift" is then given by: This function reaches its maximum at the value , : The corresponding maximum value of the dynamic shift at , is then: which is the same as eq. 28 derived in the Part I without consideration of the fraction parameter . , To obtain the joint p.d.f., , , we need to take a 3D inverse Laplace transform (ILT) of , : , ℒ , , , (11) By the change of variables ; the matrix , transfers to the matrix: The 2D ILT of , relative to expressed trough the ILT of relative to as: Where the functions ; 1; / 2√ are regularized hypergeometric functions. Thus, the joint probability of occupation times takes the form: Because of factors in the first term, this formula can be rewritten in the equivalent form in terms of a probability density function: Note that , so the second term can also be rewritten in the differential operator form: To obtain the marginal distribution of , the joint probability should be integrated over : , , , where . Making the matrix products with vectors 1 and and expressing transition rates through stationary fractions we get the pdf of occupation time : Where is the average number of transitions between two states during the observation time .
To our knowledge, Pedler was the first to obtain the same formula for the distribution of the occupation time in one state of a 2-state system 2 . Pedler has used the same 2D ILT technique as presented here but used the moment generating function of the weighted sum of occupation times proposed by Darroch et al. 3 . This formula also coincides with the one derived in Berezhkovskii et al. 4 by two methods different from ours. Kovchegov et al. 5,6 obtained the same result by another method.

Discussion
The distribution of the state occupancies is of a mixed discrete-continuous type, consisting of two discrete ( , ) and one continuous term ( ). The discrete terms describe populations of singlemolecule events in which the molecule does not leave the initially occupied state. These terms lead to the appearance of two pseudo-static populations on the 2D histograms. The corresponding discrete probabilities and tend to the equilibrium fractions as the burst duration approaches zero, and decay exponentially to zero as T increases: The continuous term of distribution describes the populations of single-molecule events in which the molecule leaves the initially occupied state at least once. This term leads to the appearance of a dynamic population on the 2D histograms. The total probability of such events becomes zero if the burst duration tends to zero and increases with increasing observation time . At → ∞ the function tends to the Dirac -function located at the equilibrium fraction of state 1 ( ).
The described behavior of the state occupancy distribution, , is illustrated in Fig. S2 for the exchange described by the transition rates constants 3, 1, corresponding to an equilibrium fraction 0.75, for different observation times.
Note that the mode (i.e. ,the most probable value) of the dynamic term, which corresponds to the maximum of the dynamic populations observed in the 2D histograms, deviates from the true equilibrium value 0.75 as the burst duration decreases. This deviation of the mode of the dynamic fraction is shown in Fig. S3 C for arbitrary transition rate constants and burst durations and illustrated in Fig. 3 of the main text. For slow dynamics (relative to the burst duration T), the mode shifts toward the limiting fraction (0 and 1) closest to the true equilibrium fraction. This deviation is absent only for 0.5. The deviation vanishes when the duration of the single-molecule events sufficiently exceeds the characteristic exchange time, i.e., when the average number of transitions is large ( 100 . Note that the sum of the transition rates is equal to the relaxation rate determined by FCS for a two-state system:

The moments over the total distribution
The mean and variance of the total distribution of the occupation time in state 1 can be derived using the moment generating function presented above (see also reference 2 ). By the change of variables / , where is the time spent in state 1, the mean and variance of total distribution are given by: Thus, the mean state occupancy is equal to the equilibrium fraction and its variance is proportional to the product of the equilibrium fractions, with a scaling factor that decays from 1 to 0 as the number of exchanges between states increases:

The moments over the dynamic part of the state occupancy distribution
Using the expressions for the total distribution (eq. (21)) and its moments (eq.(25)), we can find the moments of the state occupancy over the dynamical part ( ) only. First, the total fraction of the dynamic term is: This quantity is presented in Figure C2 A. Then, for the mean and variance of the state occupancy over the distribution / we find: , 1 1 The mean and standard deviation are presented in Fig. C2 C-D.

The mode of the dynamic distribution
The mode of the dynamic term of the state occupancy, , can be determined numerically by finding the maximum of the term . The resulting values for are shown in Fig. C2 The limiting fractions , are shown in Fig. C2 B as white curves and Fig. C3 B as a function of the average number of transitions . While the mean fraction (over the exchange term ) behaves similarly as the mode (Fig. C2 F), the slope of the dependence changes from 0 to 1 with , and the linearity is violated for intermediate values of (e.g., curves with 4, 8, 16 are obviously non-linear).
From the experiment, we can quickly determine the mode of the distribution from the peak of the distribution in the 2D histograms, from which we would like to infer the true equilibrium fraction . By rearranging eq. (29), we obtain: which enables the estimation of the equilibrium fraction from the experimentally determined value of , as described in the next section.   The equilibrium fraction is obtained by finding the intersection of the contour line of the mode with the corresponding value of . The relaxation rate determined by FCS is given by .

Qualitative description of the dynamic distribution
Three regimes can be defined based on the average number of transitions during a single molecule event, , ( ): very fast dynamics where 100; very slow dynamics where 1 ; and intermediate dynamics for 1 100. The following properties characterize these regimes: Very fast dynamics:  The distribution of the state occupancy is close to a -function.  The relative width of the distribution is <15%.  Almost all single-molecule events are located in the dynamic population ( 1). Thus, the observed static populations can be attributed to additional "truly static" species, i.e., not pseudostatic populations of a dynamic system. If it is possible to separate these static populations from the dynamic population, one can determine the total fraction of dynamic molecules, i.e., the parameter of our model. The knowledge of this parameter allows avoiding the ambiguity described in Supplementary Note 6.  The observed mode and mean of the state occupancy calculated over all dynamical population are close to the true equilibrium fractions ( ).  The kinetic exchange and diffusion contributions in the FCS curves are separated by 2 orders of magnitude. Thus, the kinetic relaxation rate can be determined accurately from FCS, and the exchange rates and can be determined by the product of the relaxation rate with the observed mode fractions , 1 , even without applying the correction described in section 3.3.

Very slow dynamics:
 The distribution of the state occupancy is highly skewed towards the edges ( 0,1 ).  The distribution stretches over the whole range of state occupancies.  A considerable fraction of single-molecule events is located in the pseudo-static populations, i.e. 0.5 and 0.5. As the dynamic distribution is broad and the pseudo-static populations are distributed narrowly, the dynamic population is challenging to detect due to its low amplitude.  The mode of the state occupancy exists only if the transition rates are close to equal, . If the rates are different, the observed mode fractions fall onto the nearest edge ( 0,1 ), i.e., to the location of the static populations. The mean dynamic fractions ( ) deviate from the center of the distribution ( 0.5).  As the kinetic relaxation time and the diffusion time are in the same order of magnitude, the corresponding kinetic and diffusional contributions to the FCS curves overlap. Thus, the reliable determination of kinetic relaxation rates from FCS becomes challenging.

Intermediate dynamics:
 If the transition rate constants are similar, , the distribution is unimodal with a pronounced maximum. On the other hand, if ≫ the distribution takes a monotonic shape and the mode of the state occupancy collapses to 0 or 1.  If the mode and the kinetic relaxation rate (from FCS) can be determined, the true equilibrium fractions can be estimated by the method described in the section 3.3.  If the dynamic and static populations are separable, the total fraction of dynamic molecules (i.e., the parameter can be estimated from the calculated value of .  In the intermediate regime, the true values of and can be estimated using the conversion of the mode into the true equilibrium fraction as described above if the difference between the transition rates is such that equilibrium fraction is within range , , 1 , .

Supplementary Note 4 -Description of FCS model 4.1 Derivation of the correlation function
The time correlation of two time-dependent variables and is defined as: , where an upper bar denotes the time average. Let suppose that observed species can be in several states (indexed by ), and each state is characterized by its pair of signal values and . We can write the observed signals and as: where is the number of species in state at time . We can write a correlation function of the total signals from N species using the definition eq.
where , , are the vectors of the corresponding quantities in the different states, is the matrix whose -th element is the correlation function of the number of species in states and and the symbol "⨂" designates the dyadic product (also called the outer product) of two vectors: ⨂ . If the species can transit from one state to another, the kinetics of this exchange is characterized by the transition rate matrix . Then, the time evolution of the numbers follows the equation: The general solution of this equation is: After sufficient time the system reaches the stationary state in which the number of species in the different states fluctuates around the stationary value .
We can thus define the fluctuation correlation functions of signals and by: ⨂ Using the result from eq. (36), the signal and fluctuation correlation functions can be rewritten as: As we consider a system in the stationary state, the averaging over time in eq. (39) can be replaced by averaging over the ensemble and the quantities , ⨂ and ⨂ acquire the meaning of first two moments and variance of the vector random variable : If the total number of observed species is stationary (i.e., if there is no diffusion of species into or out of the confocal volume), the distribution of the vector quantity is multinomial: where is the vector of the stationary fractions. The mean and variance of multinomial distribution are: Then, the second moment ⨂ can be found from the relation of the first two moments of the distribution (see eq.(40)): In the stationary state, the net flux from and to each state is zero.
Equivalently, starting from the equilibrium fractions , the system will remain in this state and thus not evolve in time, i.e.: Using this property and the expressions of the moments, the fluctuation correlation matrices can be rewritten in the form: In practice, the correlation functions are usually normalized by the product of the total values of the quantities and in the stationary state, i.e. by: Thus, the kinetic signal correlation function in the absence of diffusion takes the form: 1 , 1 1; where we have introduced the matrix ⨂ describing the variance of the fractions of the contributing species. relates to the covariance of the correlated quantities and in the stationary state by Cov , .
In the presence of a diffusion process that is slower than the kinetic exchange process, the correlation functions acquire an additional term defined by the fluctuation of the total number of species in the confocal volume. In other words, even if the kinetics reach equilibrium, the correlation function does not decay to zero. Mathematically this lead to the replacement where is the diffusion time. The parameters and are the width of the focal and the axial plane of the detection volume, respectively, where the intensity decays to 1/ of the maximum value.

The correlation function of two isolated subsets of species
Suppose that the total set of states of molecules consists of two subsets with no exchange between them, designated by the indexes d and s. (In the following, they will be considered to be dynamic or static, but for now they just represent two isolated sets.) Transitions in such a system are described by the total transition rate matrix that has block-diagonal form. The matrix exponential of a block-diagonal matrix is also block-diagonal: (50) Let and be the vectors of the stationary fractions in subsystems d and s, correspondingly. Then, the vector of the stationary fractions of the whole system and the corresponding diagonal matrix are: where is the total fraction of molecules in subsystem . As there is no exchange between the subsystems, the number of molecules (in the ensemble) in subsystems and remains constant. Therefore, the fraction is also stationary and a characteristic of the whole system, independent of the transition rate matrices and . If the correlated quantities a and b have different values in the subsystems d and s, the vectors of these values for the whole system can be written as: The correlation function of the whole system is then obtained from eq. (49) by substitution of eqs.
Then, the corresponding expression for the correlation matrix in the reduced form is: 1 where is the correlation matrix of subsystem d and is the correlation matrix of subsystem s. Thus, the correlation matrix of the whole system (i.e., the fluctuation correlation matrix for a single molecule) is the linear superposition of the correlation functions of the subsystems. The expression for kinetic part of the correlation function in eq. (49) takes the form: , Here, the matrices , and are the diagonal matrices of the corresponding fractions , and and produce the averages of the products of the correlated quantities in the final correlation functions through multiplication with the signal vectors: On the other hand, the matrices ⨂ , ⨂ and ⨂ produce the products of the averages of the correlated quantities in the final correlation functions: The normalization factor of the correlation functions can then be written as (compare eq. (47)): The matrices , and are the variance matrices of the corresponding fractions and yield the covariances of the correlated quantities in the final correlation functions: cov , As an example, cov , is the covariance of the quantities and in the stationary state of the subsystem .

The correlation matrix in the presence of static states
If the species in the subsystem are static, i.e., there are no transitions between them, the corresponding transition matrix is a zero matrix ( → ), and its matrix exponential is the identity matrix ( → ). Thus, we obtain the following expression for the dynamic system in the presence of static states: is the vector of the equilibrium fractions of the dynamic states, is the fractions of the static states, and is the fraction of molecules participating in the dynamic exchange.
If static species are absent (i.e. → 1), the correlation matrix converges to the correlation matrix of the dynamic species (see eq. (53)). Otherwise, if dynamic species are absent (i.e. → 0), the correlation matrix converges to the diagonal matrix of the static fractions . As the correlation lag time approaches zero ( → 0), the correlation matrix converges to the total fraction matrix , that is, the initial amplitudes of the correlation function are proportional to the product of the correlated quantities and averaged over all states (both static and dynamic). The asymptotic value of the correlation functions at very long correlation lag times ( → ∞ can be obtained using the eigenvalue decomposition of the matrix exponential (see the matrix defined in the next section).

Calculation of the matrix exponential using eigen-value decomposition (EVD)
The matrix is the matrix exponential of the product of the transition rate matrix and the correlation time . The matrix exponential can be expanded using the eigenvalue decomposition of the matrix : where are the eigen-matrices and λ the eigenvalues of , i.e., the solutions of the characteristic polynomial det λ 0, where det … denotes the determinate of a matrix. For a transition rate matrix, one eigenvalue is always zero, so we can set λ 0.
The matrices can be expressed in terms of the tensor (dyadic) products (⨂) of the right and left eigenvectors and : Note that eigen-matrices are orthogonal and idempotent, i.e.: The vector of the equilibrium (or stationary) dynamic state fractions is proportional to the left eigenvector corresponding to the zero eigenvalue (λ 0 and is any column (or diagonal) of the matrix : where ∀ denotes "any".
Therefore, we can write the eigen-matrix in the form: Equation (63) is convenient for computer calculations of the correlation functions, as eigenvalues and eigenvectors are effectively obtained by a numerical eigenvalue decomposition (EVD).
The alternative relation for the matrices is: This alternative expression for the eigen-matrices is more convenient for deriving analytical expressions for particular kinetic schemes (see section 4.8), as it defines the explicit relations between the elements of the eigen-matrices and the elements of the transition rate matrix.

The form of the correlation functions using EVD
Substituting eq. (62) into eq. (60) and then into eq. (48), we obtain the following expression for the kinetic part of the correlation function in the presence of static states: Note that the second form of the matrices arises from the orthogonality and idempotency of the eigen-matrices (eq.(64)).

The signal contrast factors
We require that transitions between states fulfill the detailed balance condition. The detailed balance condition states that the inward and outward flux of each state are equal, i.e., there is no net flux in the network. That is, for a given state , the sums of all inward and outward pathways must be equal: The detailed balance condition fulfills this condition: (70) In matrix notation, this condition can be expressed as: Under detailed balance, the matrices are symmetrical and each of their columns sums to zero. To prove the symmetry, that is , we obtain using eq. (67) and (68): where denotes that the detailed balance condition has been used and . designates a diagonal matrix. To prove that the columns of the matrices sum to zero, that is , we use the matrix symmetry proved above: The variance matrices introduced above have the same properties. Matrices with these properties can be decomposed in terms of the contrast matrices : ; G ; .
The elements of the contrast matrices are given by: Here, is the generalized Kronecker delta function. The product of the matrix with the correlated quantities is equal to the contrast of the states i and j, i.e., the product of the differences of and in the states i and j: Therefore, we obtain for the covariance of the correlated signals: cov , In each matrix and , there are only 1 independent elements. Let us illustrate this decomposition on the example of a three-state system, which yields the following decomposition of the covariance matrices: Note that the normalized contrast factors are not independent. The cross-correlation contrast can always be expressed in terms of the corresponding auto-correlation contrasts: . If the correlated variables and depend linearly on some common variable, but the slope of the two dependencies are opposite, then the amplitudes of the kinetic exchange, ∆ , , will be positive for autocorrelation functions and negative for cross-correlation functions (anti-correlation), as in the latter case the differences and in eq.(81) will have opposite signs. This is the case in color-FCS where ~ and ~ 1 (see section 4.7.1 below). On the other hand, if there are no dynamics in the system, both autocorrelation and cross-correlation kinetic amplitudes will vanish, and no anti-correlation will be detected. This feature can be used for the detection of dynamics in the system. In practice, for example for the fitting of experimental curves, the kinetic part of correlation functions is often described in the formal multiexponential form. Using the derivation above, we can express the exponential amplitudes in terms of physically meaningful parameters: The offset amplitudes can be expressed in terms of covariances of correlated signals: cov , cov , cov , ; The expression above is convenient to consider the asymptotic behavior of the kinetic part of correlation functions.
Using EVD, we can write the multiexponential form as (see eq.(82)): Therefore, the species correlation functions are independent of the brightness or FRET efficiencies of the states and depend only on their transition rates and fractions. This result is valid for the case that the photophysical effects are decoupled from the transitions between states.
Species-FCS functions can be obtained using filters defined by the fluorescence lifetime, anisotropy, or spectrum of the specific states. Contrary to traditional color-FCS, with filtered-FCS, it is possible to reconstruct the full correlation matrix from the species cross-correlation functions, because the optimal selection of the filters maximizes the contrast between states. In practice, a complete set of filters for each species is needed for the computation of the experimental species correlation function. The generation and the practical computation of species filters are beyond the scope of this paper and described in detail elsewhere 10,11 . Note that for fFCS, the correlation function is a pair correlation of the signal of two species-specific vectors (S a and S b ). These vectors correspond to the recorded photon stream weighted by the species-specific filters, and their time average corresponds to the fractional intensity of the respective species.

Examples
As an example, we derive expressions for the correlation functions for two-and three-state cases in the following.

Two-state case
The transition rate matrix is: The eigen-expansion of the transition rate matrix and the eigen-values (see eq.(62)) are: The stationary fractions are the columns or the diagonal of , equal to / (see eq.(67)): The eigen-matrix is obtained using eq. (67), or from the EVD above: The matrix defining the amplitudes of the exponents of the correlation functions is: The matrix above produces the product of the differences of the correlated quantities (∂ in eq.(98)).
Finally, we can account for the presence of static states and write the kinetic correlation function of the quantities and in the two-state system (see eq.(82)): In the explicit form, the correlation function is given by: Substituting the correlated signals and in the cFCS form (eq.(85)), we obtain the following expression for the contrast factors ∂ (see eq.(80)): Note that the contrast factors also depend on the total fractions of the different species through the average FRET efficiency ∑ . From eq. (100) it follows that the cross-correlation amplitude is proportional to the square root of the product of the autocorrelation amplitudes: Likewise, for the correlation functions between two species 1 and 2 by fFCS (eq.(87)), we obtain: In the case of fFCS, the kinetic part of the species cross-correlation function has the simple form: In both FCS experimental variants, the cross-correlation contrast factors can be expressed in terms of the two auto-correlation factors.
This result has important implications for the data analysis. While in total, four amplitudes can be fit, they represent only two independent parameters that can be estimated from the data.

Three-state case
The transition rate matrix is given by: The eigen-expansion of the three-state transition rate matrix is: Tr where Tr is the trace of the transition rate matrix, i.e. the sum of the diagonal elements.

⇒ ;
The equilibrium fractions are: 1 Then, the matrices defining the amplitudes of the exponents of the correlation functions are: Both the covariance matrix and the matrix can be expressed using the contrast factors: ; ; 1 0 1 0 0 0 1 0 1 ; The matrices above (one example matrix is shown) produce the product of the differences of the correlated quantities in states and : Thus, the expression for the matrices can also be expressed using the contrast factors: The complete kinetic correlation function then takes the form: Expressions for a higher number of states can be derived similarly.

Supplementary Note 5 -Ambiguity arising in the cFCS analysis
Let us consider the simplest exchange scheme -a two-state dynamic system in the absence of static states. The kinetic part of the FCS curves for such a system are given by eq. (99) by setting the parameter to 1: Phenomenologically, each FCS curve is fully described by three parameters: the relaxation rate and the auto-correlation amplitudes and , which define the cross-correlation amplitudes as . On the other hand, the full characterization of the system requires knowledge of six parameters: two transition rates ( , ) and the four components of the correlated signals ( , , , ). Therefore, the system remains underdetermined as there are more unknowns then independent experimental parameters. For cFCS experiment, the auto-correlation amplitudes are independent and defined by the FRET efficiencies of the two states. Thus, the number of parameters that can be fit is three ( , , ) while the number of physical parameters is four ( , , , ). The system remains underdetermined, and there is an infinite number of solutions (this corresponds to Ambiguity I discussed in the main text and the next section). The explicit dependence of the estimated physical parameters can be obtained by varying one of the physical parameters and solving the system of equations: The dependence of the fitted physical parameters for a particular system is illustrated in Figure S6. Thus, the analysis of FCS experiments alone does not allow to recover all physical parameters of the system, even in the simplest case. The information of the state properties (i.e., their FRET efficiencies) should be introduced into the analysis. In the case of fFCS, this is done at the step of filtering. For cFSC, the joined fitting of FCS and TCSPC can be used. However, as the number of states is rising and additional complications are introduced (for example, the presence of static states), the number of physical parameters rises faster than the number of observed phenomenological parameters, and the system becomes underdetermined again.

Supplementary Note 6 -Ambiguities arising in the cFCS+TCSPC analysis
In the analysis of the simulations, we found ambiguous solutions for fitting two-state systems in the presence of static states. In the following, we will provide a theoretical description of the observed ambiguity. Specifically, we found an infinite number of solutions that resulted in equal values for the reduced chi-square, . The solutions are undefined because the system is underdetermined, i.e., there are more parameters of the model than independent experimental observables (Ambiguity I). Additionally, the solution space is split into several branches because there is an inherent ambiguity for the system regarding the assignment of the FRET efficiencies to the individual states (Ambiguity II).

Ambiguity I -Underdetermined system
We consider a two-state dynamic system in the presence of static states of the same FRET efficiencies and one additional purely static state. We can treat such a system as a three-state system with transition rates , , , set to zero. The equations for the kinetic part of the cFCS curves can be derived from eq. (105): For now, we assume that we know which FRET states participate in dynamic exchange (states 1+2) and which are purely static (state 3). This implies that we know the correct assignment of the different FRET efficiencies to the different states, which we classify as low-FRET (L), medium-FRET (M) and high-FRET (H): Dynamic state fractions: Total state fractions: FRET efficiencies: Label: (106) The model function for such a system is described by three FRET efficiencies , three total fractions , the total fraction of dynamic molecules and two transition rate constants and . Since the fractions are normalized to one, they represent only two independent parameters of the model. We thus obtain a total of 8 independent model parameters. The analysis of the fluorescence decays obtained from TCSPC yields a set of three FRET efficiencies ( ) together with the corresponding total fractions ( ). Since the fractions sum to 1, we thus determine 5 independent parameters from TCSPC. From the FCS analysis, we obtain one relaxation time ( and three pairs of amplitudes of the correlation function of the dynamic and static contributions, and . However, as described above in section 4.8.1, the cross-correlation amplitude is related to the autocorrelation amplitudes and contains no unique information (see eq. (108) below). Additionally, the offset FCS amplitude, , is determined entirely from the quantities determined from the TCSPC analysis, that is, the FRET efficiencies and the total fractions. Thus, the only new information obtained from the FCS amplitudes is the amplitudes of the kinetic contribution of the autocorrelation functions, i.e.
A and A , given by eq.(105).
The normalized contrast factors ∂ depend only on the parameters obtained from the TCSPC analysis (i.e., the FRET efficiencies and the total state fractions and are hence known (see section 4.8.1 and eq. (100)): The cross-correlation amplitudes are related to the autocorrelation amplitudes by: Therefore, the only independent information that is obtained from the correlation amplitudes is the term 1 , which is identical for the two amplitudes A and A . Thus, we obtain only a single independent parameter from the FCS amplitudes, leaving us with a total of 7 independent observed parameters and thus an underdetermined system with respect to the 8 model parameters.
Specifically, as seen in equation (105), the ambiguity arises between the quantities and which relates to the transition rate by / . Thus, different combinations of the total fraction of dynamic molecules, , and the transition rate can result in identical amplitudes of the FCS curves, which results in ambiguous solutions subject to the boundary condition given in equation (106).
In the following, we will derive an expression to describe all possible solutions for a particular value of the FRET efficiencies and total fractions , i.e., for a given contrast factor ∂ . The fractions are functions of the transition rates , : ; 1 (109) Assuming that we know which states are dynamic (in this case states 1 and 2), the contrast factor is known from the TCSPC information and correct (∂ , ∂ , . Moreover, we know the correct value for the sum of the rates from the relaxation time of the FCS curves. Thus, we obtain the following dependency of the total fraction of dynamic molecules on the transition rate constant : Hereby, the value of the transition rate varies in the range from 0 to . An example of the possible solutions for is given in Figure S7 for equal transition rate constants of 5 (A) and 2.5 / 7.5 (B).
In practice, it may be possible to determine from other sources such as the two-dimensional histograms (see Appendix E and main text) or dynamic-PDA analysis. For the case of equal backward and forward rates (Fig. S7 A), the ambiguity is then resolved. However, in the general case of , one obtains two principal solutions for the rates which differ by the assignment of the rates to the states (Fig S7 B). This ambiguity arises because only depends on the product of the transition rate constants, which is unchanged under inversion of the state assignment.
On the other hand, one could determine the transition rate constants (or the dynamic fractions ) from the experimental data, e.g. from the mode of the dynamic distribution (see Appendix E). In this case, the ambiguity is always resolved.

Ambiguity II -State assignment
In the previous section, we have assumed that the FRET efficiencies to the different static and dynamic states are known a priori. In the experiment, this condition is not fulfilled if only the information from TCSPC and FCS is considered. However, as discussed in the main text, FRET-lines are a powerful tool to address this problem. Unfortunately, ambiguity I remains even in the case that the assignment between dynamic and static states is known because the assignment of the rates to the dynamic states remains ambiguous. The missing information about the state assignment leads to another ambiguity in the analysis. If is the number of dynamic states and is the total number of states with different FRET efficiencies, then the number of permutations to assign out of states as dynamic and as static is given by: Therefore, in principle there are , possibilities to assign FRET-states to the dynamic states. As the amplitude of the dynamic term of the correlation functions, A , depend on the contrast factor of the correlated quantities (see eq. (107) above), different state assignments may have the same covariance and correlation amplitude and may thus not be differentiable.
For the given example of two dynamic states of three FRET-states, there are 3! 6 possibilities to assign the states defined by the FRET efficiencies , , as given in Table E1. The permutations 1↔1', 2↔2' and 3↔3' differ only by the assignment of the FRET efficiencies to the dynamic states, while the FRET efficiency of the purely static is unchanged. Due to the uncertainty introduced by ambiguity 1, they have the same correlation amplitudes and are thus indistinguishable.
On the other hand, the state assignments 1, 2 or 3, i.e. the choice of the dynamic states to be L+M, L+H or M+H, results in different amplitudes . In principle, there could thus be three possible solutions for the system based on the state assignment. Not all of these solutions are physical due to the boundary conditions: The possible solutions additionally have to fulfill the condition: Thus, there is an ambiguity of the recovered parameters of the dynamic states caused by the fact that the contrasts factors ∂ defined by the different connectivity of states can be compensated by the factor defined by the different transition rates such way that resulting amplitude of correlation curves is the same. This is described by the following condition (compare eq. (110)): where ∂ is the contrast factor of the correct state assignment. Thus, we obtain three branches of solutions for the three different assignments: Where describes the general dependency given by eq. (116). The relation can be further simplified by considering that the contrast factors have the same normalization (see eq. (107)) to: For the systems described in the previous section where the true exchange is between the low-FRET and medium-FRET states (eq. (106) and Fig. F2), an additional solution is found by assigning the dynamic exchange to the low-FRET and high-FRET states (Fig. S8 A and D). The third solution for exchange between the medium-FRET and high-FRET states is non-physical and results in values of 1. Similar solutions are obtained for true exchange between low-FRET/high-FRET and medium-FRET/high-FRET states (Fig. S8 B-C and E-F), whereby all three state assignments result in physical solutions the MH case.
In the case of fFCS, the ratio of the contrast factors ∂ is equal to 1 (see eq. (101). Thus, the ambiguity of type II disappears, however the ambiguity I remains. In conclusion, the fitting of FCS together with TCSPC does not allow to recover the kinetic parameters even for a simple 2-state exchange when static species are present.

Simulations of freely diffusing single molecules
Simulations of single-molecule measurements were implemented using a Brownian dynamics approach as described previously 12 . The spatial intensity distribution of the observation volume was modeled by a 3D Gaussian distribution. To simulate the kinetic exchange, the time that a molecule spends in state is assumed to be exponentially distributed: where the rate represents the sum of all rates depopulating state : Simulated photon counts are saved in SPC-132 data format (Becker & Hickel GmbH, Berlin, Germany) and treated as experimental data. Table S1: Fluorescence parameters of the different FRET species used in the simulations for the high-FRET (HF), medium-FRET (MF/MF') and low-FRET (LF) species. Given are the FRET efficiency E, the intensity-weighted average donor fluorescence lifetime ⟨ D(A) ⟩ F and the brightnesses of the donor fluorophore and the FRET-sensitized acceptor fluorophore Q G and Q R . The MF species has different parameters in simulations 1-3 compared to simulation 4. In simulation 1, an additional MF' species is included. In all simulations the donor-only lifetime was  0 = 4 ns, the donor-only brightness was Q D0 = 50 KHz, the Förster radius was R 0 = 52 Å, the quantum yield of the donor and acceptor were  D =0.8 and  A =0.32, the spectral crosstalk was =0.017, the diffusion time was t =5 ms and the detection efficiencies ratio between the green and red detection channels was g G /g R =0.8.  Table S2: Kinetic parameters and species fractions of the different FRET species used in the simulations. Given are the average number of molecules in the confocal volume N tot , the species fractions of the different species , the total fraction of dynamic molecules , the states involved in dynamic exchange and the respective interconversion rates. The total fraction of a given species is related to the total fraction of dynamic molecules, , and the static and dynamic fractions of the species, and , by 1 . The MF species has different parameters in simulations 1-3 compared to simulation 4 (see Table S1). In simulation 1, an additional MF' species is included. In all simulations the donor-only lifetime was  0 = 4 ns, the donor-only brightness was Q D0 = 50 KHz, the Förster radius was R 0 = 52 Å, the quantum yield of the donor and acceptor were  D =0.8 and  A =0.32, the spectral crosstalk was =0.017, the diffusion time was t d =5 ms and the detection efficiencies ratio between the green and red detection channels was g G /g R =0.8. . As a fit with three FRET species results in identical compared to the ground truth model, it is impossible to identify the number of states correctly using only the TCSPC information. B) Color-FCS autocorrelation and cross-correlation functions reveal no kinetic contribution as is visible from the plateau of the correlation functions up to the timescale of diffusion ( = 5 ms). C) Photon distribution analysis (PDA) reveals the four FRET species and correctly recovers the distances. Weighted residuals (w. res.) of the fits are given above. The simulation parameters are given in Table S1 and S2. . Weighted residuals of the fits are given above. The simulation parameters are given in Table S1.