A novel fluorescence lifetime imaging system that optimizes photon efficiency.

Fluorescence lifetime imaging (FLIM) is a powerful microscopy technique for providing contrast of biological and other systems by differences in molecular species or their environments. However, the cost of equipment and the complexity of data analysis have limited the application of FLIM. We present a mathematical model and physical implementation for a low cost digital frequency domain FLIM (DFD-FLIM) system, which can provide lifetime resolution with quality comparable to time-correlated single photon counting methods. Our implementation provides data natively in the form of phasors. On the basis of the mathematical model, we present an error analysis that shows the precise parameters for maximizing the quality of lifetime acquisition, as well as data to support this conclusion. The hardware and software of the proposed DFD-FLIM method simplifies the process of data acquisition for FLIM, presents a new interface for data display and interpretation, and optimizes the accuracy of lifetime determination.


INTRODUCTION
We have developed a physical implementation and statistical model of a new method for fluorescence lifetime imaging (FLIM) data collection and analysis. Our approach falls in the category of the ''frequency domain'' approach to lifetime acquisition, yet uses a detector operating in the photon counting mode. This digital frequency domain (DFD) method overcomes the problems of duty cycle, modulation of the detector gain and expensive radio frequency synthesizers used in the classical analog frequency domain approach. In our approach we implemented all the operations performed in a frequency-domain lifetime instrument in a digital form using a single field programmable chip. Since all operations including the generation of the light modulation frequency, the generation of the cross-correlation sampling frequency and the assignment of the time of arrival of a photon to a bin are digital, there are no calibrations or adjustments to be performed. The mathematical model presented below fully accounts for all the elements of the DFD method. In addition, the mathematical model reproduces, as a limiting case, the principle of the time-correlated single photon counting (TCSPC) approach. Therefore, on a common statistical basis we can compare the two approaches and derive some general conclusions about the relative precision of the two methods. We found that with proper system design the two methods can be made to have comparable precision. More importantly, the mathematical model was used to maximize the precision of the DFD implementation and to determine which parameters are crucial to reach optimal performance. In terms of precision of the lifetime measurement, we were able to fully quantify the effect of the instrument response including the jitter of the detection system.
Fluorescence lifetime is a fundamental spectroscopic quantity that allows quantitative analysis through several approaches, including the identification of molecu-lar species based on lifetime, Förster Resonance Energy Transfer (FRET), contrast due to different ion concentrations, and measurements of chemical equilibria. For measurements done in a cuvette, in which spatial resolution is not needed, there are two major approaches for the acquisition of the fluorescence decay. One is based on TCSPC or time-sampling of the intensity decay after pulse excitation, and a second is based on the measurement of the harmonic response of the fluorescence system.
When spatial resolution is needed, such as with microscopy, different considerations come into play depending on the kind of microscope used. One major difference between the laser scanning confocal microscope and the camera-based microscopes is that in the former the detector works in the serial mode, although there are some recent scanning instruments using multiple foci in conjunction with a camera to collect the image (Grant et al., 2005;van Munster et al., 2007). For FLIM instruments operating in the serial mode, a bottleneck in the rate of data acquisition is caused by the recovery time of the TAC element that is common to the TCSPC-based instruments.
During the last few years, several new instruments were introduced for FLIM operation in the confocal microscope by several manufacturers. As a consequence, FLIM has become more common and is now available in several labs. One of the reasons for the interest in FLIM is that FRET imaging by lifetime resolved methods is generally considered to have fewer problems with background and autofluorescence than intensity ratio methods. However, FLIM methods are limited to few labs due to the high cost of the instrumentation and the difficulty of performing lifetime analysis in many pixels simultaneously (Pelet et al., 2004). The analog frequency-domain approach offers some simplification in the analysis methods and in the laser sources used (Clegg et al., 1992). However, the traditional frequency-domain electronics operating in the radio frequency range require gain modulated detectors, radio frequency amplifiers (Gratton and Limkeman, 1983), and are not a simple addition to existing laser confocal microscopes. Whatever approach is utilized, current methods are relatively expensive, require specialized electronic and modulated sources, and involve sophisticated analysis methods to extract information about the lifetime image and the processes that produce lifetime variations in different pixels of an image (Pelet et al., 2004).
In this article we describe new data acquisition hardware that requires minimal modifications to the configuration of common commercial laser confocal microscopes. The cost of the new electronics is minimal. We also use the phasor method of data analysis that is ''native'' to the proposed hardware and simplifies the calculation and presentation of lifetime images. Overall, the proposed approach has the potential to make FLIM technology more widely available. Our approach uses serial detectors in the photon counting mode, and the digital heterodyning method to acquire data which is directly analyzed in the frequency domain. The principle of the digital heterodyning method was previously described (Eid, 2002). However, the original implementation by Eid (2002) used a specialized acquisition card and had limitations in terms of the duty cycle and speed of data acquisition. In this article we develop a mathematical model of DFD-FLIM. Using this model we were able to choose optimal system parameters, which minimize the distribution of lifetime values. Using an FPGA chip, we have implemented a version of the digital heterodyning method, which has a 100% duty cycle so that no photons are lost, operates in several harmonic frequencies simultaneously, can be used in conjunction with common detectors in commercial laser scanning microscopes, only requires a modulated light source instead of a pulsed source, has uncertainty levels comparable with TCSPC methods, and has very high throughput with a very low cost. The cost of the FPGA chip and evaluation board is below $100. The mathematical model contains, as a limiting case, the description of data acquisition using the TCSPC method of data acquisition. Using the model, we were able to compare theoretically and experimentally the proposed DFD and the TCSPC methods and to conclude that they have the same statistical accuracy. We were also able to predict the effect of instrument jitter on the precision of the DFD approach and to conclude that there is an optimal jitter level that improves the precision of the determination of very short lifetimes.
We implemented the DFD approach in several confocal microscopes, both with 1-photon and 2-photon excitation, and with home built and commercial instruments. As shown below, the hardware collects data directly under the form of ''phasors'' at several harmonic frequencies simultaneously. The phasor representation of the fluorescence decay allows a very simple interpretation of the FLIM image and the calculation of FRET efficiencies without the usual translation of the decay into exponential components. (The calculation of FRET efficiencies will be described in a forthcoming publication). The phasor approach to data analysis was previously proposed by us and by others (Clayton et al., 2004;Jameson et al., 1984;Redford and Clegg, 2005). It provides a simple graphical interface for FLIM data presentation and analysis without the need of fitting the fluorescence decay at each pixel. We compared data collected with the DFD and the TCSPC method on the same sample. After optimization of the design parameters, the DFD approach gave results comparable to those obtained with TCSPC.

