This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.105.066084v1 90/6/2179    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Palo, K. Articles by Kask, P. Search for Related Content PubMed PubMed Citation Articles by Palo, K. Articles by Kask, P.

Calculation of Photon-Count Number Distributions via Master Equations

Evotec Technologies, Harku, Estonia

Correspondence: Address reprint requests to Dr. Peet Kask, Evotec Technologies GmbH, Instituudi tee 11, 76902 Harku, Estonia. Tel.: 37-2-651-1146; Fax: 372-651-1146; E-mail: peet.kask{at}evotec-technologies.com.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Fitting of photon-count number histograms is a way of analysis of fluorescence intensity fluctuations, a successor to fluorescence correlation spectroscopy. First versions of the theory for calculating photon-count number distributions have assumed constant emission intensity by a molecule during a counting time interval. For a long time a question has remained unanswered: to what extent is this assumption violated in experiments? Here we present a theory of photon-count number distributions that takes account of intensity fluctuations during a counting time interval. Theoretical count-number distributions are calculated via a numerical solution of Master equations (ME), which is a set of differential equations describing diffusion, singlet-triplet transitions, and photon emission. Detector afterpulsing and dead-time corrections are also included. The ME-theory is tested by fitting a series of photon-count number histograms corresponding to different lengths of the counting time interval. Compared to the first version of fluorescence intensity multiple distribution analysis theory introduced in 2000, the fit quality is significantly improved. It is discussed how a theory of photon-count number distributions, which assumes constant emission intensity during a counting time interval, may also yield a good fit quality. We argue that the spatial brightness distribution used in calculations of the fit curve is not the true spatial brightness distribution. Instead, a number of dynamic processes, which cause fluorescence intensity fluctuations, are indirectly taken into account via the profile adjustment parameters.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
In life science, fluorescence is a widely used reporter of the state and dynamics of the studied sample. It is an extremely sensitive reporter: even single fluorescent molecules can be detected and identified. There is a family of fluorescence methods that monitor and make use of fluorescence intensity fluctuations emitted from a microscopic volume containing a low number of fluorescent molecules on average. The name of this family is fluorescence fluctuation spectroscopy (FFS). Historically, the first FFS method is fluorescence correlation spectroscopy (FCS) (1). In FCS, the autocorrelation function of fluctuating fluorescence intensity is calculated and analyzed. It is most widely used to study diffusion. For applications where two or more species have to be resolved, like in drug screening, a few other FFS methods have been developed that make use of molecular parameters other than diffusion time.

Eighteen years after the first realization of FCS, photon-count number histogram was introduced as an alternative to correlation function for fluorescence fluctuation analysis (2), rendering fluorescence brightness as a molecular parameter. Nine years later, two independent research groups reported of experiments on successful fitting of photon-count number histograms (3,4). They have used different names for the method, photon-counting histograms (PCH) and fluorescence intensity distribution analysis (FIDA), respectively. The two approaches differ in the assumed shape of the spatial brightness distribution (which influences the shape of the theoretical count-number distribution) and in details of the computation algorithm (which may affect the speed of analysis but not the results). However, the two approaches are similar in sharing the main physical assumption: that fluorescence intensity emitted by a molecule is constant during a counting time interval.

FIDA was developed to be applied in high-throughput drug screening under a demand of short data-acquisition time (<2 s) and online analysis. It has indeed been used successfully for that purpose for years (5–7). In particular, it is most appropriate in quantifying biological assays where a single particle (such as vesicular particle) carries a great number of binding sites for a labeled ligand. In this case, the two fluorescent species (receptor carrier and unbound ligand) have a huge brightness contrast and are thus easily resolved. Both PCH and FIDA are thus ideal tools for evaluating the state of chemical equilibrium in such an assay type (6,8–10). In cases of a modest brightness contrast, it makes sense to use two parallel detection channels that monitor different polarization or spectral components of fluorescence. Fitting a joint histogram of photon-count numbers originating from two different detectors (two-dimensional FIDA) (11) has turned out to be an efficient analysis method (5,6,12–14). Two other FIDA-based methods have been reported later that distinguish molecules according to several molecular characteristics: fluorescence intensity multiple distribution analysis (FIMDA) according to diffusion time and brightness (15), and fluorescence intensity and lifetime distribution analysis (FILDA) according to brightness and lifetime (16). Two-dimensional FIDA and FIMDA have been established on the basis of the PCH algorithm, too (17,18). All the above-mentioned methods constitute a family and share the theory discussed in this article.

It is known from FCS theory and experience that fluorescence intensity of a molecule fluctuates significantly within a counting time interval of PCH or FIDA (which is typically 10–40 µs). First of all, the fluctuations occur due to diffusion into and out of the monitored volume, and due to transitions between the singlet (visible) and triplet (dark) electronic states. These two processes are routinely taken into account in FCS. Starting in 2000, the assumption of constant emission intensity of a molecule during a counting time interval has been critically inspected. First, the second factorial cumulant of photon-count number distribution could be calculated as a function of the width of the counting time interval (15). A better accuracy than that provided by the second cumulant has been a want for PCH/FIDA and related methods, but the mathematics for this have been absent. To our knowledge, Monte Carlo method that has been applied in photon arrival-time interval distribution (PAID) (19) is the only successful approach in this respect so far. This calculation method suffers from long calculation times and a cosmetic problem: the random errors characteristic to Monte Carlo approach.

