## Abstract

Fluorescence lifetime imaging microscopy (FLIM) can be used to quantify molecular reactions in cells by detecting fluorescence resonance energy transfer (FRET). Confocal FLIM systems based on time correlated single photon counting (TCSPC) methods provide high spatial resolution and high sensitivity, but suffer from poor signal to noise ratios (SNR) that complicate quantitative analysis. We extend a global analysis method, originally developed for single frequency domain FLIM data, with a new filtering method optimized for FRET-FLIM data and apply it to TCSPC data. With this approach, the fluorescent lifetimes and relative concentrations of free and interacting molecules can be reliably estimated, even if the SNR is low. The required calibration values of the impulse response function are directly estimated from the data, eliminating the need for reference samples. The proposed method is efficient and robust, and can be routinely applied to analyze FRET-FLIM data acquired in intact cells.

© 2009 Optical Society of America

## 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}*α*= 1. We assume repetitive excitation with a period

_{q}*T*and write the excitation

*E*(

*t*) as a Fourier series:

where *ω*= 2*π*/*T*. *E*
_{0} 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 *α*
_{2}
* ^{i}* =

*α*and

^{i}*α*

_{1}

*= 1-*

^{i}*α*, since

^{i}*α*

_{1}

^{i}+

*α*

_{2}

*= 1. We selected*

^{i}*α*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

^{i}*α*can be eliminated yielding a linear relation between the components of

^{i}*R*:

^{i}_{n}where

The global analysis of the data proceeds by fitting the data *R ^{i}_{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 Im

*R*against Re

^{i}_{n}*R*. All points

^{i}_{n}*R*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

^{i}_{n}*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 *α ^{i}* 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

^{i}*R*on the line connecting

^{i}_{n}*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

^{i}*et al*. [15] and yields:

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 ^{i}_{n}* that can be estimated by Fourier methods if

*E*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

_{n}*τ*

_{1}and

*τ*

_{2}, 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.

^{i}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*by Fourier analysis, and since

_{n}R_{n}*R*can be calculated from Eq. (5) using the known fluorescence lifetime, the values of

_{n}*E*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

_{n}*σ*and a delay

*t*

_{0}:

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*, we assume mono-exponential decay kinetics, and use Eqs. (5) and (14) to write the argument and the absolute value of E

_{n}R_{n}_{n}R

_{n}as:

The argument depends only on the delay *t*
_{0}, 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 _{n}R_{n}*) is nearly linear. Likewise, Fig. 2(b) shows that the logarithm of ∣

*E*∣ is a near linear function of

_{n}R_{n}*n*

^{2}, for high values of

*n*. We can obtain approximate linear expressions of arg(

*E*) and ln∣

_{n}R_{n}*E*∣, as a function of

_{n}R_{n}*n*and

*n*

^{2}, respectively, by Taylor expansion around an appropriately large harmonic number

*n*

_{0}:

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

where *a* and *b* are constant values. Thus *t*
_{0} 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*
_{0} = 2 ns and *σ* = 0.1 ns. Linear fitting of the higher Fourier harmonics yielded estimations of *t*
_{0} = 2.01 ns and *σ* = 0.188 ns. Thus, *t*
_{0} 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 _{n}R_{n}*∣ compared to the approximation of arg(

*E*). Rather than extending the range of

_{n}R_{n}*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*

_{0}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*

_{0}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].

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 ^{i}_{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 ^{i}_{n}* is physically acceptable if a straight line through

*R*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

^{i}_{n}*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*below the line connecting the left-most intersection with the zero-lifetime point on the half-circle,

^{i}_{n}*R*= (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

_{n}*R*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.

^{i}_{n}## 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:

- TCSPC curves are extracted from several data sets, from donor-only samples and from samples with acceptor. Pixels with a total number of photons less than a preset threshold of 20 counts are excluded from the analysis.
- The parameters of the IRF,
*σ*and*t*_{0}, are calculated from the Fourier transform of the average of the selected curves, as described in section 2.3. - The Fourier coefficient of the
*n*th harmonic is calculated from each TCSPC curve:where

*b*s the number of counts in the_{k}^{i}*k*th bin of the histogram of photon counts acquired at pixel*i*, and*B*is the total number of used bins. Using Eq. (14) and*σ*and*t*_{0}these Fourier coefficients are corrected for*E*to obtain_{n}*R*.^{i}_{n} - The distribution of
*R*_{n,1}of the donor-only data is estimated by a weighted mean and standard deviation of the real and imaginary parts, to obtain*R*̄_{n,1},*σ*Re*R*_{n,1}and*σ*Im*R*_{n,1}. The weight in each pixel is chosen equal to the square root of the number of photons, which is an estimator for the SNR of a Poisson process (SNR =*N*/√*N*, with N the mean and variance of the Poisson process). - All points from images with donor and acceptor are filtered according to the procedure described in section 2.5, assuming a symmetrical distribution of the donor with mean
*R*̄_{n,1}and standard deviation equal to the average of the estimated standard deviations:*σ**R*_{n,1}= (*σ*Re*R*_{n,1}+*σ*Im*R*_{n,1})/2. - A straight line Im
*R*=^{i}_{n}*u*+_{n}*v*Re_{n}*R*is fit through all remaining donor/acceptor points^{i}_{n}*R*given the constraint that this line should pass through the point^{i}_{n}*R*̄_{n,1}:$${u}_{n}=\mathrm{Im}{\overline{R}}_{n,1}-\mathit{v}\phantom{\rule{.2em}{0ex}}\mathrm{Re}{\overline{R}}_{n,1},$$and the slope