MATERIALS AND METHODS Introduction to Digital Frequency-Domain
Methods Frequency domain lifetime methods are characterized by the usage of a periodic modulated excitation, for which the finite lifetime of the fluorophore results in a phase delay and demodulates the emission. These characteristics of the response at any given modulation frequency are traditionally represented by s p and s m , respectively. A derivation of these equations is given by (Spencer and Weber, 1969): Traditional analog approaches to frequency domain acquisition use modulation of the gain of a photomultiplier tube (PMT) or of an image intensifier. For DFD lifetime acquisition, we replace the modulation of these analog signals with a digital modulation scheme inside the data collection circuitry (rather than at the detector level) by using multiple sampling windows within each excitation period, as shown in Figure 1. Because the modulation is performed digitally, the signal can be processed by an arbitrary number of sampling windows without any loss of signal, and the PMT can be operated at full gain during the entire acquisition.
To achieve sensitivity to short lifetimes and to ensure evenly distributed sampling of the fluorescence emission, we use heterodyning between the frequency of the sampling windows and the excitation frequency of the light source, as shown in Figure 1. As will be described later, this allows us to translate the fluorescence response into a cross-correlation phase histogram, as shown in Figure 2, which contains a functional form given by the convolution of the fluorescence emission with the shape of the sampling window.

DFD Hardware
We created a device called the FLIMBox, in which the DFD algorithm is implemented in a system which uses a field programmable gate array (FPGA). An FPGA is a chip which can be rapidly reprogrammed with different circuit layouts by uploading firmware, and thus serves as a convenient tool for the development of scientific hardware. The specific FPGA we used was a Xilinx 1 Spartan 1 3E, XC3S100E TM (San Jose, CA), on an Avnet Electronics Marketing Evaluation Kit (Phoenix, AZ) with a Cypress EZ-USB FX2 TM chip (San Jose, CA).
The FPGA contains two digital clock managers that provide clock synchronization to an external clock and clock multiplication services. Each digital clock manager multiplies an input clock frequency by n d /m d , where n d and m d are integers ranging from 1 to 32. To implement the heterodyning principle in a digital system, it is convenient to have a cross-correlation frequency, which is a whole integer fraction of the sampling frequency, f s . To reach the minimal uncertainty plateau described in the parameters optimization section below, at least four sampling windows are required. This can be obtained by choosing n 1 5 32, m 1 5 17, n 2 5 32, m 2 5 15, and dividing the frequency by four, such that a clock is generated, which is four times the sampling clock. The resulting cross-correlation frequency of f cc 5 f s /256, is given by: For example, for a 48 MHz excitation frequency, a 43 clock is generated at 192.8 MHz, yielding a sampling frequency of 48.188 MHz and a cross-correlation frequency (the difference between the light modulation and the sampling frequency) of 188 kHz. This fast 43 clock is then used as the input for a counter which tags incoming photons with the sampling window number (0-3) which corresponds to their arrival time.
To relate the window during which a photon arrived to a portion of the excitation period, it is necessary to know the relative phase difference between the sampling and excitation clocks. For digital stability of this phase difference measurement, this was implemented as a cross-correlation phase counter, which counts at the time scale of the sampling clock, has a periodicity equal to the cross-correlation frequency, and uses negative feedback on the measurement of the excitation clock to lock on to a consistent phase relationship. This circuit enforces the same phase relationship between the cross-correlation counter and the sampling and excitation clocks each time the device is activated. Therefore the cross-correlation phase counter provides a consistent measurement of the relative phase difference between the sampling and excitation clocks.
The circuit then outputs a value identifying the arrival window, w arrival , and the cross-correlation counter value, p ccc , for each photon count. For our implementation, these two values can be combined to form a crosscorrelation phase index as follows: where n w is the total number of sampling windows. This cross-correlation phase index is used to construct the cross-correlation phase histogram, H(p), which is a histogram of the phase indexes for each photon detected.
A requirement for imaging in a raster-scan instrument is a mechanism for scanner synchronization, so that the system can ensure that the data for each pixel is acquired at the same physical position without significant drift. A scan-enabled control line was included in the chip, which enables data collection only when it is active, allowing the scanner control mechanism to signal when each frame (or line) has started. We then use the total time of arrival of a photon with respect to the starting of the frame (macro-time) to divide the information into pixels, each of which has a corresponding phase histogram.