This work was undertaken to investigate the effect of intensity fluctuations during the counting time interval on the shape of the count-number distribution. As the outcome, we have succeeded in calculating the theoretical count-number distribution as a function of the width of the counting time interval, concentration, count rate per molecule, three profile parameters, diffusion time, triplet state lifetime and population, detector afterpulsing probability, and detector dead-time. We shall demonstrate that many factors, in addition to the brightness profile, significantly influence the count-number distributions under conditions typical to FFS experiment. It will be explained why and how a good fit quality in count-number distributions can be nevertheless achieved when using a simple theory that apparently ignores intensity fluctuations during a counting time interval.

	THEORY

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Comparison of PCH and FIDA
This section compares the theory of PCH (3) and its later modification (20) with that of FIDA (4). For simplification, we shall use a single terminology. A common aspect of all above-mentioned theories is that they all make use of the assumptions that:

1. Emission intensity of a molecule is a function of its translational coordinates only.
   
2. Intensity fluctuations during a counting time interval are neglected.
   

(1)

with In Eq. 1, V is the size of the volume where the molecule is enclosed, q is count rate per molecule in the focus, T is width of the counting time interval, and B(r) is the spatial brightness distribution (known also as the observation volume profile, or point spread function). Here and below throughout this article we use the convention that the spatial brightness is unity in the exact focus. (This convention is useful for theoretical calculations and their interpretation, but unusual in applications.) Eq. 1 is a special case of Mandel's formula (21), with as a probability that the molecule is in the given volume element. Poisson distribution alone expresses the probability distribution of the count number under a constant light intensity.

The description of the shape of the spatial brightness function (to be referred to in this article as the "brightness profile") is in fact the only part where the reviewed two theory versions significantly differ. The other steps involved may differ much in technical details and computational algorithms. However, applying a common brightness profile would yield identical count-number distributions.

In connection with two-photon excitation, the square of Gaussian-Lorentzian profile for B(r) is the logical selection (3). It exactly matches the squared intensity profile of an ideal focused TEM₀₀ laser beam,

(2)

Here, is the distance from the optical axis, z is the distance from the focal plane, w₀ is the waist radius of the brightness profile (here defined as the distance on the focal plane from the optical axis where the emission intensity has decreased by a factor of e²), and z₀ is the distance from the focal plane where the spot radius has increased by a factor of

A remarkable property of Eq. 2 is that under the assumptions 1 and 2, above, neither of the parameters w₀ and z₀ is suited to serve as a profile adjustment parameter that would affect the shape of the count-number distribution. A single-molecule count-number distribution has been analytically expressed for the brightness profile given by Eq. 2 (see Eq. 15 of (3)). Profile parameters contained in that formula simply represent a scaling factor. However, the shape of the count-number distribution of a single molecule only depends on molecular brightness. This conclusion can be also drawn by studying the unitless characteristics of the profile (such as where µ-values denote profile moments). They are independent of the parameters w₀ and z₀.

It is likely that this ideal brightness profile is closer to the real one in the case of two-photon excitation compared to the one-photon excitation. Indeed, in the case of one-photon excitation, a detection pinhole is a necessary element of the optical setup. In the context of photon-count number histograms, it serves as an additional source of spatial brightness distortions. Contrary to the study establishing PCH (3), one-photon excitation has been used in the study establishing FIDA (4). For calculations of the distribution of the number of photon-counts from a molecule, the three-dimensional spatial brightness distribution can be represented by its reduced form, a one-dimensional function expressing the relationship between volume and brightness. In reference (4), an empirical relationship between volume and the logarithm of inverse brightness has been expressed as

(3)

Here, a₁ and a₂ are shape adjustment parameters that are determined by fitting a photon-count number histogram in a calibration experiment. A₀ is a scale factor related to the size of the observation volume.

When developing the mathematical model, fit quality of count-number distributions was considered even more important than physical essence. A physically reasonable model of the brightness profile would rather look like Here, the first term would exactly describe a Gaussian profile and the terms with the adjustment parameters a₁ and a₂ would be considered corrections. However, Eq. 3 was selected because it yields a better fit to count-number histograms. Later we shall return to the question whether and why this was a good choice.

An attempt to apply PCH algorithm with an ideal (three-dimensional Gaussian) brightness profile turned out to yield unsatisfactory fit quality (20). In that work, most of the attention is concentrated on finding what the true brightness profile looks like. Electromagnetic diffraction theory was used that is a rather time-expensive means to calculate the three-dimensional profile of the equipment. Still, using the calculated brightness profile, count-number histograms could not be fitted well. However, when studying moments of the calculated profile µ_j, it was realized that products j^3/2µ_j are approximately constant for j > 1. Theoretically, such products are constants (independent of j) for the pure three-dimensional Gaussian profile. Thus, it was concluded that the shape of the calculated profile might be represented by a single parameter describing how much the first moment deviates from the ideal value extrapolated from higher moments. This parameter was used as a floating parameter when fitting count-number distributions. Fit quality was found to be good. Fit quality was further improved when the second moment of the profile was also adjusted by another floating parameter. Thus, the modern version of PCH has introduced profile adjustment parameters, like the original FIDA, only on a different mathematical basis. To what extent the adjustment parameters describe the true brightness profile and to what extent they express dynamic processes during the counting time interval is a question to be answered below.

