Resolving fluorescent species by their brightness and diffusion using correlated photon-counting histograms

Fluorescence fluctuation spectroscopy (FFS) refers to techniques that analyze fluctuations in the fluorescence emitted by fluorophores diffusing in a small volume and can be used to distinguish between populations of molecules that exhibit differences in brightness or diffusion. For example, fluorescence correlation spectroscopy (FCS) resolves species through their diffusion by analyzing correlations in the fluorescence over time; photon counting histograms (PCH) and related methods based on moment analysis resolve species through their brightness by analyzing fluctuations in the photon counts. Here we introduce correlated photon counting histograms (cPCH), which uses both types of information to simultaneously resolve fluorescent species by their brightness and diffusion. We define the cPCH distribution by the probability to detect both a particular number of photons at the current time and another number at a later time. FCS and moment analysis are special cases of the moments of the cPCH distribution, and PCH is obtained by summing over the photon counts in either channel. cPCH is inherently a dual channel technique, and the expressions we develop apply to the dual colour case. Using simulations, we demonstrate that two species differing in both their diffusion and brightness can be better resolved with cPCH than with either FCS or PCH. Further, we show that cPCH can be extended both to longer dwell times to improve the signal-to-noise and to the analysis of images. By better exploiting the information available in fluorescence fluctuation spectroscopy, cPCH will be an enabling methodology for quantitative biology.


Introduction
Absolute numbers are necessary for quantitative biology, but are challenging to obtain in single cells. Although ad hoc techniques exist [1][2][3][4][5][6], an established approach to measure the concentrations of fluorescently tagged molecules and their rate of diffusion is fluorescence fluctuation spectroscopy, which has two main branches. Fluorescence-correlation spectroscopy (FCS) uses correlations between the observed photon counts at different times to obtain a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 both the diffusion coefficient and the local numbers of a fluorescent molecule in the illuminated observation volume [7]. Analysis of photon-counting histograms [8,9] and related methods of moments [10,11] exploit fluctuations in the amplitude of the photon counts to estimate the brightness of a fluorescent molecule (and so allow an estimate of absolute numbers and the degree of oligomerization through the number of fluorescence tags per molecule). Nevertheless, information on both diffusion and brightness are often needed [12,13]. Although it is possible to employ PCH and FCS on the same data [12][13][14], the two techniques are not inherently linked, either making interpretation complicated or requiring global fits across multiple experiments to extract the parameters of interest.
Techniques developed to be sensitive to both diffusion and brightness do exist, but often are either computationally intensive or employ approximations that can limit their use. For example, higher-order FCS [15] does not consider the statistical uncertainty of the higher order moments and even later implementations are not able to fully exploit the information available on the amplitude of the photon counts [16]. To improve signal-to-noise, PCH was extended to consider longer dwell times in the observation volume [17][18][19]: molecules then have time to move to regions with different intensities of illumination so that their diffusion affects the data. By considering many dwell-times, these techniques become sensitive to both diffusion and brightness. Yet the PCH methods only approximate the dependence on diffusion, which can lead to errors in the estimation of the higher order moments [11,20], and any alternatives are computationally intensive [18,21,22].
Here we present a method that is not only sensitive to both brightness and diffusion, but also unites FCS, PCH, and moment analysis. The new technique, which we call correlated photon-counting histograms (cPCH), generalizes dual-color PCH [23] so that the two channels being compared can now be shifted in time. Like PCH, the histograms depend on brightness, but how these histograms change from one time point to another depends on diffusion. The histograms can be generated using data either from the same channel or from two different channels, and, although we will focus on single channel measurements, the theoretical expressions derived apply generally.
We also show that cPCH can be applied to the study of images, just as raster image correlation spectroscopy (RICS) extends FCS [24] and spatial intensity distribution analysis (SpIDA) extends PCH [25]. Spatial cPCH can be used to resolve combinations of both immobile and mobile proteins, and like cPCH can resolve different species based on their brightness and diffusion properties. We demonstrate its potential using a simulated system of a fluorescent ligand binding to a receptor with two distinct binding sites.