System Layout
The complete schematic layout can be seen in Figure 3, where an example is shown with the FLIMBox installed on a commercial confocal system with a modulated diode laser. We used an Olympus FluoView FV1000 (Olympus America, Center Valley, PA) with a variety of modulated diode lasers from ISS (Champaign, IL). The diode laser is driven by a 48 MHz LVTTL signal produced by the FLIMBox which is high for 458 of the repetition period. In other system configurations we have also used 80 MHz titanium-sapphire lasers in place of a modulated diode laser.
To provide an optical phase reference, the laser output is then split by a beam splitter, which deflects a portion of the light to a Gigahertz photodiode (Det200, Thorlabs, Newton, NJ). This photodiode signal is then amplified and sent to a zero-crossing trigger which produces a LVTTL signal that goes into the FLIMBox and serves as a frequency and phase reference for f ex . The main portion of the laser is sent into the microscope. Without this optical phase reference, we found that the laser can sometimes introduce an unpredictable and significant phase shift, which disrupts phase accuracy. With the optical phase reference in place, the system was stable to within 0.18 in phase.
The outputs of the PMTs for two channels are sent to Gigahertz amplifiers (ACA-4-35-N, Becker & Hickl, Berlin, Germany), which are then connected to constant fraction discriminators (Model 6915, Phillips Scientific, Mahwah, NJ) which trigger an LVTTL signal on the zero-slope of the PMT response when it goes past a set threshold.
The photon count LVTTL signals are processed inside of the FLIMBox on two fully independent channels. The photon arrival time information is then placed into a FIFO, from which the data is transferred via USB to a computer for processing.
The cross-correlation phase histogram at each pixel, H(p), is then constructed by a computer program and used for the FLIM analysis. The intensity image can be obtained by simply summing the points (n p ) of the phase histogram at each pixel as follows: HðpÞ: Phasor Calculation In this work we use the phasor representation of the fluorescence decay for FLIM data analysis. The phasor approach was originally proposed by us (Jameson et al., 1984) and subsequently expanded by others (Clayton et al., 2004;Redford and Clegg, 2005). To analyze data with the phasor approach, it is necessary to calculate the intensity-normalized phasor components for harmonic h as follows: HðpÞ cosð2php=n p Þ; ð6Þ HðpÞ sinð2php=n p Þ: From these two terms, the phase and modulation values are calculated according to the vector transformation equations: The instrument response function is accounted for by scaling and rotation of the phasor using a reference phasor from a known reference lifetime. The modulation and phase values of the known lifetime are inserted into Eqs. (24) and (25) to obtain the modulation factor and phase shift of the instrument response at each harmonic of interest. This allows one to completely account for the response of the system with a single measurement.
A global representation of FLIM images at any given harmonic frequency can be obtained by turning / F and m F back into instrument response corrected versions of the cosine and sine values of the Fourier transform.
Each pair of g F,h and s F,h values can be treated as a vector called a phasor, and these values can be accumulated on a two-dimensional histogram called a phasor plot, as shown in Figure 4a.
Presenting FLIM data in this format provides a number of advantages. Since the g and s coordinates come from Eqs. (6) and (7), they are both intensity normalized linear coordinates. The application of the instrument response only performs a scaling and a rotation, so this does not disrupt the linearity. Since g and s are linear, this means that the phasor positions of multi-exponential lifetimes are the vector sum of their single-exponential phasors. This also means that each pixel with multiple molecular species has a phasor position which is a vector sum of the phasors for each species, according to: where f i,h are the fractional fluorescence intensities of each phasor component such that: With knowledge of the relative brightness of two species, this property of the phasors allows one to extract relative concentrations with ease.
Using Eqs. (1) and (2) we can see that all single-exponential lifetimes occur when s p 5 s m , and thus when m F,h 5 cos(/ F,h ) (Jameson et al., 1984). This corresponds to a semicircle on the phasor plot called the ''universal circle'' which is given by: All multiexponential lifetimes are therefore linear combinations of single-exponential lifetimes, as given by Eqs. (12)-(14), and therefore all multiexponential lifetimes have phasors, which are inside the universal circle.