Theories distinguishing between static and dynamic brightness
A step beyond the borders defined by assumptions 1 and 2, above, was taken when introducing a method called fluorescence intensity multiple distribution analysis (FIMDA) (15). In FIMDA, a series of theoretical distributions are fitted to a corresponding series of photon-count number histograms collected simultaneously at different widths of the counting time interval. In this way, FIMDA resolves fluorescent species based on both molecular brightness and diffusion time. A key equation in this theory relates the second factorial cumulant of photon-count number distribution with the autocorrelation function of fluorescence intensity G(t) = I(0)I(t) – I²,

(4)

T denotes the width of the counting time interval, I is the detected fluorescence intensity at time t, and ... denote ensemble average. This equation can be applied to each species separately and is interpreted in terms of the dependence of the apparent particle number and the apparent count rate per molecule q(T) = K₂(T)/TK₁(T) on sampling time T, where K_j(T) denotes the j^th factorial cumulant of the count-number distribution. An expression for G(t) was taken from the FCS theory, where it has been derived for a three-dimensional Gaussian brightness profile (22,23). Brightness of a molecule if averaged over a counting time interval is a qualitatively different variable than the static spatial brightness. The form of the distribution of the dynamic brightness was approximated by a modification of Eq. 3, i.e., by a flexible empirical formula that was originally designed for FIDA to approximate a non-Gaussian spatial distribution of unknown shape. A single brightness distribution was applied to histograms of different widths of the counting time interval.

Another highly sophisticated theory is covered by photon arrival-time interval distribution (PAID) (19). PAID theory is based on simulation of a high number of diffusion pathways of a single molecule and the conversion of these random pathways into fluctuating intensities as functions of time. Since the number of random paths is finite, the calculated distributions of photon intervals or count numbers have a random noise. Nevertheless, the algorithm can even be used for fitting experimental data, as long as the same trajectories are used each time. This grants reproducibility of theoretical curves. Scaling of time is a means to recalculate photon statistics at a different diffusion coefficient. This theory can easily be modified to take account of other effects such as triplet-state induced fluctuations.

Master equations
From this section on, we will describe a new approach that is free of the cosmetic problem of randomness. We will study a molecule in a given closed volume, undergoing diffusion, singlet-triplet transitions, and emitting photons that can be detected by an optical setup. P_S(n,T,r) denotes the probability that the studied molecule is in the singlet state, is located in point r at the end of the counting interval, and has emitted n photons. P_T(n,T,r) expresses a similar probability, but with the studied molecule being in the triplet state at the end of the counting time interval. In our model, we do not distinguish between the ground and the excited state; in the singlet state, the molecule is visible, whereas in the triplet state the molecule is dark. The probability density of the triplet-to-singlet transition is constant (inverse lifetime of the triplet state), whereas the probability density of singlet-to-triplet transitions is proportional to laser intensity. Our starting point is the Master equations, which is a system of differential equations describing the evolution of the molecule:

(5)

The first term on the right side describes diffusion, where D is diffusion coefficient and is the Laplace operator. The second term describes triplet-to-singlet transitions, where _T is the lifetime of the triplet state. The third one stands for singlet-to-triplet transitions, where is the probability density of singlet-to-triplet transitions exactly in the focus. The fourth term in the first equation describes photon detection. In Eq. 5, we have distinguished the excitation E(r) and the transmission F(r) profiles.

We call the theory outlined here and below as the Master equation (ME) theory, as it is based on Master equations.

Statistics of occupation times for a two-state model
In the first approach, the molecule under study is assumed immobile. Thus, it is located in a point in space where intensity of the exciting light is constant. We assume that the molecule emits light of a constant intensity when it is in the singlet electronic state, and is invisible when in the triplet state. Thus, the brightness integrated over a counting time interval is proportional to the time that the molecule spends in the singlet state. Our problem is reduced to calculation of the distribution of time spent in the singlet state.

In this section, we study a model system with two discrete states (1 and 2) and random transitions between them with given rate constants k₁₂ and k₂₁. We first assume that the system is initially in State 1; we are interested in statistics of overall occupation times t₁ and t₂ that the system has spent in States 1 and 2 out of the full time T = t₁ + t₂. It is the easiest to express the probability of no transitions during time T:

(6)

In cases when transitions occur, the distribution of time t₁ is described by probability density functions. Let us first derive an expression for the probability density function for occupation time t₁ corresponding to a given number (n) of transitions during time interval T, W(t₁;n,T). For example, with two transitions, there are three time intervals: from time zero until the first transition t_1a, from the first until the second transition t₂, and from the second until the end, t_1b. There are two independent variables out of the three; therefore we deal with the probability density of two independent variables (e.g., t_1a and t_1b) which can be expressed as a product of the following factors: the probability that no transition occurs during t_1a; the probability density of the first transition; the probability that no backward transition occurs during t₂; the probability density of the second transition; and the probability that no forward transition occurs during t_1b:

(7)

It remains to integrate this expression over t_1a (or t_1b) that yields a probability density of a single variable,

(8)

In the same way, one can derive expressions for any number of transitions. In fact, we distinguish two cases, with odd and even number of transitions, since they yield two different final states of the system. For an arbitrary odd number of transitions (n = 2s + 1), the probability density function is

