1. Introduction

Fluorescence lifetime imaging microscopy (FLIM) can be used to map the nanosecond fluorescence lifetimes of fluorophores in single cells [1]. This provides information about the state of the fluorophores and their molecular environment. The fluorescence lifetime is sensitive to environmental conditions such as pH [2], and excited state reactions such as fluorescence resonance energy transfer (FRET) [3], properties that can be exploited to resolve physiological parameters in the cell [4, 5]. With the advent of fluorescence proteins, FRET has become the tool of choice for detecting protein-protein interactions in living cells. FLIM has proven to be one of the most reliable approaches to measure FRET, and is the only method that allows quantification of the relative concentrations of interacting species [5].

Time correlated single photon counting (TCSPC) has been employed with confocal and two-photon microscopy to measure the arrival times of the emitted photons in response to pulsed excitation [6, 7, 8]. The fluorescence lifetimes and relative concentrations of the fluorescent species can be derived by fitting a sum-of-exponential model to the acquired decay curve. However, at least 1000 photons per FLIM curve are needed to distinguish two lifetimes [9]. In biological applications typically an order of a magnitude less photons per pixel are collected, severely limiting the quantitative analysis of the data. An additional complication is that for a quantitative fit, the shape of the impulse response function (IRF) of the instrument should be know, which requires separate calibration measurements and complicates data fitting.

Global analysis methods can increase the accuracy and precision of the analysis of TCSPC data [10, 11, 12]. Such methods were first applied in microscopy to FLIM data acquired by frequency domain systems [13, 14, 15, 16]. In this approach, it is assumed that the lifetimes of all fluorescence species in the sample are invariant and that only their relative concentrations are different in each pixel of the FLIM image. It was shown that this provides sufficient information to resolve a bi-exponential decay model, and that accuracy and precision of the estimated parameters could be increased significantly. This approach has been used successfully in biological applications to quantitatively determine the fraction of interacting molecules with FRET assays in fixed and living cells [14, 17, 18, 19, 20].

Previous work has shown that Fourier domain methods originally developed for frequency domain FLIM data can be used to analyze time-domain FLIM data [21, 22, 23]. Here we show that single frequency domain global analysis methods can be applied to time-domain data to accurately fit bi-exponential models. In addition, a new filtering method optimized for FRET-FLIM data is proposed which greatly enhances the robustness of the global analysis method. We show that this yields a robust method to recover the relative concentrations of interacting species from TCSPC data in a FRET assay. In addition to these global analysis methods, we introduce a simple method to estimate the required parameters of the IRF from the data thereby avoiding the use of separate reference measurements. We validate our methods using simulated data, and demonstrate their use with experimental data. We show that the proposed method is capable of fitting bi-exponential curves to FRET- FLIM data recorded with photon counts that are one order of a magnitude lower than previously recommended for fitting TCSPC data.

2. Theory

In this section we describe the theory of analyzing time-domain data with Fourier analysis methods. The reader is referred to the excellent textbook by Lakowicz for an in-depth theory of time-domain and frequency domain FLIM [1], and we limit ourselves to the equations that are essential for the development of a Fourier-based global analysis method of time-domain data.

2.1. Fourier analysis of time-domain FLIM data

The decay kinetics given by a sum of Q exponentials with fluorescence lifetime τ_q:

where the sum of the fractional fluorescence α_q is normalized to one: ∑_q=1 ^Q α_q = 1. We assume repetitive excitation with a period T and write the excitation E(t) as a Fourier series:

where ω= 2π/T. E ₀ is real and E_n = E_n ^*, since E(t) is a real-valued signal. The time-resolved response to the excitation defined by Eq. (2) is given by

where F_T is the total fluorescence intensity. The normalized result can be written as a Fourier series:

2.2. Global analysis of bi-exponential time-domain FLIM data

Global analysis of FLIM data is achieved by fitting a set of decay curves to Eq. (4), assuming that the lifetimes τ_q that are invariant. At this point it is useful to define the quantity R_n,q:

which is the Fourier coefficient at harmonic number n, for the single-exponential component q. We assume a set of measurements, indexed with the superscript i, that all have the same fluorescence lifetimes τ_q, and rewrite Eq. (4) as

where

are the normalized Fourier coefficients of the data, corrected for the instrument response E_n. For the special case of a bi-exponential decay we can write explicitly:

where we substituted α ₂ ⁱ = αⁱ and α ₁ ⁱ = 1-αⁱ, since α ₁ ⁱ + α ₂ ⁱ = 1. We selected αⁱ to correspond to the fractional fluorescence of the second species for consistency with the analysis of FRET-FLIM data described below. The real and imaginary components of Eq. (8) form two equations from which αⁱ can be eliminated yielding a linear relation between the components of Rⁱ_n:

where

The global analysis of the data proceeds by fitting the data Rⁱ_n to the linear equation (9) and solving Eqs. (10) for R _n,1 and R _n,2. Figure 1 shows graphically the principle of this approach using a phasor plot of ImRⁱ_n against ReRⁱ_n. All points Rⁱ_n derived from mono-exponential curves will fall on a half-circle [16, 24]. All data points of a two-component mixture will, in the absence of noise, be found on a straight line. The intersections R _n,1 and R _n,2 of this line with the half-circle define the two mono-exponential lifetimes of the mixtures. The lifetimes then follow from Eq. (5):

Equation (8) suggests that we can trivially solve for αⁱ if R _n,1 and R _n,2 are known, but this solution is not optimal for use with experimental data. As shown in Fig. 1, for arbitrary noisy points a better estimation of αⁱ is found by the normalized scalar projection of Rⁱ_n on the line connecting R _n,1 and R _n,2. This turns out to be equivalent to a least squares estimation of αⁱ, as was derived before by Verveer et al. [15] and yields:

 

Fig. 1. Graphical representation of global analysis of FLIM data using a phasor plot [24]. The half-circle represents all possible values of mono- exponential lifetimes, whereas the straight-line segment represents all possible mixtures of two mono-exponential species. Global analysis proceeds by estimating the straight line from data and the intersections with the half-circle (R _n,1, R _n,2, open square and circle). For a noisy data point, (Rⁱ_n, filled circle), the fluorescence fraction is given by the projected distance to R _n,1.

Download Full Size | PDF

Both Eqs. (11) and (12) were derived for frequency domain FLIM data before by us [15, 16], and later confirmed by others [25, 26]. We can apply this approach directly to the normalized Fourier coefficients of time domain FLIM data Rⁱ_n that can be estimated by Fourier methods if E_n is known, for instance from a reference sample. Such an analysis is based on a single harmonic, even if multiple harmonics are generally available for TCSPC data. This implies that information in the other harmonics is not used, but nevertheless this is sufficient for the important case of fitting bi- exponential decay curves. In principle, any of the harmonics could be used to calculate the two lifetimes τ ₁ and τ ₂, and the relative fractions αⁱ. In this work we have used the first harmonic, since for TCSPC data with pulse-like excitation it is the strongest harmonic, and therefore is expected to have the best SNR.

We generally collect data at the highest temporal resolution available to us (8 ps, see also materials and methods, section 3.5). The temporal resolution can be lower, as long as the decay curve is sufficiently sampled to avoid aliasing artifacts. In cases where the resolution is inherently low, which is true for some systems that employ gated detection, the harmonic content of the signal should be evaluated before applying a Fourier domain analysis.

2.3. The impulse response function

To analyze time-domain data, the impulse IRF of the instrument should be taken into account, implying that the Fourier coefficients E_n of the excitation E(t) must be estimated. This can be done using a sample with a known fluorescence lifetime τ(or a scattering sample, where τ = 0). In this case, Eq. (6) is used to estimate E_nR_n by Fourier analysis, and since R_n can be calculated from Eq. (5) using the known fluorescence lifetime, the values of E_n follow. Here we propose an alternative approach that does not require a reference sample. We assume that the IRF can be approximated by a Gaussian function with standard deviation σ and a delay t ₀:

 

Fig. 2. The argument and the absolute value of the Fourier coefficients of TCSPC FLIM data. The Fourier coefficients for a simulated FLIM decay curve (τ ₁ = 3 ns, τ ₂ = 1.5 ns, and α = 0.5), were calculated. The parameters of the Gaussian IRF were t ₀ = 2.0 ns and σ= 0.1 ns. (a) Argument of the Fourier coefficient as a function of the harmonic number n. (b) Logarithm of the absolute value of the Fourier coefficient as a function of n ². The straight lines are linear fits to the higher harmonics (n > 20, dotted lines).

Download Full Size | PDF

Convolution with a train of delta pulses with a fundamental frequency ω and Fourier transformation gives the Fourier coefficients E_n of a delayed train of Gaussian pulses:

To estimate E_n from the product E_nR_n, we assume mono-exponential decay kinetics, and use Eqs. (5) and (14) to write the argument and the absolute value of E_nR_n as:

The argument depends only on the delay t ₀, and the absolute value depends only on the width of the pulse σ. In Fig. 2(a) it can be seen that for higher values of n, the curve for arg(E_nR_n) is nearly linear. Likewise, Fig. 2(b) shows that the logarithm of ∣E_nR_n∣ is a near linear function of n ², for high values of n. We can obtain approximate linear expressions of arg(E_nR_n) and ln∣E_nR_n∣, as a function of n and n ², respectively, by Taylor expansion around an appropriately large harmonic number n ₀:

If n ₀ is sufficiently large, the contribution of the terms that contain τ can be ignored compared to the terms that contain only t ₀ or σ, and we find

where a and b are constant values. Thus t ₀ and σ can be found by linear fits to the higher harmonic Fourier components. In the case of multi-exponential decays this approach is still valid, since Eqs. (16) are extended with additional terms depending on the fluorescent lifetimes, and will still tend to be linear for high values of n.

Figure 2 gives an example for simulated bi-exponential FLIM data, with the IRF parameters given by t ₀ = 2 ns and σ = 0.1 ns. Linear fitting of the higher Fourier harmonics yielded estimations of t ₀ = 2.01 ns and σ = 0.188 ns. Thus, t ₀ was reliably estimated using the linear approximation, but the estimation of σ was inaccurate. This is because the neglected contributions containing the fluorescence lifetimes are larger in the approximation of ln∣E_nR_n∣ compared to the approximation of arg(E_nR_n). Rather than extending the range of n, which might not be feasible with experimental data, we used a non-linear fit of Eq. (15b) to estimate σ more reliably. In this fit, the value of t ₀ was fixed to the value obtained from the linear fit, the initial value of σ was also taken from the linear fit, and an initial value of τ was derived using Eq. (15a). Although the actual data was bi-exponential, a value of σ = 0.099 ns was fitted. Thus, reliable calibration values of t ₀ and σ can be obtained directly from the FLIM data without additional calibration samples.

2.4. Global analysis of FRET-FLIM data

In the case of typical FRET-FLIM data of two interacting species of molecules it is possible to make use of additional information about the sample and apply global analysis methods. In a FRET-FLIM system the lifetime of one fluorophore, the donor, is measured. If the donor-tagged molecule interacts with another fluorophore-tagged molecule, the acceptor, the excitation energy of the donor can be transferred to the acceptor fluorophore, provided the conditions for FRET are met [3]. As a result of FRET the fluorescence lifetime of the donor is decreased to a lower value. If the binding geometry of the complex can be considered to be relatively constant, the fluorescence lifetime of the complex will also be constant in the sample. Thus, in this case we can describe the system as a mixture of two molecular species, free donor and donor/acceptor complexes, each with unique fluorescence decay characteristics. If the individual decays are mono- exponential, then the fluorescence kinetics are described by a bi-exponential model, a reasonable assumption for selected donor fluorophores. The fluorescence lifetimes are the same in each pixel, and only the relative fraction of each species varies, a system that is amendable to global analysis. Provided that the FLIM system is properly calibrated for the instrument response, no additional measurements with references samples are needed in this case. This is an advantage over intensity based FRET methods that require separate calibration samples [27]. Moreover, unlike other FRET methods, the recovered lifetimes yield an estimate of the quantum yield of both free and complexed species, making it possible to renormalize the estimated fractional fluorescence to true relative concentrations [13].

 

Fig. 3. Geometrical representation of the filtering of FRET-FLIM data. The imaginary parts of the first harmonic of the FLIM data are plotted against the corresponding real parts. From donor-only data, a Gaussian distribution of the donor Fourier components is estimated, plotted as a blue contour at a preset distance from its mean. Points outside the two red dotted straight lines are rejected. Points within the donor distribution contour are assumed to be donor-only points and are also not used. Black circles: rejected points, blue circles: donor points, red circles: accepted points.

Download Full Size | PDF

Global analysis has been applied to FRET-FLIM data collected with wide-field frequency domain systems successfully in biological applications [14, 17, 18, 19, 20]. One reason for this success is that the signal-to-noise ratio in such systems is sufficiently high. As discussed above, these global analysis methods are also applicable to TCSCP data, but these are much more likely to suffer from a low SNR, since at most a few hundred photons can be collected per pixel in a reasonable time-frame, due to the requirement to record data at rates well below the saturation value of the detector [1]. Global analysis alleviates this problem since all photon counts of one or more images are considered simultaneously. Nevertheless, compared to wide-field data, the low SNR of TCSPC data poses a problem.