DFD-FLIM Theory
To maximize the ability of DFD-FLIM hardware to determine lifetimes, we must reduce the uncertainty spread of phasors in the phasor plot. In Figure 4, where the plot is an experimental histogram of phasors at each pixel, a decrease in uncertainty would correspond to more pixels having phasors closer to the center of the phasor distribution, or a decrease in the uncertainty of the phase and modulation. To achieve this, we developed a theoretical model of the data acquisition  and processing, which shows which experimental parameters are important and how to change these parameters to minimize the width of the phasor distribution.
Consider an arbitrary normalized fluorescence excitation function E(t) modulated at a frequency f ex , which can be written as a Fourier series, where m E,h is the modulation of the excitation at each harmonic, and / E,h is an arbitrary phase shift. The probability distribution for the fluorescence emission is given by the convolution between the excitation function and the fluorescence response of the fluorophores.
For an arbitrary combination of fluorophores, the fluorescence response function can be written as: where s i are the fluorescence lifetimes, f i are the fractional intensities, and f 0 is the contribution of the uncorrelated background. This emission is then detected by a PMT (or other photon counting device), a discriminator, and a logic gate, with a total time response jitter which can be approximated as a Gaussian: where s j is the standard deviation of the detection time jitter, composed of a transit time spread, discriminator jitter, and triggering jitter for the digital logic. The arrival time is then resolved within the circuit into one of n w arrival time windows, each of which is shaped as a periodic boxcar function.
where w arrival is the index of the arrival time window being considered (ranging from 0 to n w 2 1), and f s is the sampling frequency. Since this sampling process is performed digitally, every photon is counted in one of the sampling windows, resulting in a 100% duty cycle. We note that if the number of sampling windows is very large, this model corresponds to the data produced by the TCSPC method. However, with TCSPC the assignment of the photon to a given bin is obtained by measuring the time delay between emission and excitation by the TAC method (time-to-amplitude converter).
In the DFD method the assignment is done by tagging the photon with a number equal to the phase shift between the sampling and excitation clocks. We calculate the probability distribution of the detected photons by convoluting the excitation, fluorescence response, jitter of the system, and sampling window: This equation can be greatly simplified by using the convolution theorem, so that: Since E(t) and S w,arrival (t) have two different frequencies, the frequency of H w,arrival (t) is actually the crosscorrelation frequency given by (Spencer and Weber, 1969): We can see that each arrival time window represented by S w,arrival (t) will result in an identical probability distribution in H w,arrival (t), but with a phase offset w arrival /(n w f cc ). This allows us to recombine the separate sampling windows during data processing into a single probability distribution H(t). Because our analysis technique specifically utilizes the harmonics of the Fourier series, we can perform the analysis on the harmonics of H(t) given by: The resulting modulation of H(t) for each harmonic h is given by the product of all the component modulations, as: Therefore, if we know the modulation of the instrument response, m IR,h 5 m E,h m J,h m S,h , we can extract the modulation of the fluorescence lifetime response, m F,h .
The resulting phase of H(t) for each harmonic corresponds to the relative phase difference between the sampling period and the excitation period, plus the phase offset provided by the fluorescence lifetime response, as follows: Therefore, if we know the relative phase difference between the sampling period and the excitation period, / IR,h 5 (/ S,h 2 / E,h ), then we can extract the phase off-set provided by the fluorescence lifetime response. Both the modulation and the phase of the instrument response can be obtained by measuring a reference sample with a known lifetime value.
Derivations similar to this section have been previously given, such as in Jameson et al. (1984). However, the digital sampling windows and jitter terms were not included in previous derivations.