(9)

For an arbitrary positive even number of transitions, the corresponding expression is

(10)

Now it remains to sum up contributions from all odd and from all even numbers of transitions separately. One should recognize a power-series of Bessel functions, yielding

(11)

(12)

Note that I₀ and I₁ denote Bessel functions.

Naturally, one may interchange indexes 1 and 2 in the expressions of this section, getting the solution for the case when the system is initially in state 2.

Distribution of the time that a molecule spends in the singlet electronic state
Now we may use Eqs. 6, 11, and 12 to express the distribution of the time that an immobile molecule spends in the singlet state. As initial conditions, we apply the steady-state occupation probabilities, p_S = k_TS/(k_TS + k_ST) and p_T = k_ST/(k_TS + k_ST). The analytic solution is

(13)

Fig. 1 illustrates the distribution of the fraction of time spent in the singlet state for k_ST = 0.125 x 10⁶ s^–1, k_TS = 0.5 x 10⁶ s^–1, and T = 10 µs. This set of parameters models a tethramethylrhodamine molecule with fixed coordinates monitored by a fluorescence microscope. As a closer to an experiment case, on Fig. 2, a distribution of the dynamic brightness is graphed that involves not only a single but all possible coordinates of the (immobile) molecule. One may conclude from visual inspection of the graphs that singlet-triplet transitions are sufficient to cause significant distortions in the shape of the count-number distribution, compared to the corresponding pure profile-defined distribution.

View larger version (10K): [in this window] [in a new window]   FIGURE 1  An exemplary distribution of the relative time that a molecule spends in the singlet state. This graph is calculated for kST = 0.125 x 106 s–1, kTS = 0.5 x 106 s–1, and T = 10–5 s. The distribution consists of a continuous function plus two -functions. The left -function covers an area of only 0.0013, which is the probability that the molecule is in the triplet state throughout the entire counting time interval. The right -function covers an area of 0.229, which is a probability that the molecule is in the singlet state all the time. Overall distribution is normalized to unity.

View larger version (14K): [in this window] [in a new window]   FIGURE 2  Distribution of the average brightness in a counting time interval accounting all possible coordinates of the molecule. Here both the fraction of time that the molecule spends in the singlet state and spatial coordinates of the molecule are random variables that determine the average brightness, but the molecule is assumed to be immobile during the counting time interval. The Y axes measure relative contribution to emitted intensity. The distribution has been graphed for a case without singlet-triplet transitions and for cases with singlet-triplet transitions at two different widths of the counting time interval. The area under the first graph is unity, while the two other graphs correspond to fluorescence of 6.5% lower intensity. The three distributions are similar at low brightness values having a singularity at zero brightness, but we have selected a scale for Y axes that best demonstrates the difference between the three graphs.

Another special case: pure diffusion
Another special case of Eq. 5 is when the terms describing singlet-triplet transitions are neglected. In this case, our search for analytical solutions has been less successful, compared to the content of the previous section. We have succeeded in numerical calculation of a limited number of moments of the desired function. Details of the calculations will not be presented here because the results are not used in the rest of our studies. However, as an indirect outcome of pure diffusion studies, we stopped attempts to solve Eq. 5 analytically in the general case and concentrated our attention on numerical algorithms.

Principles of the numerical algorithm
In the original model described by Eq. 5, spatial coordinates and time are continuous variables. For numerical solution of Master equations, we replace this model by another one that is better suited for computer calculations. We assume that the molecule under study may have only a finite number of locations that may be changed by random jumps between neighbor locations in short time-steps of typically t = 0.1 µs. Note that the time-steps of numerical calculations are very short compared to the width of the counting time interval T. Instead of a straightforward Euclidean three-dimensional spatial grid, we select a grid in two-dimensional cylinder coordinates. This reduces the number of grid points tremendously. Furthermore, the density of the grid need not be uniform, and it makes sense to select the highest density in the focus. We have used a constant grid step for the first few grid points, and a geometrically progressing grid step for representing high-numeric coordinates (i.e., out-of-focus periphery). Parameters of the grid (such as the initial step divided by the waist radius of the beam, and the progressing factor) have been selected as a tradeoff between calculation accuracy and speed. The results of calculations must not significantly depend on grid parameters. We have used 36–40 grid points per dimension, making 1300–1600 grid points altogether.

We keep a careful track of the evolution of the probability distribution of discrete variables, P(i,i_z,e,n). Here i and i_z are grid coordinates, e is electronic state, and n is the number of photons already emitted by the molecule. We start from time zero when no photons have been emitted. The molecule is located in each grid point with the probability that is proportional to the volume that the particular grid point represents. We recalculate the distribution in many time-steps, until the end of the counting time interval is reached.

Spatial brightness distribution
For our calculations, we have selected the spatial brightness distribution of

(14)

The second and the third factors on the right-hand side describe the ideal intensity profile of a focused laser beam, while the first factor is a Lorentzian approximation of the transmission function of a microscope with a pinhole. The first term of Eq. 14 resembles the square of Gaussian-Lorentzian function (Eq. 2), but there is an additional parameter h. It takes account of the fact that along the optical axis, the excitation profile may have a different characteristic length than the transmission profile.