Deriving cPCH
To maintain sensitivity to both the brightness and diffusion of the fluorophores, we need to simultaneously employ information on both the correlation and amplitude of the photon counts. FCS-like methods are sensitive to diffusion through measuring correlations between fluctuations in the fluorescence signal over time; PCH-like methods are sensitive to brightness through analyzing fluctuations in the amplitude of the signal. Our approach is to calculate histograms of the number of photons detected in time bins that occur a time τ apart.
We allow each channel to be stimulated by a different laser (the observation volumes in each channel can therefore be of different sizes) and refer to the two channels as channel A and channel B. Each diffusing molecule has a different brightness profile in each channel, I A (x) and I B (x), with x representing a three dimensional coordinate. We assume that the sampling time T is smaller than τ d , the time to diffuse through the observation volume, so that molecules can be considered stationary during each sampling event. For a stationary molecule, the emission and detection of photons in each channel will be uncorrelated (shot noise), and each will follow a Poisson distribution with an average rate that depends upon the fluorophore's position in the observation volume.
We aim to calculate p N (n x , n y |τ), the probability of detecting n x photons at time 0 and n y photons at a time τ later for a sample containing on average N molecules in the observation volume. To begin, we consider the case of a single molecule diffusing in the observation volume [8,26] and assume that the binning time T is much smaller than the diffusion time τ d of the molecules through the observation volume.
For a single molecule, the probability of interest is p 1 (n x , n y |τ). As we do not know the position of the molecule, we integrate the probability of detecting n x photons from a molecule at position x and n y photons at a position y over all possible positions: p 1 ðn x ; n y jtÞ ¼ Z p 1 ðn x ; n y ; x; yjtÞdxdy assuming that shot noise is uncorrelated. The term p(x, y|τ) in Eq 1 is an important difference between this equation and the singlemolecule probability functions for PCH and dual-color PCH [26]. This term defines the probability of being in position x at time 0 and position y at time τ. For free diffusion (without boundaries) and if D is the diffusion coefficient of the fluorophore, then pðx; yjtÞ ¼ pðyjx; tÞpðxjtÞ where the probability p(x|τ) is the probability density of finding the molecule at x at time τ and is 1/V, where V is the experimental volume (which is larger than the observational volume V AB ). Rather than performing the integration in Eq 1, we instead find a solution using Taylor series. If the probability of obtaining n photon counts from any given position depends only on the illumination intensity at that position, I(x), then Eq 3 implies that Eq 1 contains exponentials of I(x) + I(y), which we Taylor expand and then use the binomial theorem to expand further into terms in I(x) and in I(y). Reversing the order of integration and summation gives Without the diffusion term, p(x, y|τ), the integral in Eq 4 gives the moments of the brightness profile. Eq 4 therefore reduces to a series solution for the single-molecule case of dualcolor PCH. With the diffusion term, the integral gives a higher-order two-point correlation function, G m,n (τ): where V AB is the volume of the observation region (found from the point spread function). Consequently, Eq 4 can be written as G n x þkÀ m;n y þm ðtÞ: ð6Þ To find explicit expressions for the two-point correlation functions (Eq 5), we use a Gaussian profile for the intensity of illumination [16,18,27]: where r A and z A denote the radial and axial beam waist parameters, and � A is the constant brightness of the single molecule in channel A in units of photon counts/bin (of duration T seconds). A similar expression holds for I B (x). Eq 7 implies that the two-channel observation volume is (as in fluorescence cross-correlation spectroscopy [28]): Independently treating each spatial dimension in Eq 5 gives where κ m,n (τ, D) is The shape factors of the observation volume are defined as and obey for the Gaussian intensity profile. If the two channels have the same observation volume (A = B), Eq 10 becomes the correlation function derived by Melnykov et al. [16]. The shape factors simplify to γ m,n = γ m+n and become the shape factors for a single channel (g n ¼ n À 3 2 for a Gaussian intensity profile). Extending to N molecules. We have considered the probability of detecting photons from single molecules, but the molecules of interest will diffuse in and out of the laser's observation volume, and the actual number of molecules in the observation volume at any given time will fluctuate around the average value N. These fluctuations are also driven by a Poisson process and so the probability of detecting n x photons from channel A and n y photons from channel B will be given by a compound Poisson distribution.
First, we briefly discuss a technical point. The fraction V AB /V in Eq 6 is the ratio of the observation volume to the larger experimental volume, but can be set to unity when we consider more than one molecule in the observation volume because this ratio then cancels in the generating function [29]. Although we should consider the number of molecules in the larger experimental volume, N v , we will assume V AB = V and employ N = N v throughout because a constant concentration of molecules implies N v /V = N/V AB .
The probability of detecting n x photons in channel A and n y photons in channel B at a time τ later when there are N molecules on average in the observation volume can be found by marginalizing the probability of detecting those photons given there are m molecules exactly in the observation volume multiplied by the probability of detecting these m molecules if there are N in the observation volume on average: where p(m|N) is given by a Poisson distribution, The probability p(n x , n y |τ, m) is given by the m-fold convolution of the single molecule probabilities. The probability density of a sum of independent variables (n x and n y are both sums over the photons emitted by each molecule) is the convolution of the individual probability densities [30] (these densities determine the individual n 0 x and n 0 y emitted by each molecule). Therefore we have pðn x ; n y jt; mÞ ¼ p 1 ðn x ; n y jtÞ � p 1 ðn x ; n y jtÞ � . . . � p 1 ðn x ; n y jtÞ ðm timesÞ: To proceed, we use the bivariate probability generating function given by: t n x s n y p N ðn x ; n y jtÞ ð16Þ and note that the m-fold convolution of a probability distribution is equivalent to raising its generating function to the m'th power [30]. If g 1 (t, s|τ) is the generating function for the single-molecule probability distribution, we can write t n x s n y X 1 m¼0 pðn x ; n y jt; mÞpðmjNÞ We have an exponential generating function in Eq 17, and we will show that this property translates into a sum of the different single molecule generating functions, weighted by the concentration of each species, because the generating function of a sum of independent random variables is the product of their individual generating functions [30]. If there are w independent, fluorescent molecular species in the observation volume (which do not have substantial interactions over the time scale of the experiment), the generating function can be written as For better comparison with experiments, we also include a constant (Poisson) background noise in each channel. This noise has an average of λ A counts per bin time T in channel A and λ B counts per bin time T in channel B. The generating function is then Although later we will derive a recursion relation for the cPCH probability distribution from Eq 19, it is often more efficient to use the Fourier transform to evaluate the generating functions (replacing t and s in Eq 16 by e it 0 and e is 0 [30]) and directly calculate the cPCH distributions using the inverse Fourier transform, F À 1 : where the second term describes the background noise. The generating functions for each single molecule probability distribution, g ðqÞ 1 ðt 0 ; s 0 jtÞ, are calculated using the Fourier transform of the single molecule probability distributions p 1 (n x , n y |τ). The joint probability in the background noise is and, using the properties of the discrete Fourier transform, where L s 0 and L t 0 are the numbers of elements in s 0 and t 0 . We note that different species are included in Eq 20 by summing their single-molecule generating functions weighted by their concentrations rather than the successive convolutions of the PCH approach [8,26]. Recursive solution. The probability generating function can also be used to derive the probabilities. By definition [30] pðn x ; n y Þ ¼ lim and writing p 1 i ðn 1 ; n 2 jtÞ for the single-molecule probability, Eq 1, for species i, we find the following recursion for the cPCH probability distribution: • for n x > 0 and n y > 0: • for n x > 0 and n y = 0: • for n x = 0 and n y > 0: • for n x = n y = 0: To calculate the recursion we start at n x = n y = 0, followed by n x = 0 and n y = 1, and then by n x = 1 and n y = 0 before continuing with successively higher values of n x and n y .
Eq 25 is also applicable to dual color PCH if τ = 0. Analytical expressions for single-channel PCH can be derived similarly [31,32].