Optimization of the DFD Parameters
To determine the uncertainty in phasor values, we must find the standard deviation of the measured phase and modulation for each harmonic. We can perform the analysis on the phase histogram H(p), which represents one period of H(t), and where 0 p n p 2 1. To perform this error analysis, we first plot the crosscorrelation phase histogram in polar coordinates for the harmonic of interest. For the first harmonic, this is a plot with H(p) as the radial value, 2pp/n p as the phase value, and where p 5 n p represents a complete orbit. An example is given in Figure 5. In this coordinate system, the Fourier transform for each harmonic can be viewed graphically as the vector sum of the values plotted. M h 5 m H N and / H are the magnitude and angle of this vector sum. From this graphical representation we can observe the effect of fluctuations in H(p) on the values of m H and u H . By the principles of vector addition, u H will only be affected by fluctuations in H(p), which occur perpendicularly to the vector sum. The variances introduced by each of these components at every value of p can be added linearly. So if we consider the Poissonian error in each bin of H(p) (this is the phase histogram in which every bin contains independent events), we obtain the standard deviation in phase as follows: The derivation of the expression for the standard deviation of modulation is similar, except using the fluctuations which are parallel to the vector sum, and taking into account that modulation is normalized by N, and thus fluctuations will also affect m H via the contributions of N. These uncertainty values on phase and modulation can then be propagated through Eqs. (1) and (2) to obtain the uncertainties on s p and s m : From this analysis we can learn several insightful properties about the behavior of the phasor distribution. Firstly, the distribution of both phase and modulation decreases by the square root of the number of counts in the phase histogram acquired at each pixel. Secondly, the distribution of both phase and modulation is inversely proportional to m IR,h , the modulation of the instrument response for each harmonic. Since m IR,h ranges from 0 to 1, the value which produces the optimal distribution is when it is equal to 1. From consideration of m IR,h 5 m E,h m J,h m S,h we can then immediately see that the optimal configuration places each of those component modulation values close to 1. For a square wave excitation of width y, the modulation values as shown in Figure 6a are found by taking the Fourier transform of that square wave, which yields: Similarly, for the sampling window, as seen in Figure  6b, the modulation values are given by: For Gaussian system jitter, as seen in Figure 6c, the modulation values are given by: This analysis allows us to quantitatively evaluate the design parameters required to minimize the distribution of phasors. For lifetime values obtained from the first harmonic, we can see that the modulation, and thus the information content, has reached a plateau with an excitation square wave of 908 or smaller (Fig.  6a), four or more sampling windows (Fig. 6b), and system jitter less than 2 ns (Fig. 6c). Improvements beyond this point will make only marginal improvements to the statistical uncertainties of a measurement with the first harmonic. In our specific system configuration with a 458 square wave excitation, this produces an m ex,1 of 0.9745, an m ex,2 of 0.9003, and an m ex,3 of 0.7842 at the three lowest harmonics. So even though a diode laser driven in this manner has a pulse width as wide as 2.6 ns, it will still produce results of equivalent quality to a femtosecond pulsed laser for the first two harmonics.

Multiexponential Analysis
An extensive amount has been written about the extraction of component lifetimes by multi-frequency fitting Jameson and Gratton, 1983;Jameson et al., 1984;Lakowicz et al., 1984). In our imaging system used in FLIM we only excite at a single laser repetition frequency with a pulse sufficiently narrow to contain several harmonic frequencies. Each harmonic yields two independent variables, the phase and modulation, so a direct relationship can be established between the number of harmonics and the number of lifetime components. For example, using two harmonics, two lifetime components can be calculated using an exact formula (two lifetime values, one fractional component and one background term) using the principles outlined by (Weber, 1981).

