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

Summary

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.

Keywords:

Probability distribution – Lasers – Mass diffusivity – Fluorescence – Photons – Dwell time – Fluorescence spectroscopy – Generating functions

Introduction

Absolute numbers are necessary for quantitative biology, but are challenging to obtain in single cells. Although ad hoc techniques exist [1–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 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–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–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.

Theoretical foundations

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₁(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:

The term p(x, y|τ) in Eq 1 is an important difference between this equation and the single-molecule 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

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 dual-color PCH. With the diffusion term, the integral gives a higher-order two-point correlation function, G_m,n(τ):

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]:

Independently treating each spatial dimension in Eq 5 gives

The shape factors of the observation volume are defined as

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 (γ 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:

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 x ′ and n y ′ emitted by each molecule). Therefore we have

To proceed, we use the bivariate probability generating function given by:

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′ and e^is′ [30]) and directly calculate the cPCH distributions using the inverse Fourier transform, F -⁠ 1:

Recursive solution

The probability generating function can also be used to derive the probabilities. By definition [30]

• 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:

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₁(t, s|τ) is

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:

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

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⁻¹ with N_d being the number of data points (typically on the order of 10⁶), we can write the variance of either the cumulant or factorial cumulant of interest, X_p,r, as

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:

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

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 | τ ) = ∑ n y p N ( n x , n y | τ ).

For example, using the single-molecule probability distribution p₁(n_x, n_y|τ) (Eq 3 in Eq 1):

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₁(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₁(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:

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₂(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:

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.,

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

cPCH

Deriving a general equation for B_m,n(T, τ_d) is challenging, and we instead use FIMDA’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₁₁(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:

Using Eqs 31 and 9, we can now define the factorial moments in Eq 32 as

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 raster-scanned 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

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 line-retracing 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 τ_pd_x; the time for the laser to reach the position d_y by scanning the previous lines is τ_pL_yd_y; and τ_Ld_y is the corresponding line-retracing time. Consequently, τ in Eq 55 becomes τ_p(d_x + L_yd_y) + τ_Ld_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.

Methods

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 (γ₃, γ₄).

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 | τ , T ), can then be converted to the empirical probability distribution, p N e m ( n a , n b | τ , 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 k s. 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 k s = 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 τ + as b τ , i + = b i + τ k +.

To determine p N * ( n a , n b | τ , T ), we first find p N * ( n a > 0 , n b > 0 | τ , T ) by constructing a histogram for all pairs of bins with non-zero photon counts from the intersection of b⁺ and b τ +, with the first element of the pair, n_a, from b⁺ and the second element, n_b, from b τ +. 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 τ + do not intersect: for example bins that are in b⁺ but not in b τ + 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 | τ , 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

Writing L n x , n y ( τ ) as the total number of histogram entries observed at a given τ, we define the fit energy—or loss function—for the cPCH distribution as

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 ( τ ) = ∑ n x ∑ n y n x m n y n p N em ( n x , n y | τ ). 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, σ K 2 ( τ ), 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

Using analytical solutions for the factorial cumulants to sample parameter space

To resolve two species by brightness alone using factorial cumulants, we need to find {N₁, ϵ₁, N₂, and ϵ₂} that solve, within experimental error:

Eq 59 can be inverted. We find that

There is a single, unique solution because the two solutions are symmetric. Statistical errors in the empirical factorial cumulants can result in non-physical (negative) values for these analytical solutions.

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:

• Sample (typically 100) initial sets of parameters from empirical factorial cumulants using the inverted solutions (Eqs 60–62);

• Rank each set of parameters by its energy using Eq 58;

• Record the parameter set with the highest (worst) energy and replace these parameters with a randomly chosen copy of another parameter set;

• 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;

• 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;

• Remove and record the parameter set with the next highest energy and update this parameter set using the steps above;

• 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 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 γ₃ and γ₄ 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).

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 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₁), and receptors with two bound ligands (R₁₁)—and one non-fluorescent species—unbound receptors (R). FCS cannot differentiate between the singly -⁠ and doubly-bound receptors (R₁ and R₁₁) 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₁) 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₁, and R₁₁; Sec. S5). The concentration of R can be estimated if R₀, the initial concentration of free receptor, is known, via R₀ = R + R₁ + R₁₁.

Spatial cPCH is able to accurately estimate all four concentrations, provided we use concentrations of ligand below the concentration at which R₁, the singly bound receptors, reaches its maximum (Fig 6A). At higher concentrations of ligand, the doubly-bound receptor, R₁₁, and free ligand, L, dominate the signal, and R₁ is difficult to estimate. Errors in R₁ propagate through to errors in the estimates of R, and both errors affect estimates of the dissociation constant K 1 * = L R 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 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].

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 γ₁ and γ₂ [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–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 γ₃ and γ₄.

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.

Supporting information

S1 Fig [eps]
There are many solutions that have energies low enough to be considered within the experimental uncertainty of the data.

S2 Fig [eps]
By measuring ϵ , there is a minimum in the energy as a function of the parameter values.

S3 Fig [a]
As immobile and mobile fluorophores show different decays in their joint factorial cumulants, spatial cPCH can distinguish one from the other.

S4 Fig [eps]
Spatial cPCH can resolve two species with different diffusion coefficients.

S5 Fig [eps]
cPCH can outperform alternative methods: A mixture of fluorophores, with one species twice as bright as the other.

S6 Fig [eps]
cPCH can out perform alternative methods: Two dim species.

S7 Fig [eps]
Dead-time corrections are necessary when dead-time is present.

S8 Fig [eps]
Correcting for afterpulses.

S1 Table [pdf]
Summary of the simulated calibration results including dead-time effects, triplet states, and binning effects.