Moments
Similar to PCH [11,27], the moments of the cPCH distributions efficiently summarize the relevant information. The moments are easier to compute than the full distribution, and we can express all of the moments, and even the single-molecule probability distributions, in terms of the factorial moments of the single-molecule distributions.
For moments of order m for n x and n for n y , we use the notation: M m,n for the raw moment (about zero); (M) m,n for the factorial moment; K m,n for the cumulant; and (K) m,n for the factorial cumulant. All these quantities can be found through generating functions. As well as the probability generating function, Eq 16, the moment generating function is required: e tm e sn pðm; nÞ: The raw moments, M m,n , and cumulants, K m,n , are found from this moment generating function. Differentiating g (M) (t, s), m times with respect to t and n times with respect to s, and then setting s = t = 0 gives the raw moments; differentiating the logarithm of g (M) (t, s) and setting s = t = 0 gives the cumulants. The factorial moments, (M) m,n , and cumulants, (K) m,n , are found from the probability generating function, Eq 16. Differentiating g(t, s), and setting s = t = 1 gives the factorial moments; differentiating the logarithm of g(t, s) and setting s = t = 1 gives the factorial cumulants.
Factorial moments for a single molecule. Using Eq 3 in Eq 1 and the definition of the probability generating function, Eq 16, we find that g 1 (t, s|τ) is after performing the sums over n x and n y . Differentiating Eq 30, we find that the factorial moments of the single molecule distribution are with G m,n (τ) being given by Eq 9. We note that comparing Eq 6 and Eq 31 allows the single-molecule probability distributions to be expressed as an infinite sum over the factorial moments: ðMÞ n x þkÀ m;n y þm ðtÞ: ð32Þ Factorial cumulants for N molecules. For N molecules, the factorial cumulants provide the simplest expressions: The factorial cumulants for the N molecule case are therefore the higher-order two-point correlation function, Eq 5, multiplied by the average number of molecules N.
We can extend these results to w different species by noting that the logarithm of the generating function for the multiple molecule case, Eq 18, is the sum of the individual generating functions weighted by their concentrations. We can therefore sum the factorial cumulants of each individual species. For better comparison with experiment, we include a Poisson background noise of λ A counts per bin in channel A and λ B counts per bin in channel B by noting that only the first factorial cumulant of a Poisson distribution is non-zero (and is identical to the mean of the distribution). We find for m, n > 1. Here κ m,n and γ m,n are defined in Eqs 10 and 11. Further, we can derive recurrence relations to convert the factorial cumulants into cumulants and raw moments, which is necessary for comparing theory with experimental results (Sec. S1).

Moments of moments
The method of moments of moments is used to estimate the experimental uncertainty of higher order factorial cumulants [11]. Rather than using k-statistics [27], we note that the variance of a random variable can be expressed as a Taylor series of its raw moments [33]. To first order in N d −1 with N d being the number of data points (typically on the order of 10 6 ), we can write the variance of either the cumulant or factorial cumulant of interest, X p,r , as where Cov[M x,y , M u,v ] is the covariance of the two raw moments M x,y and M u,v and is given by [33] Using the recurrence relations in Sec. S1, we can express the factorial cumulants in terms of the raw moments and so evaluate the derivatives in Eq 36. We provide examples of variances of the factorial cumulants in Sec. S9.
Finally, we note that we have assumed that the photon counts are independent of each other, which is not true in practice, and a correction is required for the variance of the first factorial cumulant [18].
Obtaining previous methods in fluorescence fluctuation spectroscopy FCS. We can calculate an FCS curve from the moments of the cPCH distribution. The normalized autocorrelation function of FCS is: or using the definitions of the raw moments of the cPCH distribution. Via the relations in Sec. S1, we can also express Eq 38 in terms of the factorial cumulants: Fitting FCS curves requires the error at each point of the curve [34]. With cPCH, not only do the analytical solutions to the factorial cumulants provide expressions for fitting, but also the the variances calculated by the moments of moments technique estimate the errors. In addition, we can include dead-time and afterpulsing effects in the cPCH distribution and so also in the FCS curve [35].
Care is needed when τ = 0. Then n x = n y by definition because we are comparing each time bin with itself: the empirical probability distribution p(n x , n y |τ = 0) will be zero everywhere except along the diagonal where n x = n y , and the data are one (not two) dimensional.
For FCS, we therefore define replacing the relationship for the joint moment by the equivalent relationship for the second moment of a one dimensional distribution when τ = 0. PCH. We can obtain a single-channel PCH at any τ by summing over all counts on one of the two axes of the cPCH distribution, p N (n x , n y |τ), because p N (n x ) = p N (n x |τ) and p N ðn x jtÞ ¼ P n y p N ðn x ; n y jtÞ. For example, using the single-molecule probability distribution p 1 (n x , n y |τ) (Eq 3 in Eq 1): 3 2 e À ðyÀ xÞ 2 4Dt dxdy carrying out first the sum over n y and then the integral over y. This integral over x is the definition of the single-molecule probability distribution used in PCH [8].
Obtaining a PCH curve at any τ is useful: in cPCH, it can be beneficial to use larger bin sizes T as τ increases to increase the signal-to-noise ratio, a common practice in FCS. Using Eq 42, it is therefore straightforward to obtain PCH curves at different bin sizes, p 1 (n x |T), increasing the effective counts/bin to improve the signal statistics, following FIMDA [17,19]. Using several different T also provides additional constraints on the brightness in the fitting routines. As in FIMDA, the effect of molecular diffusion during longer bin times on the cPCH distributions must also be included and will be discussed below. To calculate the theoretical p 1 (n x |T), we can take diffusion into consideration during the bin time T by using Eq 51 for the factorial moments in the single molecule series solution (Eq 32).

FCA.
We can obtain the fluorescence cumulants used in FCA [11] by setting m (or n) to zero in Eq 33: ðKÞ n ðtÞ ¼ N� a n g n;0 : We note that the definition of γ m,n in Eq 12 includes any differences in the size of volume between the two channels so that (K) n,0 (τ) will yield the same numerical value as (K) n from FCA.
Dual-color PCH. The dual-color PCH distribution is the dual-color cPCH distribution with a zero time-shift:

Diffusion effects at large bin sizes
We have assumed that the binning time T is smaller than the diffusion time τ d of the molecules, but, as T increases towards τ d , molecules can no longer be considered stationary in the observation volume and will experience a range of intensities, I(x), changing the statistics of the photon counts. FIMDA and TIFCA. To allow larger T, we follow Fluorescence Intensity Multiple Distribution Analysis (FIMDA) [17] and introduce a second order binning function B 2 (T, τ d ) so that the factorial cumulants have the same form as Eqs 33 and 43, but with effective brightness and concentration parameters. These effective parameters are: where � o is the count rate of the fluorescent molecule per second rather than per bin time T, and B 2 (T, τ d ) is [17]: bð1 þ aÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À b p atanh ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 À b p ð ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 þ ba p À 1Þ ðb þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 þ ba p À 1Þ assuming a 3D Gaussian observation volume and with β = (r a /z a ) 2 and α = T/τ d . Moment analysis was also extended to include binning effects using Time Integrated Factorial Cumulant Analysis (TIFCA) [18]. All factorial cumulants above first-order are weighted by a binning function of similar order, i.e., which should be compared with Eq 43. Although the higher order binning functions in TIFCA are exact, rather than the approximations used in FIMDA, binning functions with n > 2 require more complex numerical integration. By comparing FIMDA's approximation with the TIFCA results, we can show that FIMDA's approximation for the higher order binning functions is Replacing the concentration and brightness parameters in the factorial cumulant expressions (Eq 43) with FIMDA's effective concentration parameters, we find that from which follows the approximation in Eq 48. cPCH. Deriving a general equation for B m,n (T, τ d ) is challenging, and we instead use FIM-DA's approximation (Eq 48). From Eq 32, the cPCH probability can be expressed as an infinite series of the factorial moments (M) m,n (τ, T).
We separately consider the order for each channel: the first channel is of order m and the second of order n, and binning is only required if the order of a given channel is either 2 or higher. Each moment will then be scaled by its binning function B m,n (T, τ d ). Although dual color PCH [26] and TIFCA [27] (equivalent to cPCH when τ = 0) use a second order binning function for K 11 (T), we have found that this approach applied to single color distributions matches neither simulated nor experimental data for τ > 0.
For cPCH, we thus approximate the binning functions using Eq 48:

Extending to flow and images
We can extend cPCH to the study of images. Correlations between different regions of an image should only occur if the same molecule has visited both regions. We focus on detectors run in photon-counting mode and raster-scanned images. We therefore need to include the probability of detecting the same molecule at two different positions at two different times. This movement of the scanning laser is analogous to a directed flow and, similar to rasterscanned image correlation spectroscopy [24,36], we include a flow term in the diffusion part of the correlation function (Eq 9). This change is the only one needed.
If there is a flow v x in the x-direction and v y in the y-direction, we can write the correlation function as where κ m,n (τ, D) is the diffusion term of Eq 10 and S m,n (τ, D) describes the contribution from flow: assuming a 3D-Gaussian observation volume in each of the two channels.
The times τ that we can consider depend on the time taken by the laser to travel from one pixel to another. To translate the distance between pixels into a time [24], we use the lineretracing time, τ L , the pixel dwell time, τ p , the pixel size, δ x by δ y , and the number of pixels in a line, L y , so that v x τ = d x δ x and v y τ = d y δ y in Eq 55 for a pixel located d x pixels away in the horizontal direction and d y pixels away in the vertical direction. The time for the laser to reach the position d x in the current line is τ p d x ; the time for the laser to reach the position d y by scanning the previous lines is τ p L y d y ; and τ L d y is the corresponding line-retracing time. Consequently, τ in Eq 55 becomes τ p (d x + L y d y ) + τ L d y .
With Eq 53 for G m,n (τ), we can calculate the probabilities p N (n A , n B |d x , d y ) and the associated moments as before (for example, using Eq 6). We note that for images of cell membranes, 2D diffusion is more appropriate, and Eqs 8, 10 and 11 should be replaced by their 2D equivalents.

Triplet states and detector dead-time and afterpulsing
As with other methods in fluorescence fluctuation spectroscopy, cPCH can be modified to correct for non-ideal behaviour from fluorophores and detectors. We correct for triplet states and fluorophore blinking following FCS [37] and FIMDA [17] (Sec. S6) and for detector dead-time following Melnykov and Hall [16] (Sec. S7). For afterpulsing by the detector, we build on previous work [38], but develop corrections that include the probability of detecting an afterpulse at a time τ after the detected photon (Sec. S8). This approach allows us to remove the increase in amplitude due to afterpulsing commonly seen in FCS at small correlation times.

Simulations
We simulated both the diffusion of different species of fluorescent molecules through the observation volume and the photon emission process using two different simulators.
Continuous space and discrete time. To simulate experiments when the binning time is much less than the diffusion time, we developed a bespoke simulator in Matlab (Mathworks, MA), which uses a discretized time but a continuous representation of space to better describe the observation volume. Diffusion within time bins does not affect the signal statistics because photons from all molecules are generated at the end of a time bin once their positions have been updated. Given an initial concentration, the simulator randomly places molecules in a box of x = 4μm by y = 4μm by z = 8μm with the laser focused at the box's centre. At each Δt, we sample the emitted photons for each molecule i in the observation volume from a Poisson distribution with mean I i (x i ). We then independently update the x, y, and z positions of each molecule using Gaussian random variables with variance 2D i Δt. The number of molecules in the volume remains constant: if any of the molecules extend beyond the boundaries of the simulation, they are wrapped to the other side of the simulated volume. We consider two detection channels, with each species potentially having a different brightness in each channel.
Discrete space and continuous time. To simulate experiments when the binning time approaches the diffusion time, we modified the reaction-diffusion simulator MesoRD [39]. Molecules can diffuse and fluoresce at any time, and so effects such as diffusion within a bin time, saturation of the excited state, triplet states, detector dead-time, and detector after-pulses can be explored. The simulator uses a discrete space, but continuous time representation. We therefore discretize the intensity profile, and the extent of this discretization affects accuracy.
We split the simulation volume into two compartments: an inner compartment containing the observation volume and an outer compartment for which the laser intensity is low. Each compartment is further divided into thousands of subvolumes (cubes of 25-33 nm per side). Each subvolume has a laser intensity dependent on its position relative to the center of the simulation volume. Following their concentrations, molecules are randomly distributed among the subvolumes. Each species undergoes a set of molecular reactions, such as excitation to the singlet state and relaxation to the ground state. Some reactions result in a photon being emitted and measured in a specified detection channel if a random number generator produces a number smaller than the channel's detection efficiency. In the outer compartment, molecules undergo a birth-death process, and the number of molecules in the simulation volume fluctuates. Typically, data was recorded in arrival-time mode, with 50 ns resolution and with a simulation time of 60 s. Several simulation runs (usually 5) are used.
We tested the quality of the simulations by fitting the data, similar to a calibration experiment, to ensure that we obtained the correct concentrations, brightness values, diffusion times, structure factor (s = z A /r A ), and shape parameters (γ 3 , γ 4 ).
Images. We extended our simulator to simulate photon detection using a raster-scanned laser scanning microscope. For a Xμm by Yμm image, the simulation volume is a region of x = 4Xμm by y = 4Yμm by z = 8μm. Molecules are placed randomly, and the laser is initially positioned at the top left corner: the region (−X/2: X/2, −Y/2: Y/2, 0) of the simulation volume. Each image is assigned x n by y n pixels. All pixels have a defined dwell time. At each dwell time, we count photons from all molecules in the simulation volume with the point-spread function centred at the current pixel. These photon counts are entered into the corresponding pixel of the current image frame. The position of each molecule is then updated by diffusion, and the laser moves to the next pixel. At the end of each line, the laser shifts to the next line, but the molecules first have their positions updated over the line retracing time (although in multiple steps with each step having the duration of the laser dwell time). Similarly, the molecules' positions are updated over the frame retracing time at the end of each frame. The updates are also performed in multiple steps to allow boundary conditions to be considered at each step (molecules that have moved outside the simulation volume are wrapped to the other side).