Solution Measurements
We evaluated the FLIMBox hardware, phasor analysis, and the above error analysis by performing solution lifetime measurements on a FLIM setup. We prepared fluorescein in a pH 10 solution, which is known to have a lifetime of 4.05 ns as determined in cuvette measurements. Then we used this fluorescein solution as our lifetime reference for determining the instrument response, and also as a sample for examining the pixel uncertainty. We then prepared a solution of rhodamine B in water as an additional sample. Figure 4a shows under the form of an image, the histogram of the phasor values for a 256 3 256 pixel image of fluorescein taken at 48 MHz, excited with a 470 nm diode laser, in the Olympus FV1000 with FLIMBox described earlier, with that same fluorescein sample used as the reference. The measurement was acquired at 200,000 counts per second, with 40 ls per pixel. Pixel phase histograms were accumulated until there were on the average of 538 counts in each pixel. Figure 4b shows the rhodamine B phasor plot taken on the same instrument, but with 50,000 counts per second and an average of 202 counts per pixel.
By using the error Eqs. (28) and (31), we calculated the expected uncertainties in both phase and modulation, and we compared these with the experimental standard deviations for the fluorescein and rhodamine B measurements. We found very high correspondence between the theory and experimental results, as shown in Table 1. We repeated this comparison with many measurements under various conditions and system configurations, including for data taken with a TCSPC system, and found the equations holding in each case. This confirms that the theoretical uncertainty equations completely account for the precision of the system, and confirms that under real physical conditions, the statistical errors of a FLIM measurement can be correctly evaluated and effectively improved by the optimized parameters derived by the mathematical model of the DFD method.
We also evaluated the phase and modulation lifetime values and uncertainties of the two measurements. Since the fluorescein measurement is used as a reference it is not meaningful to consider its lifetime value, since these are fixed at 4.05 ns by the referencing process. However, we can still consider the uncertainties in these. The fluorescein measurement had a s p per pixel of 4.05 6 0.522 ns, and a s m per pixel of 4.05 6 0.379 ns. The predicted values for the uncertainties according to Eqs. (32) and (33) are 0.505 and 0.376 ns. To avoid skew in the lifetime values, the mean of all lifetime values is not used to determine rhodamine B's mean lifetime. Instead, the lifetime of the center of mass of all the phasors is used, as the phasor addition is linear. For the rhodamine B measurement, the measured s p for each pixel was 1.723 6 0.287 ns, and the measured s m was 2.001 6 0.448 ns. The predicted uncertainties per pixel are 0.281 and 0.428 ns.
Note that the phasor for rhodamine B is not exactly on the universal circle, although we expected a single exponential according to our measurement in cuvette. The phasor analysis provides an immediate and intuitive explanation for this result. In this case, the system on which the measurement was taken had background counts of 1500 per second out of a total 50,000 counts per second, which corresponds to a 3% background for the rhodamine B measurement. Since uncorrelated background counts occur at a phasor position at the (0,0) coordinate, this phasor is linearly combined according to Eqs. (12)-(14) with the phasor for rhodamine B's actual lifetime. This linear combination with the uncorrelated background reduces m and therefore increases s m , but not the s p value. When the phasor contribution of the uncorrelated background is subtracted from the rhodamine B phasor, the phasor coordinate will move to the universal circle, and s m becomes 1.7 ns.
If instead of examining the above measurements at each pixel, we use the same data to evaluate the uncertainty of the mean, rather than the uncertainty of each pixel, we find much more statistical precision. The fluorescein measurement then has a phasor point which is determined by 35 million counts, and the rhodamine measurement has a phasor point which is determined by 13 million counts. The uncertainties of the means for the fluorescein lifetimes are then 2.0 ps and 1.5 ps for s p and s m . Similarly for rhodamine B the uncertainties of the means are 1.1 ps and 1.8 ps for s p and s m . This emphasizes that the precision of lifetime measurement is not constrained by the precision of individual photon arrival times. The precision of the average value is determined by statistics, and therefore is a function of the system parameters and the number of counts. The accuracy of the phasor position determination depends on the stability of the system after it has been referenced with a known lifetime.
To evaluate the long-term stability of our system, we examined the trend of phase and modulation values over the course of hours and days using an excitation modulation of 48 MHz. We determined that after an initial 45 min warm-up period, the phase remained stable within around 0.18, and the modulation remained stable within less than 0.001 over time periods of hours to days. These correspond to lifetime stability around 10 ps for fluorescein, as given by Eqs. (32) and (33).
For comparison with an existing standard for lifetime measurement hardware, we measured fluorescein with a 443 nm diode laser modulated at 48MHz and the FLIMBox on our FV1000 setup, and with a Becker & Hickl TCSPC card (model 830, Becker & Hickl, Berlin, Germany) connected to a home-built 2-photon system using an 80 MHz Titanium-Sapphire laser. According to our model, when the data are processed in the same way and for the same number of counts, the precision of the determination of the phase and modulation values should only depend on the instrument modulation factor. Since at 80 MHz, the modulation factor for fluorescein is less than at 48 MHz, we anticipated that the precision of the data acquired with FLIMBox (at 48 MHz) should be better than the precision of data acquired with the B&H card (at 80 MHz). As can be seen in Figure 7, the phase uncertainties obtained for the measurement of fluorescein with the FLIMBox is less than that of the TCSPC system under these two sets of conditions. The relevant point here is the demonstration that both data acquisition systems reach the statistical limits in the precision of the measurement of lifetime values, and that the modulation factor and the number of photons fully account for the precision of the measurement. Gerritsen et al. (2006) define a count-independent metric for the precision of a lifetime measurement as follows: By inserting Eqs. (32) and (33) in Eq. (37), we are able to obtain the minimum F value (representing the highest sensitivity) for the highest modulation of the instrument response, as shown in Figure 8.

Multiexponential Analysis of Enhanced Cyan
Fluorescent Protein We measured the lifetime of enhanced cyan fluorescent protein in a 20 mM Tris buffer solution with a 443 nm diode laser for excitation. Using the first and second harmonic to determine two exponential components and the amount of unmodulated background (we use an exact formula since we have two phases and two modulations and four unknown parameters), we obtained a relative contribution of (12 6 1)% for a lifetime of 0.75 6 0.03 ns, and 88% for a lifetime of 3.17 6 0.03 ns. The standard deviations of the lifetime determinations are obtained from eight independent measurements of this system. We also noticed a slight spectral dependence of the relative contributions. These results correspond well with the literature, confirming the ability to extract multiple exponential components by using DFD-FLIM (Tramier et al., 2004).