2.5. Global analysis of TCSPC data with low SNR

For low SNR FRET data it is advantageous to measure samples that only contain donor. Donoronly samples are not strictly needed for the global analysis, but they can provide information on the lifetime of the free donor, and on the SNR of the data. Donor-only samples are often easy to obtain by omitting the acceptor during sample preparation. In some cases, this is not straightforward, for instance in the case of sensors that consist of linked donor and acceptor fluorophores. But even then it is often possible to construct a control sensor without acceptor, or to remove the acceptor before the measurement by photobleaching [28]. Thus, from a number of donor-only images, the distribution of the normalized Fourier coefficients of the donor, R _n,1, can be estimated, and plotted in a phasor plot as a two-dimensional distribution. The normalized Fourier coefficients Rⁱ_n from a sample with donor and acceptor must be found on the straight line that intersects the half-circle of single-exponential values, as illustrated in Fig. 1.

To increase the reliability of the global fit, noisy data points can be filtered by applying physical constraints. A value of Rⁱ_n is physically acceptable if a straight line through Rⁱ_n can be found that passes through the donor distribution and that intersects the half-circle at any point to the right of the donor distribution. Strictly, we should use the probability distribution of the mono-exponential donor coeffiecients on the half-circle to calculate a confidence value. In practice, this is difficult with the available distribution of R _n,1, which includes multi-exponential points due to noise variation. Instead we adopt a simple geometrical scheme that is illustrated in Fig. 3. We plot a contour line from the Gaussian distribution of the donor coefficients at a given confidence level (for instance at one standard deviation) and find the intersections of this contour with the half-circle of mono-exponential values. Any point Rⁱ_n below the line connecting the left-most intersection with the zero-lifetime point on the half-circle, R_n = (0,1), cannot belong to a linear mixture unless its donor lifetime is unlikely to occur, or unless its second lifetime is negative. Likewise any point Rⁱ_n above the line through the left-most intersection, tangent to the half-circle, cannot be part of linear mixture unless its donor lifetime is unlikely or its second lifetime is not intersecting the half-circle (i.e. is not mono-exponential). Thus, we filter the donor/acceptor data by rejecting points that not fulfill these criteria. Furthermore, we exclude points that belong to the donor distribution (i.e. within the donor probability contour) from the set of donor-acceptor points. The result of the filtering is illustrated in Fig. 3, where black circles indicate rejected points. For simplicity, we select the donor probability contour as a circle at a given probability level of a symmetrical 2D Gaussian distribution, which is estimated from the data. The procedure is easily adapted for arbitrary donor distributions, but we found this to be unnecessary. By adjusting the probability level, the filtering can be made less or more stringent. The fluorescence lifetimes can be found from the remaining donor/acceptor points as described before. However, the additional knowledge of the donor distribution can be used to further improve the fit by constraining the mixture line to pass through the estimated mean of the donor coefficients. For an appropriately chosen mono-exponential fluorophore, the distance of the donor mean to the mono-exponential half-circle will be well within the error defined by its distribution, thereby validating its choice as a constraint for the linear fit.

3. Materials and methods

3.1. Simulations

Simulated data were generated using the Matlab data processing package (The MathWorks, Natick, MA, USA). To generate a set of simulated curves, a Gaussian pulse was generated and convolved with a bi-exponential decay curve with two fixed lifetimes and a relative fraction that was drawn from a random distribution with a mean of 0.4 and a standard deviation of 0.3. To avoid border and sampling artifacts, the curve was generated with a length of 100 ns and a resolution of 0.8 ps, and then cropped and rebinned to a length and resolution equal to our experimental data (25 ns and 8 ps, respectively). After normalizing each curve to a preset total of counts, noise was added by replacing the counts in each bin by a value drawn from a Poisson distribution. Donor-only data were simulated in a similar fashion, assuming single-exponential decays. The lifetimes and the relative fractions were then recovered using the data analysis procedure described below. Where appropriate, simulations were repeated 500 times with different noise realizations, and the means and standard deviations of the recovered parameters were calculated and reported.

3.2. Data analysis