Calculating the experimental cPCH
To extract parameters, we must fit the experimental data.
Formatting data for cPCH is straightforward if the data is binned using a single time bin, T. For each τ, we construct a two-dimensional histogram with photon counts forming one axis of the histogram and photon counts observed a time τ later forming the other. Systematically considering all time bins separated by τ then allows us to fill the histogram's bins. The empirical histogram of photon counts for cPCH, denoted p � N ðn a ; n b jt; TÞ, can then be converted to the empirical probability distribution, p N em ðn a ; n b jt; TÞ by normalizing so that its sum is unity. This method can, however, become memory intensive at the smaller bin times of around 200 ns used in FCS.
Using arrival-time data. We use arrival-time data to apply cPCH for all the time ranges employed in FCS. It is simplest to work with bin times that are integer multiples of the detector's sampling time.
Given a list of photon arrival times t, we first convert the arrival times and the desired bin time, T k , into integer multiples of the sampling time, i.e. to t s and T s k . We then divide the integer arrival times by the integer bin time, rounding up to the closest integer, to give a list b, which indicates in which bin time photons arrive. For example, if t s = {1, 2, 3, 6, 11, 17, 24, 46} and T s k ¼ 10, then b = {1, 1, 1, 1, 2, 2, 3, 5}, indicating that there were 4 photons in the first bin time, 2 in the second, 1 in the third, 0 in the fourth, and 1 in the fifth. Numerically, this formulation allows the accumulation property of sparse matrices to be used to determine the total photons detected per bin. For the above example, we have that the bins with positive counts are b + = {1, 2, 3, 5} with corresponding counts of total photons of {4, 2, 1, 1}. To find the bins at time τ, we normalize τ by the bin time T k to obtain τ k . We define b þ t as b þ t;i ¼ b þ iþt k . To determine p � N ðn a ; n b jt; TÞ, we first find p � N ðn a > 0; n b > 0jt; TÞ by constructing a histogram for all pairs of bins with non-zero photon counts from the intersection of b + and b þ t , with the first element of the pair, n a , from b + and the second element, n b , from b þ t . Second, when n a = 0 and n b > 0 and when n a > 0 and n b = 0, we construct a histogram of the photon counts for bins where b + and b þ t do not intersect: for example bins that are in b + but not in b þ t have n a > 0 and n b = 0. Finally, for n a = n b = 0, we count the number of bin pairs with zero photon counts, which is the difference between the total number of bin pairs and the number of bins with a positive photon count in at least one channel. The empirical probability is calculated by normalizing p � N ðn x ; n y jt; TÞ by the total number of bin pairs.

Calculating the fit energy
We calculate the error of a particular model for a given set of parameters similarly to PCH [8] and model each pair of photon counts {n x , n y } as a binomial experiment with a chance of success of p N (n x , n y |τ) and a chance of failure of 1 − p N (n x , n y |τ). Fitting the cPCH distribution. If N d is the number of data points at a particular τ and a particular n x and n y , the standard error, σ p , for the estimate of p(n x , n y , τ) obeys where p N em ðn x ; n y jtÞ is the empirical probability distribution. We note that Eq 56 ignores any correlations between different segments of the dataset. Writing L n x ;n y ðtÞ as the total number of histogram entries observed at a given τ, we define the fit energy-or loss function-for the cPCH distribution as where we have summed the normalized square error of each entry of the histograms over all τ.
If the model fits the data within the empirical variance, we expect the energy value to be close to one. Fitting the factorial cumulants. It is often more efficient to fit the factorial cumulants of the data rather than use Eq 57. First, we find the moments M m,n (τ) from the empirical cPCH distribution using M m;n ðtÞ ¼ P n x P n y n m x n n y p N em ðn x ; n y jtÞ. Second, we convert these moments into the cumulants and then into the factorial cumulants (Sec. S1). Third, we determine the empirical variance of each measured factorial cumulant, s 2 K ðtÞ, using moments of moments (Eq 36). Specific expressions, which need only be calculated once per dataset, are in Sec. S9 and are in principle identical to those for dual color FCA and TIFCA [27]. Writing L m,n (τ) now as the total number of moments used at each τ, the fit energy is which is the sum squared error of each factorial cumulant for each τ normalized by the total number of moments used.
Nevertheless, the analytical solutions do provide an efficient method to sample brightness and concentration values that are within the experimental error of the data. We use the method of moments of moments to estimate the error (variance) for each factorial cumulant in terms of the empirical raw moments following Eq 36. Once this variance is estimated and assuming the empirical factorial cumulants obey a Gaussian distribution with a mean equal to their empirical value and a variance given by Eq 36, we can resample each moment within the expected experimental error and use Eqs 60 to 62 to estimate the corresponding values of the parameters.
Not only do these samples speed up fitting by providing initial values for any optimization algorithm but they also allow us to visualize both the parameter space and the 'energy' landscape that the fitting algorithm explores (Sec. S2 and [35]). This visualization can indicate, for example, that more data is needed to reduce the size of the relevant region of parameter space (by reducing the error on the fourth factorial cumulant) or that cPCH should be used to resolve the two species by their diffusion too.