By the term Q(f₁) we have denoted a spatially undefined function that contributes to the first moment of the profile (as a fraction f₁) but not to the higher moments. It takes account of fluorescence from peripheral parts of the sample contributing to fluorescence above the ideal level defined by the first term of Eq. 14. A similar approach has been used earlier (20), but in respect to three-dimensional Gaussian profile. We have experienced that the first term alone would not yield a good fit quality. This is the reason why we have introduced the second term. We consider it as a kind of perturbation; but details of the shape of the perturbation are unknown. We do not know how adequately the second term describes deviations of the true profile from that described by the first term, but it improves fit quality indeed. Artifacts due to this term cannot be excluded, but it is not likely that their magnitude is high. In our examples, the fraction f₁ is typically 6%.

Note that we have defined the spatial brightness distribution as a product of the excitation and the transmission profiles. It does not include the saturation of fluorescence intensity due to uneven triplet population; these effects are taken into account by other parameters of the model.

Precalculated data
To recalculate the distribution P(i,i_z,e,n) in a loop of time-steps, we need the following data that characterize a single time-step. First, each grid point has four neighbors where the molecule may jump to; thus, the four jump probabilities have to be precalculated for each grid point. Secondly, for each grid point and each pair of initial and final electronic state of the molecule, the distribution of the photon-count number is precalculated.

The jump probabilities must fulfill the following conditions. First, the mass center must be conserved. The distribution of a molecule that is located in a particular grid point before a jump, must be centered on the original grid point after the jump. (This condition is relaxed only at the borders of the reservoir.) Secondly, the variance of the distribution after the jump is determined by the diffusion coefficient D and is given by 2Dt in each of the three spatial dimensions (i.e., in each Cartesian coordinate). In case of a constant grid density using Cartesian coordinates, all jump probabilities would be equal. However, the problem becomes more complicated in the case of an uneven grid density and radial coordinates. For the axial z-coordinate at the i^th grid point, the two jump probabilities, up and down, must fulfill the following set of equations (they are mathematical expressions of the above-mentioned conditions):

(15)

The solution of this set of equations is

(16)

In the radial -coordinate, a set of equations similar to Eq. 15 can be written down as

(17)

and the solution is

(18)

From jump probabilities, one may calculate relative volumes that are represented by each grid point. The ratio of volumes represented by two neighbor grid points equals to the ratio of jump probabilities between them. Relative volumes are needed when determining the initial probability distribution of the molecule: the probability that a molecule is in a particular grid point is proportional to the volume represented by the particular grid point.

Further, we precalculate photon-count number distributions emitted by the molecule in a given grid point and a given pair of initial and final electronic states. In calculations of count-number distributions during each time-step, time is a continuous variable, in contrast to our model of diffusion. If we could omit singlet-triplet transitions, the classical intensity of emission would be constant in each grid point, equal to qB(_i,z_j) where (_i,z_j) are grid point coordinates. Any count-number distribution in a time-step t would be Poissonian, with the mean value of qtB(_i,z_j). However, we consider also singlet-triplet transitions. We introduce variable u = t_S/t (which is the fraction of time spent in the singlet state) and apply Eqs. 6, 9, and 10 for each initial state:

(19a)

(19b)

(19c)

(19d)

Expressions above are pairwise-normalized, e.g., With each of the expressions in Eq. 19, we apply Mandel's formula using u as the integration variable. For example,

(20)

We calculate the integrals of the type of Eq. 19 numerically. When using 1400 grid points, we precalculate 4 x 1400 count-number distributions before launching the calculation of cumulative count-number distributions. As an illustration, in Fig. 3 four count-number distributions are plotted that have been calculated for a single grid point near the focus under a realistic set of parameter values.

View larger version (15K): [in this window] [in a new window]   FIGURE 3  An exemplary set of photon-count number distributions calculated for a grid point near the focus and time-step of 0.0690 µs. A set of parameter values that are close to those of the fit curve (that will be described later) has been used in this example.

The first and all consecutive time-steps are similar. First, because our model assumes no movement in space during a time-step t, we calculate a new count number and electronic state distribution for each grid point that takes account of additionally detected photons and singlet-triplet transitions during the time-step:

(21)

Here, S and T denote singlet and triplet states, respectively. At the end of each time-step, diffusion jumps occur:

(22)

It makes sense to keep the first term on the right-hand side of Eq. 22 positive; in fact, this is the criterion for selection of the time-step t.

When the number of time-steps that correspond to the counting time have passed, we sum up the count-number distributions from both final electronic states in each grid point. The outcome is the photon-count number distribution of a single molecule.

From single-molecule to overall count-number distribution
It is possible to convert the photon-count number distribution of a single molecule into an overall count-number distribution of many molecules using the approach which uses a convolution of many single molecule distributions (3). We present here another conversion algorithm.

The probability distribution P(n;T,r) of photon-count numbers emitted by a molecule located in point r at the end of the counting time interval T already represents an average over all possible paths of the molecule that end at r. The number of molecules in a volume element dV is Poisson-distributed,

(23)

The number of photons emitted by a random number of molecules from a volume element is

(24)

Here P_(m)(n;T,r) is the distribution of count numbers from m independent molecules that are all in the same volume element by the end of the counting time interval. This may be expressed as a m-fold convolution of P(n;T,r). However, expressions with convolutions are clumsy and calculation of convolutions is time-expensive; therefore, we avoid them. Let us select the representation of generating functions now. Generally, the generating function of the distribution P(n) is defined as