Fluorescence Lifetime Images
To demonstrate the ability to acquire lifetime images, and the linear combination of phasors we obtained a series of letter-shaped microchannels with a fluorescein solution in the letter ''L,'' a rhodamine B solution in the letter ''D,'' and a mixture of the two solutions in the letter ''F.'' Figure 9a shows the lifetime image generated from s p , and Figure 9b shows the phasor plot obtained from this FLIM image. Taking advantage of the intrinsic spatial correlation in images, the phasors for this image were resampled by averaging with their nearest neighbors to improve statistical accuracy.
By using a graphical calculator inside of the Globals for Images program (LFD, Irvine, CA) which applies Eqs. (12)-(14), we calculated that the ''F'' point on the phasor plot is 25% away from the ''L'' phasor, and 75% away from the ''D'' phasor. This shows that the phasor for the letter ''F'' is composed of 75% of the fluorescein solution by relative intensity.

DISCUSSION
We have developed a mathematical framework that provides an understanding of the effect of different parameters in the DFD-FLIM approach. We have shown that the experimental realization of the DFD principle in the FPGA chip gives values in very good accord with the model. This implies that the there are no additional factors (systematic or random) in the hardware implementation which are not accounted for by the mathematical model. As a consequence of the theoretical description, we have determined the system parameters which optimize the precision of phasor distributions and corresponding lifetime values. We have reached these optimal parameters in our implementation of the FLIMBox for a hypothetical system in which the lifetimes values range between 0.1 and 10 ns. We have also shown that using an optical reference for the excitation, the drifts of the physical implementation are very small, giving better than 10 ps stability for the measurement of the lifetime of fluorescein at 48 MHz.
In this work we provide a statistical framework that correctly predicts the errors in lifetime determination for any given sample and system configuration. Equations (28) and (31) show that the best precision of the position of the phasor (and of the lifetime value) is achieved when all the modulation factors are close to 1. This corresponds to the well-known principle that given a certain lifetime value, the best precision is obtained when the excitation pulse is narrow and repeats with a period, which is much longer than the lifetime, the time windows (of the photon delay histogram) are in a large number and the instrument jitter is small. However, the value of this equation is that it allows us to quantitatively predict how much each of these factors will ultimately affect the precision of the lifetime determination. For example, we can predict that for a lifetime of about 3 ns (typical of a GFP) quasi optimal measurement conditions are obtained at excitation frequencies around 40-50 MHz with pulses of about 1-2 ns, using four time windows in the DFD implementation and with instrument jitter on the order of 1-2 ns. Improving on these conditions will only marginally decrease the error in the lifetime determination. All factors in Eqs. (34)-(36) have very broad shapes (Fig. 6) showing that the system performance degrades slowly when we move away from the optimal conditions. Therefore it is possible to design a system that is quasi-optimal over a relatively broad range of lifetime values. This evaluation resulted in the design of the very simple and cheap system of data acquisition for FLIM presented in this article.
A striking difference between the DFD and the TCSPC approach is the small number of time windows used in the DFD (four windows) with respect to the relative large number of time bins used in the TCSPC when applied to FLIM (64 bins or more). At first sight it would appear that the small number of windows in the DFD approach will limit the capability of measuring very short lifetimes since each window ultimately is several nanoseconds wide. However, the presence of jitter and the sliding window principle makes the measurement of very short lifetime possible and the errors are comparable to the state-of-the-art TCSPC systems.
In our mathematical model of DFD, jitter was presented as a phenomenon which reduces modulation, and therefore reduces statistical accuracy. While this is true, there is also a substantial benefit of having a small amount of jitter in a system. To maintain the ability to measure very short lifetimes, it is necessary for a short lifetime response to be oversampled, so that even very small lifetime changes show up as a response in several bins of the phase histogram as shown sche-matically in Figure 10. This oversampling can be provided by the shape of the excitation pulse and/or by the jitter of the detection electronics. Depending on settings (laser, PMTs, etc.), we measured the jitter on the FV1000 system with the FLIMBox to range from 0.82 to 2 ns. This jitter, according to Eqs. (34)-(36), results in modulation values ranging from 0.95 to 0.84 for the first three harmonics. The FLIMBox produces 256 bins in the phase histogram, but for memory conservation this is usually binned down to 64. This amount of jitter corresponds to an oversampling by around 3-6 bins, providing increased sensitivity to small lifetime values, with minimal expense in uncertainty. If the jitter is too large, resulting in oversampling by a large number of   10. With pulsed excitation, in the absence of jitter, a lifetime value much shorter than the bin size gives ''delays'' all in the same bin so that small changes in lifetime result in identical histograms. The solid curve at bin 1 represents this case in which the center of mass of the distribution cannot be determined with a precision better than the bin size. In the presence of jitter, the broadening of the distribution of delays gives different bin contents even if the changes in lifetime are very small. The dotted curve, which is the convolution of the solid curve with a Gaussian jitter, shows that bin 4 contains less counts than bin 6 so that the center of mass of the distribution can be determined with precision better than the bin size. [Color figure can be viewed in the online issue, which is available at www.interscience. wiley.com.] bins, then the precision of the measurement will reduce according to the modulation value [Eq. (36)] for that jitter. The effect of jitter in the determination of very short lifetime values is more complex than it will appear solely from Eq. (36) because of the effect of limited count statistics in any real measurements. According to Figure 10, in the absence of jitter, there is always a probability that some counts will appear in bin 2 although the majority of the counts will be in bin 1. The amount of counts in bin 2 will also depend on the precise position of the time of photon arrival inside bin 1. Therefore, if we have only a limited amount of counts, statistically we could have either no count in bin 2 or just a few counts. These discrete statistics produce a ''pixelation'' (in the phasor plot) of phasors corresponding to very short lifetime values. The presence of jitter reduces this ''pixelation,'' providing a continuum of phasors centered at the correct phasor position (an equivalent effect is also present in the TCSPC approach). The important point here is that the resolution of very short lifetimes is improved and the error is reduced by the presence of a small amount of jitter and that the ultimate lifetime resolution of the DFD method is only a matter of count statistics. The effect of the count statistics on the resolution of lifetimes can be precisely determined using the mathematical model and therefore a direct quantitative comparison between the DFD and TCSPC method could be done for any given implementation of the two methods.
This insight about jitter reveals new possibilities for FLIM. PMT-based systems usually have a quantum efficiency of around 25% (up to 40% for the GaAs modules), while avalanche photodiode (APD) based systems can have quantum efficiencies of over 60%. These APDs are usually avoided for lifetime measurements because they have a time response jitter on the order of a nanosecond. According to our model, the additional jitter of the APDs will have very little influence on the precision of lifetime measurements. This indicates that DFD-FLIM using APDs may actually yield an overall improvement in lifetime resolution due to the larger number of photons detected.
Our mathematical model shows that an excitation source with a pulse narrower than 908 of the repetition period is near optimal. Further reduction of the excitation pulse width only results in marginal improvement of the precision of the measurement. The DFD method, by virtue of being a frequency domain method, only works with periodic excitation. This is usually the case using pulsed lasers, however the model shows that the narrow pulses from picosecond or femtosecond lasers are not a requirement for FLIM and do not provide substantial benefit over properly modulated diode lasers in the subnanosecond regime.
The DFD approach increases the duty cycle of the detection to 100% as compared with the 50% duty cycle of the analog FD method. Furthermore, the four windows design reaches the plateau in terms of measurement precision for the first harmonic and provides a modulation of 64% for the second harmonic. There is a practical reason why we cannot increase the number of time windows to a much larger number. This is due to the maximal internal frequency of operation of the FPGA chip. For example, to achieve four windows operation at about 50 MHz, we need to generate a fre-quency in the chip which is factor of four larger than the modulation frequency, i.e., 200 MHz. To increase the number of windows by another factor of 2 we will need to generate a frequency of 400 MHz, which is above the maximum limit of the specific FPGA we have used. Other faster chips are commercially available and eventually the precision of the DFD method could be improved. However, according to Figure 6b, the improvement obtained by the eight-window design with respect to a four-window design is marginal. Of course, if the laser repetition frequency is much lower than 50 MHz (for example 20 MHz), we could increase the number of time windows in the FPGA chip accordingly.
In a recent paper (Medine et al., 2007), it is stated that, ''a principle strength of TCSPC-FLIM is statistical accuracy, and this is reflected in the ability to fit two or sometimes three exponential functions to a particular decay.'' We show in the section about the resolution of the decay of enhanced cyan fluorescent protein in two exponential components that DFD-FLIM has comparable statistical accuracy to TCSPC-FLIM for multicomponent analysis, but with substantially reduced hardware expense.
The TAC-based system we used (B&H model 830) has a 125 ns dead time. The dead time limitation becomes particularly relevant in FLIM when the laser samples a very bright pixel. Saturation not only results in lost counts, but can induce significant skew on the lifetime values. Using several cards in parallel removes this limitation but at a very high cost. The FLIMBox hardware collects data at a maximum speed of one count per sampling period. The fundamental limit in count rate is due to the dead time of the discriminator (which is about 10 ns) rather than the processing speed of the FPGA chip. One bottleneck of the particular FPGA we used is the size of the FIFO buffer connected to the USB interface, which limits the average (not the burst) count rate to about 2 MHz.
In conclusion, DFD-FLIM provides a powerful new technique for obtaining lifetime images with a signalto-noise and lifetime measurement ability comparable to the leading techniques in the field, yet with significantly reduced hardware expense. In addition, the theoretical model we developed reveals the system parameter values which must be chosen to optimize the efficiencies of photons for phasor analysis. When DFD-FLIM is combined with the phasor analysis methods, the resulting pair can make lifetime techniques accessible and practical for a wide variety of researchers and applications.