The fitting algorithm
In fluorescence fluctuation spectroscopy, the fit energy is often minimized using algorithms for nonlinear least-squares [8,11,28]. In our hands, such algorithms (Levenberg-Marquardt) perform poorly for brightness analysis, including cPCH, probably because the parameter space has multiple local minima.
Given the rugged energy landscape, we minimize the fit energy using nested sampling. Nested sampling is a Markov chain Monte Carlo method where multiple sets of parameters are simultaneously optimized, with the worst solutions being replaced at every iteration by randomly perturbing another solution [40,41].
Briefly (for more details see [35]), we: 1. Sample (typically 100) initial sets of parameters from empirical factorial cumulants using the inverted solutions (Eqs 60-62); 2. Rank each set of parameters by its energy using Eq 58; 3. Record the parameter set with the highest (worst) energy and replace these parameters with a randomly chosen copy of another parameter set; 4. Update the copied parameter set, one parameter at a time, using 20 iterations from an update distribution (either a normal or, for the diffusion coefficient, log normal distribution). If the parameters are still within their bounds, always accept if the new energy is less than E min or otherwise accept using a Metropolis step; 5. Either decrease the step size for this parameter set by the exponential of the reciprocal of the number of accepted steps if the number of accepted steps is greater than the number of rejected steps or increase the step size by the exponential of the reciprocal of the number of rejected steps if the number of rejected steps is greater than the number of accepted steps; 6. Remove and record the parameter set with the next highest energy and update this parameter set using the steps above; 7. Terminate if the desired number of iterations has passed, the energy has decreased to a specified value, or the parameters no longer change beyond a certain threshold.
We use the average of each of the recorded parameters to find the best-fit values (after an appropriate burn-in).
We adjusted the acceptance step to have an E min (typically set to 2) to explore parameter sets that are compatible with the data within the expected error rather than only lying in the region with the lowest energy. Note that because parameter sets can have symmetric solutions (by swapping the parameter sets for species 1 and 2), we enforced one parameter in one species to be greater than or equal to the same parameter in the other.

Results
The shape of the cPCH distribution changes with the concentration, brightness, and diffusivity of the fluorophores and with the time difference, τ, between the two detection channels (Fig 1). Let τ d be the time taken for a fluorophore to diffuse across the observation volume.
Then for small τ/τ d , there is a strong correlation between counts in the τ = 0 channel and counts in the time shifted channel (Fig 1B). Molecules have not had time to diffuse to a region of different intensity, and correspondingly there is a higher probability of identical photon counts in each channel. As the τ/τ d ratio increases, the correlation between n x and n y counts decreases, and the histogram becomes more symmetric.

FCS, moment analysis, and PCH are special cases of cPCH
From the same data, cPCH can recover the results of FCS [7], moment analysis [11], and PCH [8,9]. Similarly, scanning FCS [42,43] and RICS [24,36] can be considered special cases of spatial cPCH. We generate a single point on an FCS curve from the joint first moment of the cPCH distribution, M 1,1 (τ), normalized by the mean of each channel (Eq 39 and Fig 2a), and the PCH distribution by summing over all the counts in either channel of the cPCH distribution (Eq 42 and Fig 2b). We can determine the factorial cumulants obtained using FCA by   Fig 1B). B) Obtaining a PCH distribution: summing over one of the channels of the cPCH distribution gives the PCH (parameters are � a = 0.7, N = 6, τ d = 50 μs, and s = 5, with a time shift of τ = 25 μs). https://doi.org/10.1371/journal.pone.0226063.g002 Correlated photon-counting histograms setting either m or n to zero in the factorial cumulants (K) mn (τ) (Eq 33). Finally, the dual color PCH is equivalent to the dual color cPCH distribution at τ = 0 (Eq 44).

Spatial cPCH
Extending FFS techniques to images enables the analysis of membrane-bound and immobile proteins, and we compare the photon counts from one pixel in the image to those of another in spatial cPCH. For raster-scanned images, these pixels are temporally separated by the time the laser needs to move between the two. The probability of detecting a correlation between photons from the first pixel and photons from the second pixel depends on the fluorescent molecule being able to move from the first pixel to the second in this time. The probability therefore depends on the distance between the two pixels, the size of the pixels, the molecule's diffusion, and the scan rate of the laser (Eq 54).
As in cPCH, the spatial cPCH distributions show higher probabilities along the n x = n y axis at short correlation times (corresponding to short distances between pixels) because there is a high probability that the molecules will occupy a region of similar intensity in the later pixel ( Fig 1C). As the distance between pixels increases, the molecules are no longer likely to be in regions of similar intensity, and the distribution becomes more symmetric (for the single color case). Like scanning FCS [44] and RICS [24,36], spatial cPCH can be used to directly estimate the diffusion coefficient, rather than the diffusion time. Using spatial cPCH, we can also measure the brightness of immobile molecules (S3 Fig) and so distinguish between the image's mobile and immobile fractions.

Validation using a single species
To validate cPCH, we simulated a single fluorescent species diffusing through a Gaussian observation volume. To ensure that interpreting the analysis is unaffected by sampling errors, we simulated a long (14 minute) experiment, but we emphasize that in practice cPCH requires similar amounts of data as other FFS methods. We fit both the cPCH and the factorial cumulants to the simulated data. Similar to a calibration measurement, we also estimate parameters that characterize the observation volume (the shape factors γ 3 and γ 4 from Eq 11 with γ m+n = γ m,n and the structure factor s = z A /r A from Eq 7).
Inference of both the factorial cumulants and the biophysical parameters is accurate ( Table 1). The minimum values of the fit energies are close to 1, which is similar to their values evaluated with the true parameters. We therefore conclude that the approximations behind cPCH adequately describe the behavior of the simulated fluorescent molecules. Correspondingly, the empirical factorial cumulants match those calculated from Eq 33 using the inferred values of the parameters (Fig 3A). As expected, the higher order factorial cumulants have more statistical uncertainty. We note that we do not normalize the factorial cumulants (unlike in higher order correlation analyses [16]) because we find that the un-normalized factorial cumulants have better sensitivity to brightness, allow more cumulants to be fit, and permit the use of analytical solutions to speed-up the fitting algorithm (Methods). A) Fitting the factorial cumulants of the cPCH distribution gives strong agreement between the empirical factorial cumulants estimated from simulated data and those calculated with the best-fit parameters (using Eq 33 and the parameters in Table 1). For the lower moments, the empirical and estimated factorial cumulants overlap. B) Fitting the factorial cumulants of the spatial cPCH distribution gives similar strong agreement (crosses denote the empirical factorial cumulants and the dashed lines show the fits). We consider a mixture of a slowly (D = 0.1 μm 2 /s) and a quickly (D = 200 μm 2 /s) diffusing species. https://doi.org/10.1371/journal.pone.0226063.g003

Correlated photon-counting histograms
To test spatial cPCH, we simulated molecules diffusing in a three dimensional box whose centre was scanned with a simulated laser-scanning microscope. Tests using a single species achieved similar accuracy to the results for regular cPCH ( Table 2). The secondary peaks in the factorial cumulants occur when slowly moving molecules are sampled again on subsequent lines of the image and can be seen in both the theoretical and simulation results (Fig 3B).