(25)

The generating function of the distribution given by Eq. 24 can be expressed as

(26)

In the first step, we have used the property that the generating function of m-fold convolution of a distribution equals the m^th power of the generating function of the distribution. In the second step, we have used the identity

We have denoted the generating function of P(n;T,r) by G(;T,r). However, it is useful to take into account that is a generating function of a function that is a modification of P(n;T,r):

(27)

Now we may extend Eq. 19 to the case when all volume is considered. Since all volume elements are independent sources of emission, the contributing generating functions are multiplied. We get

(28)

Note that is the generating function of They do not depend on the selection of the integration volume, provided the volume is large enough to include all parts of space contributing to emission.

Fourier-transform is a special selection of the generating function. To calculate the overall count-number distribution, we start with the calculation of a function describing a single molecule, Afterwards, we apply the Fourier-transform, apply Eq. 28, and finish with the reversed Fourier transformation.

We have assumed a single fluorescent species above. If more than a single species is involved, then their contributions are simply summed up in the exponent on the right-hand side of Eq. 28.

Representation of generating functions has been a cornerstone when deriving theoretical expressions and developing calculation algorithms for FIDA (3), two-dimensional FIDA (11), FIMDA (15), and FILDA (16).

Corrections for afterpulsing and dead-time of the detector
Distortions in photon-count number distributions due to nonideal detectors had been studied long before the first fluorescence fluctuation experiments (24,25). Such distortions have an impact in FFS, too; in particular, the dead-time distortions are stronger at higher count rates (26). Afterpulsing and dead-time are expected to have a strong impact on the shape of the count-number distribution of extremely short photon-counting time intervals, because the characteristic time of these processes is in the range of tens of nanoseconds.

Afterpulsing denotes a nonideal behavior of the detector where an artificial count is generated soon after a real photon pulse. This occurs with a probability that is almost independent of variations in the count rate but characteristic for a given detector. According to manufacturer's specification, the afterpulsing probability of the detector SPCM-AQR (EG&G, Vaudreuil, Quebec, Canada) that we use is 0.3%. If we ignore the delay between the original and the artificial pulse then it makes no difference whether to apply afterpulsing corrections to the single-molecule or multi-molecule count-number distribution. Otherwise, it is preferable to apply it when information about the intensity fluctuations within a counting time interval is available, i.e., when we calculate a single molecule count-number distribution. In this study, we have ignored the delay between the original and artificial pulse and applied afterpulsing correction to the overall theoretical distribution. Provided the number of photon-counts is n and the probability that a photon-count is followed by an artificial count is , the number of afterpulses j follows a binomial distribution,

(29)

As a correction formula, the distribution of true photon-counts P₀(n) and the distribution of counts including afterpulses P_c(n) are related by a convolution,

(30)

Dead-time of our photon detector is 50 ns, according to producer's specification. Compared to afterpulsing, dead-time correction is more difficult in theory. In this study, we have applied an approximated formula,

(31)

Our approximation means that each count number is redistributed by a binomial distribution, Eq. 29, with the probability of missing a count depending on the original count number.

Dead-time causes a drop in the count rate. If we applied Eq. 31 without a modification to a series of count-number distributions corresponding to different T, then the mean count rate after the correction would not be constant but a function of sampling time. To keep all count-rates after the corrections equal, we have replaced _dead in Eq. 31 by effective dead-time as a function of sampling time. The dead-time of the shortest counting time interval T₀ is considered as a primary model parameter, whereas effective dead-time at any other sampling time T is adjusted to yield the same count rate.

If we graphed the calculated photon-count number distributions, then one could hardly notice differences between the results when using different models. It is more informative to plot a characteristic of the shape of the count-number distribution, such as as a function of sampling time. (K_i denotes the i^th factorial cumulant of the count-number distribution.) A better known characteristic of count-number distributions is the apparent count-rate per molecule, Fig. 4 illustrates how the graphs of q_a and ₃ versus the width of the counting-time interval T are modified when singlet-triplet transitions, dead-time, and afterpulsing are taken into account. Diffusion time is kept at 30 µs in all cases.

View larger version (15K): [in this window] [in a new window]   FIGURE 4  Apparent count-rate per molecule (top panel) and a shape characteristic of count-number distribution 3 (bottom panel) versus sampling time calculated for four different sets of parameters. (Dotted line) A model taking account of diffusion but no singlet-triplet transitions, no afterpulsing, and no dead-time. (Solid line) The model taking account of diffusion and singlet-triplet transitions. (Dashed line) The model taking account of diffusion, singlet-triplet transitions, and detector afterpulsing. (Dash-and-dot line) The model taking account of diffusion, singlet-triplet transitions, and dead-time. Values of parameters are selected close to the corresponding values of the final fit curve to experimental data, except the background count-rate is set to zero.

	TEST EXPERIMENTS

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

In experiments, we have saved the photon-count numbers in 2-µs counting intervals, using an ISA-bus counter card built at Evotec. Raw data could be later converted into a set of histograms, each corresponding to a specific counting time interval of 2, 4, 6, 8, 10, 12, 16, 20, 30, 40, 50, and 60 µs. In parallel, FCS data were acquired.

FCS data were fitted using a three-dimensional Gaussian model for the intensity profile. We do not think it is the appropriate model, but this is, by far, the most widely used model in FCS. Thus, we have applied the following formula to fit the normalized autocorrelation function of photon detection:

(32)

Here, m is the mean number of molecules per confocal volume, is the ratio of constant background to fluorescence intensity, is the apparent singlet-to-triplet rate constant, _D is the diffusion time, and s is longitudinal/transverse axial ratio of the brightness profile. Furthermore, from the overall fluorescence count rate and the estimate of m, one may calculate the apparent count rate per molecule,

(33)

However, the apparent count-rate per molecule is numerically different from the count rate of a molecule exactly in the focus (r = 0). Indeed, the overall fluorescence count rate can be expressed as the product of concentration, count rate in the focus, and the first moment of the profile, cq(0)µ₁, whereas the number of molecules in the confocal volume is Brightness in the focus is greater than q_a by a factor of µ₁/µ₂, which is for the three-dimensional Gaussian model.

Furthermore, the amplitude of the triplet term of the correlation function ( in Eq. 32) is not a straightforward triplet/singlet population ratio. Rather, the population ratio is a function of spatial coordinates, and amplitude of the triplet term is a certain average over spatial coordinates. Contribution to correlation function from a volume element is proportional to the square of brightness, therefore

Assuming a three-dimensional Gaussian profile, the conversion formula is

(34)

Now it remains to express the brightness characterizing a molecule in the focus and in the singlet state,

(35)

Note that other profile formulas would yield other values of the numerical coefficients.

We have not corrected FCS data against afterpulsing and dead-time, except that in fitting we have ignored the correlation function at delay times below 0.4 µs. (Distortions are tremendous between time zero and 0.1 µs.) In fact, count rates are affected rather modestly. Under conditions of experiments described below, afterpulsing is expected to increase the count rate by 0.33%, whereas dead-time is expected to reduce count rate by 0.5% at the lowest excitation power up to 1.1% at the highest one.

The results of FCS fitting serve as reference data to evaluate the results of fitting the series of count-number histograms.

When fitting count-number histograms, we have applied a number of different models. A model that has been used as a reference is the original FIMDA model as described in the year 2000 (15). We denote it FIMDA-2000. The ME-theory described above has been used in its full length (ME-FIMDA), but we have applied it also without afterpulsing and dead-time corrections, as well as on the basis of three-dimensional Gaussian brightness profile modified by f₁-parameter.

Whatever model was used for fitting the set of histograms, fit quality has been evaluated using the estimator

(36)

where i is index of the data point (each corresponding to a different histogram and/or different count number), W_i = is weight of a particular data point, is a data point on a histogram (an estimate of probability), is the corresponding data point of the fit function, N^(E) is the sum of the length of histograms, and M is the number of counting time intervals per experiment (which is different for different histograms). The value of ² as defined above is slightly below unity for a good fit quality (the value of 0.75 has been obtained when fitting simulated data), and much above unity for a poor fit quality.

	RESULTS

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
The series of measurements at different concentrations but a single excitation power was undertaken as a test of reasonability of the applied models and theories: ideally, fitting should yield estimates of brightness and diffusion time that are independent of the sample concentration. All three methods that were applied (FCS, FIMDA-2000, and ME-FIMDA) passed this test well. The most remarkable difference is the fit quality between FIMDA-2000 and ME-FIMDA, but this difference is also represented in the measurement series on a single sample but varying the excitation power. This second series has a number of other remarkable results too; therefore, we shall describe its results in detail.

Fit quality when fitting a series of count-number histograms
The value of ² as a quantifier of the fit quality significantly depends on the excitation power, but even more significantly on the model that is used in fitting. On Fig. 5, the mean value of ² is plotted against the excitation power for FIMDA-2000 (as described in (15)) and ME-FIMDA with detector afterpulsing and dead-time corrections (as described in this article). In both cases, some parameters of the model of identical or a similar meaning were fixed in fitting. Background count rate was fixed to 0.333, 0.395, or 0.470 kHz, depending on the excitation power; these values were determined from measurements on pure water. In FIMDA-2000, the axial ratio was fixed to the mean value from fitting FCS data; this is 5.6, 5.8, and 6.3 for the lowest, medium, and highest excitation power, correspondingly. Correspondingly, in ME-FIMDA, the ratio z₀/w₀ was fixed to a rather arbitrarily selected value 2.5. (Fit quality only very weakly depends on this parameter.) Lifetime of the triplet state was fixed to 2.2 µs, which is a mean value obtained by FCS fitting. Afterpulsing probability has been fixed to 0.0033 and dead-time to 0.039 µs. Both these values were estimated from fitting count-number histograms measured under constant illumination of the detector, using the same formulas with the same approximations as when correcting theoretical count-number distributions for fluorescence. At the end, FIMDA-2000 and ME-FIMDA had an equal number of free parameters in fitting: number of molecules, specific brightness, diffusion time, population of the triplet state, and two adjustment parameters of the brightness profile.

View larger version (11K): [in this window] [in a new window]   FIGURE 5  Fit quality parameter 2 for measurements at different excitation power (X axis) and two different models. The lowest power of excitation has been selected as a unity; it corresponds to 60 µW power of the laser beam focused on the sample. To our estimation, this yields a 525-kHz count-rate from a tetramethylrhodamine molecule if in the singlet state and exactly in the focus.