*v*is estimated by least squares estimation:$${v}_{n}=\frac{{S}_{\mathit{RI}}-{S}_{R}\mathrm{Im}{\overline{R}}_{n,1}-{S}_{I}\mathrm{Re}{\overline{R}}_{n,1}+{S}_{w}\mathrm{Re}{\overline{R}}_{n,1}\mathrm{Im}{\overline{R}}_{n,1}}{{S}_{\mathit{RR}}-2{S}_{R}\mathrm{Re}{\overline{R}}_{n,1}+{S}_{w}{\left(\mathrm{Re}{\overline{R}}_{n,1}\right)}^{2}}$$where

*S*= ∑_{R}Re_{i}w^{i}*R*,^{i}_{n}*S*= ∑_{I}Im_{i}w^{i}*R*,^{i}_{n}*S*= ∑_{RR}(Re_{i}w^{i}*R*)^{i}_{n}^{2},*S*= ∑_{RI}Re_{i}w^{i}*R*Im^{i}_{n}*R*, and^{i}_{n}*S*= ∑_{w}. The weights_{i}w^{i}*w*are chosen equal to 1/^{i}*I*, where^{i}*I*is the total number of counts in pixel^{i}*i*. The calculation of these sums does not require all*R*to be stored in memory simultaneously, and memory requirements are therefore low.^{i}_{n}

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.

## 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*
_{0} 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.

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.

#### 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*
_{0} 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*
_{1} calculated from Eq. (14) corresponded well to the value obtained by Fourier transformation of the measured IRF (0.895 + 0.446*j* vs. 0.892 + 0.450*j*). 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)).

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 *τ*
_{1} = 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.

## Acknowledgment

We thank Philippe I. H. Bastiaens for useful comments on the manuscript.

## References and links

**1. **J. R. Lakowicz, *Principles of Fluorescence Spectroscopy*, 3rd ed. (Springer, 2006). [CrossRef]

**2. **K. Carlsson, A. Liljeborg, R. M. Andersson, and H. Brismar, “Confocal pH imaging of microscopic specimens using fluorescence lifetimes and phase fluorometry: influence of parameter choice on system performance,” J. Microsc. **199**, 106–114 (2000). [CrossRef] [PubMed]

**3. **R. M. Clegg, “Fluorescence resonance energy tranfer,” Fluorescence Imaging Spectroscopy and Microscopy **137**, 179–251 (1996).

**4. **P. I. H. Bastiaens and A. Squire, “Fluorescence lifetime imaging microscopy: spatial resolution of biochemical processes in the cell,” Trends Cell Biol. **9**, 48–52 (1999).' [CrossRef] [PubMed]

**5. **F. S. Wouters, P. J. Verveer, and P. I. H. Bastiaens, “Imaging biochemistry inside cells,” Trends Cell Biol. **11**, 203–211 (2001). [CrossRef] [PubMed]

**6. **A. Schönle, M. Glatz, and S. W. Hell, “Four-dimensional multiphoton microscopy with time-correlated single-photon counting,” Appl. Opt. **39**, 6306–6311 (2000). [CrossRef]

**7. **W. Becker, A. Bergmann, M. A. Hink, K. König, K. Benndorf, and C. Biskup, “Fluorescence lifetime imaging by time-correlated single-photon counting,” Microsc. Res. Tech. **63**, 58–66 (2004). [CrossRef]

**8. **M. Peter and S. M. Ameer-Beg, “Imaging molecular interactions by multiphoton FLIM,” Biol. Cell **96**, 231–236 (2004). [CrossRef] [PubMed]

**9. **E. Gratton, S. Breusegem, J. Sutin, Q. Ruan, and N. Barry, “Fluorescence lifetime imaging for the two-photon microscope: time-domain and frequency-domain methods,” J. Biomed. Opt. **8**, 381–390 (2003). [CrossRef] [PubMed]

**10. **S. Pelet, M. J. R. Previte, L. H. Laiho, and P. T. C. So, “A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation,” Biophys. J. **87**, 2807–2817 (2004). [CrossRef] [PubMed]

**11. **P. Barber, S. Ameer-Beg, J. Gilbey, R. J. Edens, I. Ezike, and B. Vojnovic, “Global and pixel kinetic data analysis for FRET detection by multi-photon time-domain FLIM,” Proc. SPIE **5700**, 171–181 (2005). [CrossRef]