The proposed method was implemented using the Matlab data processing package (The Math-Works, Natick, MA, USA). The process consists of the following steps:

In this work we selected n = 1 in steps 3–7.

3.3. Sample preparation

MCF7 breast tumor cells (American Type Culture Collection, Manassas, VA, USA) were grown in DMEM (Sigma-Aldrich Biochemie GmbH, Hamburg, Germany) supplemented with 10% FCS and 2 mM glutamine (Sigma-Aldrich Biochemie GmbH, Hamburg, Germany). Cells were transiently transfected with an EGFR-YFP plasmid using Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA). Cells were allowed to express the protein for 24 hours after which 1 mg/ml of the antibiotic G418 (Serva electrophoresis Gmbh, Heidelberg, Germany) was added in order to select stable transfectants. Before sample preparation, transfected cells were seeded on glass bottom Lab-Tek chambers (Nalge Nunc International, Rochester, NY, USA), grown overnight, and then starved for 5 hours in DMEM supplemented with 0.5% FCS. The cells were stimulated with 100 ng/ml of epidermal growth factor (EGF, Cell Signalling Technology, Inc. Danvers, MA, USA) for 5 minutes, fixed with 4% paraformaldehide in PBS, washed three times with TBS for 5 minutes, permeabilized with 0.1% TX-100 for 5 minutes and washed twice in PBS. The samples were then incubated with 15 μg/ml of PY72, a generic antibody against phosphotyrosine (In vivo Biotech Service, Hennigsdorf, Germany) labeled with Cy3.5. Donor-only samples were prepared by skipping this last step.

3.4. Measurement of the impulse response function

The impulse response function of the PicoQuant LSM Upgrade kit was measured directly by measuring the light reflected from a mirror surface, using the same settings of the system as used in the experiments but replacing the dichroic mirror by a 20/80 beam splitter and the emission filter by a OD 3 neutral density filter.

3.5. TCSPC measurements

Confocal TCSPC images were collected using a LSM Upgrade Kit (PicoQuant, Berlin, Germany) attached to an FV-1000 microscope (Olympus Deutschland GmbH, Hamburg, Germany). A 470 nm pulsed diode laser (LDH 470, PicoQuant, Berlin, Germany) was used as an excitation source. A UMYFPHQ dichroic (Olympus Deutschland GmbH, Hamburg, Germany) and a 525/50 HQ filter were used to detect the emitted photons using a Single Photon Avalanche Photodiode (SPAD) with a temporal resultion of 8 ps. Pixel-by-pixel and tail fits of the data were done with the SymPhoTime software (PicoQuant, Berlin, Germany). The relative concentration of phosphorylated EGFR-YFP in the pixel-by-pixel fit was obtained by dividing the amplitude of the short lifetime component by the sum of the two amplitudes.

 

Fig. 4. Recovery of IRF parameters from simulated TCSPC FLIM data. (a) Mean value and standard deviation of the estimated value of the offset t ₀ as a function of mean photon count per pixel, for different total numbers of pixels. (b) Mean value and standard deviation of the estimated value of the pulse width σ as a function of mean photon count per pixel, for different total numbers of pixels.

Download Full Size | PDF

4. Results

4.1. Simulations

Numerical simulations were performed to validate the proposed method. To study the ability to recover the simulated parameters as a function of the noise, homogeneous intensity images where generated with different mean number of photons per pixels. To investigate the role of the total number of pixels, the number of images used in the analysis was varied. The data sets used in the simulations presented here consisted of 1, 4, or 16 pairs of donor and donor/acceptor images of 32×32 pixels. Figure 4 shows the recovered parameters of the IRF as a function of the mean count per pixel, for three different image sets. Shown are mean value and the standard deviation of t ₀ and σ, calculated from 500 simulations with different noise realizations. The recovery of the IRF parameters shows good agreement with simulated parameters if the mean photon count per pixel is sufficiently high and improves if the number of pixels increases. Thus, these simulations indicate that it is possible to reliably recover the parameters of the IRF from the data, without additional calibration samples.

 

Fig. 5. Global analysis of simulated TCSPC FRET-FLIM data. (a) Mean value and standard deviation of the short lifetime estimated by global analysis, as a function of mean photon count per pixel, for different total numbers of pixels. (b) Fractional fluorescence error as a function of mean photon count per pixel, for different total numbers of pixels. The fractional fluorescence error is defined as the mean of the absolute difference between the recovered and simulated α.