View larger version (31K): [in this window] [in a new window]   FIGURE 6  Weighted residuals of fitting a set of count-number histograms using FIMDA-2000 theory (15). The interval between horizontal lines is worth six units of standard deviation. Residuals of the six histograms that are graphed with solid circles are characterized by 2 above 2.

View larger version (34K): [in this window] [in a new window]   FIGURE 7  Weighted residuals of fitting of the set of count-number histograms using the theory outlined in this article and the model with diffusion, singlet-triplet transitions, afterpulsing, and dead-time.

View larger version (21K): [in this window] [in a new window]   FIGURE 8  A characteristic of the shape of photon-count number distributions, 3, as a function of sampling time. (Circles) Experimental data. (Dotted line) The best fit distributions when using FIMDA theory. (Dashed line) The best fit distribution when using the theory outlined in this article and the model with diffusion and singlet-triplet transitions, but assuming an ideal detector. (Solid line) The best fit distribution by the theory of this article and the model with diffusion, singlet-triplet transitions, afterpulsing, and dead-time.

In summary, all three aspects—sufficiently adequate brightness profile, account of molecular dynamics during the counting time interval, and corrections for detector nonideality—are responsible for a good fit quality of our ME-FIMDA model.

Estimated mean number of molecules
On Fig. 9 the estimated mean number of molecules per confocal volume is plotted as a function of excitation power for FCS as well as for ME-FIMDA. In both cases, confocal volume is defined as

but in the context of the particular FCS model, B(r) is the profile distorted by saturation due to singlet-triplet transitions, while in the context of ME-FIMDA, B(r) is the unsaturated profile. To our understanding, this is the reason why the numerical value of the mean particle number estimated by FCS significantly depends on excitation power, while the dependence is at least three times weaker for ME-FIMDA. It should be pointed out here that this difference is not a sole property of the ME-approach described in this article: a more adequate FCS model taking account of profile distortions (27,28) might yield similar results as well.

View larger version (12K): [in this window] [in a new window]   FIGURE 9  Number of molecules per confocal volume estimated from measurements at different excitation power. FCS and ME-FIMDA data have been plotted. Dashed lines to guide the eye.

View larger version (11K): [in this window] [in a new window]   FIGURE 10  Diffusion time estimated from measurements at different excitation power. FCS and ME-FIMDA data have been plotted. Lines are linear regression curves.

View larger version (12K): [in this window] [in a new window]   FIGURE 11  Triplet/singlet population ratio at r = 0 estimated from measurements at different excitation powers. FCS and ME-FIMDA data have been plotted. Lines are linear regression curves.

View larger version (15K): [in this window] [in a new window]   FIGURE 12  Count-rate per molecule in the singlet state at r = 0 estimated from measurements at different excitation powers. FCS and ME-FIMDA data have been plotted. Lines represent linear functions through data points of the lowest excitation power.

	DISCUSSION

TOP ABSTRACT INTRODUCTION THEORY TEST EXPERIMENTS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
In earlier studies, the theoretical count-number distributions have been expressed as functions of mean particle numbers, mean count-rates per molecule, and shape parameter(s) of the spatial brightness function. Here, a theory has been introduced taking account also of diffusion time, singlet-triplet transition parameters, afterpulsing probability, and dead-time of the detector. We have demonstrated that under experimental conditions typical to FCS studies, all four additional physical mechanisms have a significant influence on count-number distributions. At this point, it is worth discussing how and why in the earlier studies where the relevant processes have been ignored, a rather good fit between a count-number histogram and a theoretical curve has been achieved.

We think that the basic reason for a good fit quality has been that the model describing the spatial brightness profile has involved adjustment parameter(s). In fact, the adjustment parameters of the brightness profile serve as a means to compensate for the absence of some other physically relevant processes in the model. The fact that the profile parameters describe more than just the profile became very obvious when fitting a series of histograms corresponding to different sampling times (15).

An important outcome of this study is the realization that under typical conditions of fluorescence fluctuation experiments, diffusion, singlet-triplet transitions, afterpulsing, and dead-time are all significant determinants of the count-number distribution. It seems that attempts to derive analytical expressions for the count-number distributions are fruitless, but numerical methods are of a practical value.

A drawback of this study is a low speed of the introduced calculation algorithm. Fitting a series of count-number distributions lasts tens of minutes. We have presented a theory that has been successfully tested in first experiments, but it is not yet a calculation algorithm suited for high throughput applications. Slow calculation is reported to be an issue in Monte Carlo calculations of PAID, too (19). Presently we believe that further modification of calculation algorithms will much more efficiently reduce calculation times than the steady increase in the speed of computers.

A message of this study is that simple theories, which ignore the intensity fluctuations during counting time interval, can be successfully applied, but attempts to use the true spatial brightness function may turn out even less fruitful than representing the function by an empirical formula. Parameters of the apparent brightness profile are, in fact, functions of a number of different physical properties of both the equipment and the sample. In a number of cases, a simple theory is good enough for data analysis, but not always. Simultaneous fitting of a series of count-number histograms is an example when profile-driven theories are unsatisfactory. In these cases, the theory introduced here is particularly useful.

In the approach that we have introduced, the theoretical curve is calculated numerically. Generally, numeric calculation of the fit curve is slower and more sophisticated compared to fitting data to analytical formulas, but it has its favorable properties, too. For