Resolving different species using cPCH
cPCH is better able to resolve two fluorescent species than both PCH and FCA because cPCH uses information on both diffusion and brightness. To determine the cPCH's ability to resolve, we simulated a mixture of two different species with the brightness of one species three times that of the other. When the two species have different diffusion coefficients and different brightnesses, not only does cPCH outperform PCH and FCS (Fig 4A) but also allows more accurate inference of all parameters. We see a similar improvement in accuracy for a mixture of fluorophores with one species twice as bright as the other (S5 Fig). When the two species have close diffusion coefficients but different brightnesses, cPCH performs as well as PCH and FCA (Fig 4B), although also allowing inference of the diffusion times.
cPCH performs at least as well as PCH and FCA in other tests too. We simulated a mixture of two different species with the brightness of one species double that of the other. In this case, both fluorophores are relatively dim, but the brighter fluorophore outnumbers the dimmer fluorophore, a scenario that is challenging for FFS methods (S6 Fig). With cPCH, we could resolve the concentration of the dimmer species, but with less precision than before. Nevertheless, cPCH outperforms PCH, FCA, and FCS (S6 Fig). We also determined the ability of spatial cPCH to resolve two species (Fig 3B and S4 Fig), with similar performance.

Binning effects
It is sometimes necessary to use larger binning times to improve the signal to noise ratio when the molecular brightness is low. By using larger binning times, however, the individual molecules have time to diffuse to locations of different intensity within a bin time, and so diffusion effects must be included in the analysis.
In cPCH, we apply the FIMDA binning approximation (Eq 48) to higher order binning functions (Eqs 49 and 50). This approximation works well up to the third order binning function (Fig 5A). The greater statistical uncertainty in the measurements of the higher order moments (we consider up to fourth order) mean that these moments receive less weight when Table 2. Spatial cPCH can achieve similar accuracy as cPCH. The simulation has a duration of 25 frames and uses an image size of 256 × 256 with a pixel size of δ x = δ y = 11.7 nm, a pixel dwell time τ p of 10 μs, a line retracing time τ L of 1 ms, and a frame retracing time of 5 ms.

Parameter
Actual Correlated photon-counting histograms fitting the factorial cumulants in cPCH, and the errors caused by the approximation have little impact.
In simulations, including the binning functions improves both empirical estimates of the factorial moments ( Fig 5B) and inference of biophysical parameters (Fig 5C), even when T < τ d .

Resolving mixtures of species with multiple binding sites
Finally, we demonstrate that cPCH can work in a scenario where FCS, PCH, and FCA are all unable to resolve the different species. Consider a receptor with two different binding sites, each with its own affinity for a fluorescent ligand, and with the receptor diffusing slower than the ligand. There are three different fluorescent species-unbound ligand (denoted L), receptors with a single bound ligand (R 1 ), and receptors with two bound ligands (R 11 )-and one non-fluorescent species-unbound receptors (R). FCS cannot differentiate between the singlyand doubly-bound receptors (R 1 and R 11 ) because their diffusion coefficients do not differ (being determined by the more massive receptor); PCH and FCA cannot differentiate between the unbound ligand and the singly-bound receptors (L and R 1 ) because both have the same fluorescence. We simulated this system at equilibrium using an initial amount of ligand L o and analyzed the data using spatial cPCH. From the empirical histograms of spatial cPCH, we calculated and then fit the factorial cumulants using a model with three fluorescent species (L, R 1 , and R 11 ; Sec. S5). The concentration of R can be estimated if R 0 , the initial concentration of free receptor, is known, via R 0 = R + R 1 + R 11 . Spatial cPCH is able to accurately estimate all four concentrations, provided we use concentrations of ligand below the concentration at which R 1 , the singly bound receptors, reaches its maximum (Fig 6A). At higher concentrations of ligand, the doubly-bound receptor, R 11 , and free ligand, L, dominate the signal, and R 1 is difficult to estimate. Errors in R 1 propagate through to errors in the estimates of R, and both errors affect estimates of the dissociation constant K � 1 ¼ LR R 1 (Fig 6B). An alternative is to directly estimating R using a dual-color assay. We note that we could also use cPCH (rather than spatial cPCH) because both receptors and Correlated photon-counting histograms ligand are freely diffusing. If receptors are confined to a membrane, spatial cPCH is the preferred technique.
If either FCS or PCH are used instead of cPCH, estimates of the different concentrations would have to be deduced from the overall behaviour of the fluctuations in fluorescence as a function of concentration, making the task more difficult [45]. The initial concentration of receptors is R 0 = 10 nM. We specified the (apparent) dissociation constants for the singly and doubly bound receptors as K � 1 ¼ 2 nM and K � 2 ¼ 12 nM and find the actual dissociation constants from Eq S9. These apparent dissociation constants arise because we cannot distinguish using either fluorescence or diffusion which site the ligand has bound on a receptor that has bound only one ligand (Sec. S5). The simulation employed a laser dwell time of 10 μs with images of 512 × 512 pixels and an image size of 4 μm by 4 μm, so that the pixel size is just under 12 nm. The line retracing time is 1 ms and the frame retracing time is 50 ms. We generated 2 sets of data with 25 images each for all of the concentrations used, except for the lowest initial value of ligand, which used 1 set of 70 images. B) The apparent dissociation constants (K � 1 ¼ LR=R 1 ¼ 2 nM and K � 2 ¼ LR=R 11 ¼ 12 nM) can be estimated using the measured concentrations provided the concentration of singly labeled receptors R 1 is greater than the concentration of doubly labeled receptors R 11 . https://doi.org/10.1371/journal.pone.0226063.g006