Download Full Size | PDF

Figure 5 shows the result for the recovery of the short lifetime and fractional fluorescence error for the proposed global analysis method compared to a standard global fit. In the first case, the proposed method includes the filtering of data points as described above while using the estimated parameters of the IRF. The second case simulates application of standard frequency domain global analysis, without any filtering, using the known parameters of the IRF. The fractional fluorescence error is defined as the mean of the absolute difference between the recovered and simulated α. As expected, the fit improves as the mean number of counts per pixel increases. At a mean count of a 1000 photons or higher, the filtering procedure is not necessary. With filtering, a reliable fit can be obtained even if the mean photon count per pixel is much lower, removing the systematic error that is observed in the absence of filtering. As expected for a global analysis method, the quality of the result improves as the number of pixels increases. Note also that the errors in the estimated IRF parameters are small and do not prohibit successful global analysis.

 

Fig. 6. Recovery of IRF parameters from experimental TCSPC FLIM data. (a) Argument of the Fourier coefficient of the averaged FLIM data as a function of the harmonic number n. (b) Logarithm of the absolute value of the Fourier coefficients as a function of n ². The red lines are linear fits to the higher harmonics (n > 23, blue dotted lines). (c) Plots of the measured IRF, the Gaussian pulse with parameters derived estimated from the estimated Fourier coefficients, and a Gaussian pulse with parameters derived from the first harmonic of the measured IRF.

Download Full Size | PDF

4.2. Application to FRET in cells

We applied the new global analysis methods to FRET data obtained by a well-established biological FRET-FLIM assay [14, 28]. MCF7 cells expressing the epidermal growth factor receptor tagged with yellow fluorescent protein (EGFR-YPF) were fixed and incubated with Cy3.5-labeled PY72, a generic antibody against phosphotyrosine. In our previous work we used wide-field frequency domain FLIM, combined with global analysis to quantitatively image the relative concentration of phosphorylated EGRF in each pixel [14]. Figures 6 and 7 show the results obtained with a confocal TCSPC system and analyzed with the new methods described in this paper. A total of 8 images from donor-only samples and 8 images from samples with acceptor were acquired. Figure 6(a) shows the argument of the Fourier coefficients of the average curve obtained from the data as a function of the harmonic number n. For sufficiently high value of n these data can be well approximated by a linear model. Likewise Fig. 6(b) shows that the logarithm of the absolute value as a function of the square of the harmonic number n is linear for sufficiently high values of n. Figure 6(c) displays the Gaussian pulse shape using the estimated values for delay t ₀ and the width σ. Also shown in Fig. 6(c) is the experimentally measured IRF of the system. The measured shape of the IRF is distinctly non-Gaussian, and a fit of a bi-exponential model should account for this. The global analysis only uses the first harmonic of the data and therefore we checked if the first harmonic of the measured IRF corresponded to that of the estimated Gaussian pulse. Indeed the value of E ₁ calculated from Eq. (14) corresponded well to the value obtained by Fourier transformation of the measured IRF (0.895 + 0.446j vs. 0.892 + 0.450j). To visualize this we reconstructed a Gaussian pulse from the first harmonic of the measured IRF using the inverse of Eq. (14). The Gaussian pulse estimated from the data corresponds to the reconstructed Gaussian curve with high accuracy (Fig. 6(c)).

 

Fig. 7. Global analysis of experimental TCSPC FRET-FLIM data. (a) Fluorescence intensity of EGFR-YFP. The colorbar indicates the number of counts per pixel. (b) Relative concentration of phosphorylated EGFR-YFP estimated from a pixel-by-pixel bi-exponential fit. (c) Tail fit to the average of all FLIM curves for two different regions (5-20 ns and 6-21 ns). (d) Phasor plot of a random subset of all donor-only data, with the estimated Gaussian distribution indicated by a contour at one standard deviation from the mean. (e) Phasor plot of a random subset of all data from donor/acceptor samples, with points that are accepted by the filtering procedure indicated in red color. The straight line represents a linear fit to the accepted data from which the two lifetimes are calculated. (f) Relative concentration of phosphorylated EGFR-YFP estimated by global analysis.

Download Full Size | PDF