**12. **P. Barber, S. Ameer-Beg, J. Gilbey, L. Carlin, M. Keppler, T. Ng, and B. Vojnovic, “Multiphoton time-domain fluorescence lifetime imaging microscopy: practical application to protein-protein interactions using global analysis,” J. R. Soc. Interface **6**, S93–S105 (2008). [CrossRef]

**13. **P. J. Verveer, A. Squire, and P. I. H. Bastiaens, “Global analysis of fluorescence lifetime imaging microscopy data,” Biophys. J. **78**, 2127–2137 (2000). [CrossRef] [PubMed]

**14. **P. J. Verveer, F. S. Wouters, A. R. Reynolds, and P. I. H. Bastiaens, “Quantitative imaging of lateral ErbB1 receptor signal propagation in the plasma membrane,” Science **290**, 1567–1570 (2000). [CrossRef] [PubMed]

**15. **P. J. Verveer and P. I. H. Bastiaens, “Evaluation of global analysis algorithms for single frequency fluorescence lifetime imaging microscopy data,” J. Microsc. **209**, 1–7 (2003). [CrossRef] [PubMed]

**16. **A. H. A. Clayton, Q. S. Hanley, and P. J. Verveer, “Graphical representation and multicomponent analysis of single-frequency fluorescence lifetime imaging microscopy data,” J. Microsc. **213**, 1–5 (2004). [CrossRef]

**17. **T. Ng, M. Parsons, W. E. Hughes, J. Monypenny, D. Zicha, A. Gautreau, M. Arpin, S. Gschmeissner, P. J. Verveer, P. I. H. Bastiaens, and P. J. Parker, “Ezrin is a downstream effector of trafficking PKC-integrin complexes involved in the control of cell motility,” EMBO J. **20**, 2723–2741 (2001). [CrossRef] [PubMed]

**18. **A. R. Reynolds, C. Tischer, P. J. Verveer, O. Rocks, and P. I. H. Bastiaens, “EGFR activation coupled to inhibition of tyrosine phosphatases causes lateral signal propagation,” Nat. Cell Biol. **5**, 447–453 (2003). [CrossRef] [PubMed]

**19. **O. Rocks, A. Peyker, M. Kahms, P. J. Verveer, C. Koerner, M. Lumbierres, J. Kuhlmann, H. Waldmann, A. Wittinghofer, and P. I. H. Bastiaens, “An acylation cycle regulates localization and activity of palmitoylated Ras isoforms,” Science **307**, 1746–1752 (2005). [CrossRef] [PubMed]

**20. **G. Xouri, A. Squire, M. Dimaki, B. Geverts, P. J. Verveer, S. Taraviras, H. Nishitani, A. B. Houtsmuller, P. I. H. Bastiaens, and Z. Lygerou, “Cdt1 associates dynamically with chromatin throughout G1 and recruits Geminin onto chromatin,” EMBO J. **26**, 1303–14 (2007). [CrossRef] [PubMed]

**21. **M. Digman, V. R. Caiolfa, M. Zamai, and E. Gratton, “The Phasor approach to fluorescence lifetime imaging analysis,” Biophys. J. (2007). [PubMed]

**22. **R. A. Colyer, C. Lee, and E. Gratton, “A novel fluorescence lifetime imaging system that optimizes photon efficiency,” Microsc. Res. Tech. **71**, 201–213 (2008). [CrossRef]

**23. **S. Padilla-Parra, N. Audugé, M. Coppey-Moisan, and M. Tramier, “Quantitative FRET analysis by fast acquisition time domain FLIM at high spatial resolution in living cells,” Biophys. J. **95**, 2976–2988 (2008). [CrossRef] [PubMed]

**24. **D. M. Jameson, E. Gratton, and R. Hall, “The measurement and analysis of heterogeneous emissions by multi-frequency phase and modulation fluorometry.” Appl. Spec. Rev. **20**, 55–106 (1984). [CrossRef]

**25. **A. Esposito, H. C. Gerritsen, and F. S. Wouters, “Fluorescence lifetime heterogeneity resolution in the frequency domain by lifetime moments analysis,” Biophys. J. **89**, 4286–4299 (2005). [CrossRef] [PubMed]

**26. **G. I. Redford and R. M. Clegg, “Polar plot representation for frequency-domain analysis of fluorescence lifetimes,” J. Fluoresc. **15**, 805–815 (2005). [CrossRef] [PubMed]

**27. **G. W. Gordon, G. Berry, X. H. Liang, B. Levine, and B. Herman, “Quantitative fluorescence resonance energy transfer measurements using fluorescence microscopy,” Biophys. J. **74**, 2702–2713 (1998). [CrossRef] [PubMed]

**28. **F. S. Wouters and P. I. H. Bastiaens, “Fluorescence lifetime imaging of receptor tyrosine kinase activity in cells,” Curr. Biol. **9**, 1127–1130 (1999). [CrossRef] [PubMed]