Discussion
We have introduced cPCH, which unifies several FFS techniques (FCS, PCH, and FCA) to resolve different molecular species through both their brightness and diffusion. PCH and dual-color PCH are special cases of the cPCH distribution; FCS and FCA are special cases of the moments of the cPCH distribution.
Although data can be fit to the cPCH distribution, we usually find fitting its factorial cumulants to be more computationally efficient. Using simulated data, we have shown that the factorial cumulants can resolve multiple species at higher concentrations than analyses using the distributions [35]. Nevertheless, there are situations where the distribution can be useful, such as correcting for dead-time and afterpulsing (Sec. S7, S8 and [35]) and determining the shape parameters γ 1 and γ 2 [31]. Further, it is straightforward to sample from the discrete probability distributions providing a convenient method to generate empirical distributions of the photon counts at a given τ and T for testing statistical properties of the cPCH distributions, numerical implementations, etc.
In addition, we have developed a novel fitting algorithm based on a Markov chain Monte Carlo scheme for nested sampling. To improve the algorithm's efficiency, we select initial parameters by inverting and sampling from the analytical solutions to the factorial cumulants (using theoretical expressions for the variance of each moment for the sampling). With this method, we can visualize the parameter space and initialize optimization algorithms in the relevant regions of parameter space.
We show that cPCH resolves two differently diffusing species. cPCH reduces the uncertainty in fitting if the two species have different diffusion and brightness and performs as well as both PCH and moment analysis when the two species differ only in their brightness. For example, resolving two different species by brightness alone is challenging because fourth order moments are required, which are difficult to estimate. cPCH circumvents this problem by also being sensitive to differences in diffusion for which only second-order moments are needed. As PCH and FCS are special cases of cPCH, we can also always choose to use one of these techniques instead of cPCH without having to reformat the data.
Our work is related to that of Melnykov et al. [16]. Their expressions for the higher-order factorial cumulants, although derived differently, correspond to the factorial cumulants of the cPCH distribution. Nevertheless, they normalize these bivariate expressions by the univariate factorial cumulants to allow better comparison with FCS. Our un-normalized factorial cumulants allow us both to more efficiently explore parameter space by exploiting the analytical solutions for the univariate factorial cumulants during fitting and to constrain the diffusion time if we apply a series of bin sizes up to those large relative to the diffusion time (following FIMDA [46] and TIFCA [18]).
Abdollah-Nia et al. have recently employed normalized higher-order FCS to study reaction kinetics [47,48]. They show that by normalizing the higher order factorial cumulants of the reacting system with those of a non-reacting system, the effects of diffusion and the illumination profile can be removed, isolating the reaction kinetics. This method is potentially powerful for studying reaction kinetics, but for species with similar diffusion. We have focused on resolving different species based on differences in their brightness and diffusion and have demonstrated cPCH with a ligand-binding system with two distinct binding sites, where the ligand and receptor diffuse differently. Such a system would be unamenable to the analysis of Abdollah-Nia et al.
Another recent development is 2D fluorescence lifetime correlation spectroscopy (2DFLCS), which differentiates molecules based on their diffusion and fluorescence lifetimes and can analyse molecular conformations on the microsecond scale [49][50][51]. Although similar to cPCH because 2DFLCS uses histograms of the fluorescence lifetimes of photons detected at a time τ apart, it, unlike cPCH, requires a setup for time-correlated single-photon counting. Indeed, cPCH is able to differentiate molecules based on their brightness even if their fluorescence lifetimes are identical, as expected, for example, for fluorescent oligomers. As Fluorescence Intensity and Lifetime Distribution Analysis (FILDA) combines fluorescence lifetime analysis with photon counting histograms [52], we expect that cPCH could be combined with 2DFLCS to allow molecules to be differentiated by fluorescence lifetime, brightness, and diffusion.
With a simple modification, cPCH can be extended to the study of images. Spatial cPCH, like RICS [24] and scanning FCS [44], but unlike FCS and cPCH, is able to directly estimate the diffusion coefficient of the fluorophores (rather than their diffusion time, which depends on both the diffusion coefficient and the beam waist parameter r A ). Knowing diffusion coefficients is useful for calibrating optical systems and for studying fluorescent dyes and proteins.
Although we have focused on raster scanning, the derivation of spatial cPCH would apply too to scanning FCS and directed flow. By scanning the laser over a range of positions, spatial cPCH is also able to sample more molecules in a shorter period of time. Spatial cPCH can therefore be used to quantify the concentrations of immobile fluorophores-a challenge for FCS and cPCH because such fluorophores quickly photobleach with a stationary illumination profile and only produce photon counts that are Poisson-distributed. Future work could include the thickness profile of cells by using z-scanning [53].
Both cPCH and spatial cPCH require photon-counting detectors. Spatial intensity distribution analysis (SpIDA) [25] is an extension of brightness analysis to images, but employs a calibration method to allow the use of analog photomultiplier tubes. Spatial cPCH can also be extended to analog detectors [35].
We have used simulations to verify cPCH and spatial cPCH, but, in practice, detectors often introduce artefacts because of dead-time and afterpulsing, and experiments can be undermined by triplet states and blinking of fluorophores. The theory of cPCH and spatial cPCH can include such effects (Secs. S6-8). Although we assume a 3D Gaussian observation volume (a similar assumption was made in TIFCA [18] and higher-order FCS [16] and adequately describes diffusion but not the shape parameters of the observation volume), our derivation is general, and it is possible to use the shape parameters of a Gaussian-Lorentzian observation volume for γ m,n and use numerical correlation functions for G m,n (τ, τ d ). In our simulated calibration, we optimized the shape parameters γ 3 and γ 4 . cPCH is valid for both single-channel and dual-color analysis. The combined use of two colour channels and information on diffusion can substantially increase accuracy in resolving multiple species by further reducing the requirement for the third and fourth order moments and by causing the joint factorial cumulants to be determined by only interacting species (labeled with a fluorophore from each channel). We therefore expect dual-color cPCH to be the natural choice for studies that consider interactions between labelled molecular species.
In summary, we anticipate that cPCH and spatial cPCH will become powerful techniques for quantitative, single-cell biology. Simulated calibration data used to estimate p(τ) and a demonstration of the removal of afterpulsing effects between time bins for an afterpulse probability q � = 0.01 and an average arrival time τ ap = 200 ns. A) Histograms of arrival times and afterpulses. B) Fit to the measured p(τ). C) A comparison of the calculated FCS curves for a 1 nM solution with a diffusion coefficient of 1 × 10 −10 m 2/ s, with afterpulsing, without afterpulsing, and with afterpulsing removed (using simulated data). The black dotted lines show the FCS curves of the original data without afterpulsing; the green curves show the FCS curves with afterpulsing added to the data; the dashed red curves show the afterpulse-corrected FCS curve using the measured p(τ); and the blue curves show the afterpulse-corrected FCS curve using the ideal p(τ). (EPS) S1 Table. Summary of the simulated calibration results including dead-time effects, triplet states, and binning effects. Both triplet and diffusive binning are necessary. In the third column, the minimum bin size was T = 200 ns and the maximum was T = 3.2 μs. In the fourth column, the minimum bin size was T = 1 μs and the maximum was T = 64 μs. � Both diffusion and triplet binning included in the fit. �� Only diffusion binning included in the fit. R denotes the triplet fraction; τ t denotes the triplet relaxation time. (PDF)