Figure 7(a) shows the donor fluorescence of a single cell from a sample that was fixed 5 minutes after stimulation with soluble epidermal growth factor (EGF) and incubated with Cy3.5-labeled PY72. A pixel-by-pixel fit failed to recover the proper distribution of the relative concentration of phosphorylated EGFR (Fig. 7(b)), due to the low photon counts in the individual pixels (Fig. 7(a)). It is possible to recover the two lifetimes using a tail fit of the averaged curve, but the results depend strongly on the region of the fit that is chosen as is demonstrated in Fig. 7(c).

We therefore applied the new global analysis methods to these data. Figure 7(d) shows the phasor plot of the Fourier coefficients of the donor-only data together with a contour at one standard deviation of the estimated Gaussian distribution. The mean of the donor distribution is located close to the half-circle of mono-exponential points, well within the variation of the data, and yields a lifetime estimation of τ ₁ = 3.04 ns. This is comparable to values obtained by wide-field frequency domain FLIM (2.84 ns and 2.87 ns, from phase and modulation, respectively). A standard frequency domain global analysis [16] of the data, using the estimated parameters of the IRF, yielded lifetime values of 3.28 ns and 1.84 ns. This is not consistent with the estimated donor lifetime and does not correspond well to a global analysis of wide-field frequency domain FLIM (2.92 ns and 1.34 ns). As the photon count per pixel was extremly low, we applied the improved global analysis method described in section 2.5 to the data. Figure 7(e) shows a phasor plot of a subset of the donor/acceptor data, with the fit to the accepted data points, yielding fluorescence lifetimes of 3.04 ns and 1.51 ns. Figure 7(f) shows the corresponding relative fractions of phosphorylated receptor, demonstrating that most receptors are active near the periphery of the cell, as expected under these conditions, and in agreement with previous work [14, 28].

5. Discussion

Most data analysis methods for time-domain FLIM directly fit the data curve. Although that makes it theoretically possible to fit complex multi-exponential decay models, this is in practice difficult with the noisy data encountered in biological applications. Our approach uses the fact that in many biological applications, such as FRET assays, prior knowledge about the system is available. This allows us to limit our model to bi-exponential curves and apply single frequency global analysis to the first harmonic of the data. This implies that the information content in the higher harmonics is not used, but our results show that the use of the first harmonic is sufficient for the purpose of analyzing FRET data.

Even with global analysis methods, the analysis of confocal TCSPC data is problematic as the SNR tends to be low. For this reason we developed a new filtering method for FRET data. This filtering method requires donor-only samples, which can be obtained with little effort in most FRET assays, under the same conditions as the donor/acceptor data. Simulated and experimental data confirmed the effectiveness and robustness of this filter. This approach is equally suitable to data acquired by frequency domain FLIM instrumentation.

We developed a new method to estimate a Gaussian approximation of the IRF directly from the data. The first harmonic of the estimated pulse matched that of the experimental IRF with high accuracy even if the latter was distinctly non-Gaussian. Since the global analysis only uses the first harmonic of the data, the higher harmonics of the IRF do not affect the result and need not to be estimated accurately. This is a crucial advantage over deconvolution in the time-domain, which depends on all harmonics and is a notoriously ill-posed problem. This is especially important when employing pulsed diode lasers that have a distinctly non-Gaussian IRF. Eliminating the need for a separate calibration substantially simplifies data analysis, and avoids problems introduced by changes in experimental or environmental conditions.

Previously published global analysis methods for time-domain data are based on non-linear fitting of the data [10, 11, 12]. These non-linear methods require initial guesses for the solution, a problem that was solved in different ways by these authors. In our approach, which is linear at its core, no initial guesses are necessary, considerably simplifying implementation and use. Non-linear minimization methods are time- and memory consuming and therefore computationally expensive. In our approach, all calculations are simple operations and the data do not need to be kept in memory simultaneously. The authors of the non-linear methods report applications to data with relatively high photon counts. Pelet et al. do not show FRET data, but stress the need for high photon counts in FRET applications [10]. Barber et al. report FRET data with 640-3840 counts per pixel [12]. In contrast, our results show that the number of photons may be an order of a magnitude lower than these reported values (Fig 7, 20-250 counts per pixel).

Our method enables the routine use of confocal TCSPC systems to measure FRET with low acquisition times and minimal photon counts, an essential advantage in live cell assays. The current method covers many applications of biological interest, such as quantitative FRET assays that are used to detect protein interactions in intact cells. One topic of future research is to exploit the information from the higher harmonics to further improve the quality of the global fit. Furthermore, future developments should make it possible to measure multi-exponential decays and explore complex interactions between two or more proteins in the cell.