Estimation of Fluorescence Lifetimes Via Rotational Invariance Techniques

Estimation of signal parameters via rotational invariance techniques is a classical algorithm widely used in array signal processing for direction-of-arrival estimation of emitters. Inspired by this method, a new signal model and new fluorescence lifetime estimation via rotational invariance techniques (FLERIT) were developed for multiexponential fluorescence lifetime imaging (FLIM) experiments. The FLERIT only requires a few time bins of a histogram generated by a time-correlated single-photon counting FLIM system, greatly reducing the data throughput from the imager to the signal processing units. As a noniterative method, the FLERIT does not require initial conditions, prior information nor model selection that are usually required by widely used traditional fitting methods, including nonlinear least square methods or maximum-likelihood methods. Moreover, its simplicity means it is suitable for implementations in embedded systems for real-time applications. FLERIT was tested on synthesized and experimental fluorescent cell data showing the potentials to be widely applied in FLIM data analysis.


I. INTRODUCTION
LUORESCENCE lifetime imaging (FLIM) is a powerful tool to study the micro-structure and micro-environments of molecules. FLIM has been widely used throughout modern microscopy, including in the material sciences, biology, chemical analysis and even for clinical diagnosis. Different from traditional fluorescence intensity imaging, which only provides geometric information of tissues or materials, FLIM measures the inherent lifetime of a fluorescent molecule (fluorophore) as it undergoes radiative absorption and subsequent fluorescent relaxation. As the fluorescence lifetime is sensitive to the environment, FLIM can be a good indicator to show how the fluorophore interacts with its microenvironment. Examples include imaging physiological or electrochemical parameters such as Ca 2+ , pH and pO2 [1][2][3][4][5]. When combined with fluorescence resonance energy transfer (FRET) [6][7][8] FLIM is the most robust method to study proteinprotein interactions [1, 2], premalignant lesions [3], molecular metabolism [4] and drug-targeting efficacy [5]. However, despite the potential and significant impact of FLIM, primarily in the biological sciences, estimation of the fluorescence lifetimes remains a significant challenge, particularly with low photon counts systems such as in rapid, live cell imaging. This is becoming increasingly demanding with the development of novel CMOS SPAD array based widefield FLIM systems, which can generate significant volumes of data [9][10][11]. In this paper we present a new, rapid and robust method of extracting lifetime information that requires no prior information on the lifetime components.
There are different FLIM algorithms, mainly in two categories, the time domain (TD) and frequency domain (FD) approaches. FD FLIM mostly uses intensified CCDs synchronized to a modulated excitation source for widefield imaging [12][13][14][15]. The acquisition time is typically a few seconds, but fitting methods are required to extract the lifetimes which can take several seconds to minutes depending on the accuracy requirements. FD lifetime analysis software are usually iterative based. Furthermore accuracy is limited by the CCD array modulation, with the number of phase images (the time bins) typically between 2 to 20. For the TD systems, on the other hand, a pulsed laser is typically used in conjunction with a single photon counting detector, such as a PMT or SPAD Typical TD FLIM instruments either use 1) a time-correlated single-photon counting (TCSPC) module or 2) a time-gated CCD or SPAD [1, 9,16]. For a TCSPC system, the measurements of the time delay between the laser pulses and the detected photon are repeated, and a histogram of time delays is accumulated in which the lifetimes are extracted using fitting algorithms [17]. TCSPC has been the gold standard FLIM technique due to its high timing resolution (typically < 100 ps), and recent developments in multi-channel TCSPC systems [18][19][20]  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TBME.2015.2491364, IEEE Transactions on Biomedical Engineering TBME-00300-2015.R2 2 throughput accordingly demands faster fitting strategies [21]. For a time-gated camera, a series of intensity images at different delays are recorded to extract lifetimes. Similar to TCSPC systems, curve-fitting software are used to calculated lifetimes when the number of gates is larger than 5. The limitations of either FLIM systems are robust and rapid lifetime extractions. This work aims to provide rapid lifetime extractions from the lowest photon signal data possible (highest noise) to enable increasingly rapid imaging solutions. There are two kind of algorithms mainly employed to obtain accurate fluorescence lifetimes. The first type is the widely used fitting methods, including Bayesian [22], maximum likelihood or maximum entropy [23], method of moments [24][25][26] and promptness ratio method [27]. Although fitting methods are precise, they are limiting as 1) they are computationally intensive, 2) they require prior information about how many lifetime components are contained in the data and therefore model selection is required, and 3) they easily converge to local minima. For realistic experiments, particularly in the low photon regime, it may be difficult to precisely know how many lifetime exponents are in every pixel, and researchers usually need to try or choose a proper data model for accurate fitting. The second method for extracting lifetimes are the non-fitting methods, including the phasor algorithm [28,29], Prony's method [30], the integral equation method (IEM) [31], the center-of-mass method (CMM) [32][33][34][35], rapid lifetime determination (RLD) [36][37][38][39]. The phasor method and Prony's method are based on the first order model. The main criticism of all these techniques seems to be that they are all only good for mono-exponentials apart from phasor which may solve biexponentials if you know one component [40]. Gating techniques can be merged with iterative fitting techniques for resolving bi-exponential decays [27]. Kim et al. highlighted the limitations of traditional bi-exponential maximum-likelihood estimation (MLE) fitting and discussing proper gate width and how the IRF affects the estimations. This approach combined gating methods and MLE, and was only demonstrated on data sets of (1, 2) = (1.0ns, 3,9 ns).
In this paper, we propose a new method of lifetime extraction based on a classical algorithm called estimation of signal parameters via rotational invariance techniques (ESPRIT) [22]. ESPRIT has been widely used in array signal processing for direction-of-arrival estimation of emitters, wireless communications, sonar and speech signal processing [41][42][43]. Inspired by ESPRIT, we proposed a new signal model and applied this model for estimation of fluorescence lifetime based on rotational invariance techniques, a system we term FLERIT. FLERIT is 1) non-iterative, 2) capable of resolving multiexponential decays, 3) able to resolve lifetimes when the measurement-window-to-lifetime ratio is less than two and 4) suitable for implementations with embedded hardware for realtime applications. This paper presents the theory (section II), and demonstrates the potential through application to both simulated (section III) and experimental (section IV) data.

II. THEORY
Similar to previously published literature, we suppose the number of the exponential decays is P without considering the instrumental response function (IRF) [44]. Following the model proposed by Hall and Selinger [24], the fluorescence intensity density can be expressed in continuous time domain as where Dj f and j  are the coefficient and the lifetime of the j-th decay component respectively, and   n t is the shot noise.
The i -th bin of the TCSPC histogram is where h is the timing resolution of the TCSPC system.
To reduce the computational complexity and noise, we merged B consecutive bins in the original histogram to create a new histogram with K bins as shown in Fig. 1. The histogram in Fig. 1(a) may be obtained by a TCSPC, whereas the one in Fig. 1 We can arrange the counts from all time bins in the merged histogram as mentioned in [45] as (2) The covariance matrix of (2) [46] can be obtained as where   H  represents the Hermitian transpose and Applying SVD decomposition to (4) As where D is a PP  non-singular matrix. Rewriting the right-handed side of (5) as where 1 C is the last row of A with the dimension 1 P  , 2 C is the first row of A with the dimension 1 P  , and the dimension of the 1 A and 2 Similarly, the left-handed side of (5) can be arranged as then 12 Ψ UU  .
(12) From (11), the eigenvalues of Ψ are equivalent to those of Φ according to the similarity transformation rule [46].
We can calculate Ψ from (12) by applying LU factorization [46] and its eigenvalues , by SVD decomposition. The lifetimes can be estimated accordingly by (9). Finally, we can use the following equation to estimate , 1, , (13) The above outlines the derivation of FLERIT, which can summarize to: Step 1: Reduce the original histogram into K bins; Step 2: calculate the correlation matrix of the new histogram using the covariance matrix (3); Step 3: apply SVD decomposition to X R , (4), to obtain the signal subspace S U ; Step 4: obtain Ψ from (10) and (12); Step 5: calculate the eigenvalues of Ψ by SVD decomposition; Step 6: obtain the lifetimes by (9) and

A. Dynamic Range
We set K = 8 and the photon number in the first bin is 1000. Fig. 2(a) and (b) show the normalized bias ( /   ) and F-value, respectively, in terms of the lifetime. Simulations show that FLERIT has the lowest bias among the four methods. Although CMM shows the lowest F-value with the best photon efficiency, its bias is significant when  > 3ns and a bias correction measure is required to reduce the bias [35]. IEM has the least optimized range in bias. Phasor shows much less efficient when  > 3ns.
The optimized region (F < 4) is from 0.4 to 14ns for FLERIT, 0.4 to 8.8ns for IEM, and 0.4ns to 4.23ns ns for Phasor.

B. Photon Efficiency
The normalized bias and F-value plots in terms of the photon count in the first bin (100~5000) for different methods are shown in Fig. 3(a) and (b). Here we use K = 8 and  = 3ns.
Again, FLERIT shows the lowest bias. It is interesting that the F-value of FLERIT is similar to IEM, larger than CMM and less than Phasor. But CMM has the worst bias performance, unless a bias correction is carried out. Figure 3(a) shows that both the bias and the normalized F-value should be independent of the photon count as expected.

C. Performances in terms of K
The normalized bias and F-value plots in terms of K (4 < K < 32) for different methods are shown in Fig. 4(a) and (b). Here we set  = 3ns, the photon number in the first bin is 1000, and the number of the time bins is 1024. For FLIRET, the F-value degrades as K increases. FLERIT has similar bias performances with IEM, whereas Phasor and CMM are significantly biased. Figure 4(a) shows that Phasor and CMM favor a larger K. This means FLERIT can be used to resolve histograms obtained by gated FLIM systems as well as TCSPC systems.

IV. LIFETIME RESOLVABILITY ANALYSIS
To test whether FLERIT can resolve bi-exponential decays robustly, Monte Carlo simulations were carried out with 1  = 1.5ns, 1.5ns < 2  < 6ns,     Figure 5 also shows that the threshold of c N is about 10000.

V. ANALYSIS ON SYNTHESIZED MULTI-DECAY FLIM DATA
In TCSPC FLIM experiments the photons collected at each image pixel is limited, either due to the time taken to obtain a viable histogram, limiting for live cell imaging, or due to photobleaching. Pixel binning is therefore often applied to improve the signal-to-noise (SNR) of FLIM images at the sacrifice of spatial resolution. To illustrate the FLERIT methodology and advantages, we first use synthesized data with the number of photons limited both A) without and B) with pixel binning.

A. Without pixel binning
When there is only one lifetime, the number of eigenvectors of the signal space S U is one. However, an interesting feature found in FLERIT is that when there are multiple lifetimes, the eigenvector corresponding to the biggest eigenvalue of X R is actually a linear combination of all lifetimes when the photon count is limited, i.e. The laser repetition rate is set to be 80MHz ( T = 12.5ns).
In the first simulated case we define the primary and secondary lifetimes 1  = 2ns and 2  = 5ns for all areas and and C, respectively. The image size is 256×256 pixels. The photon count at the first bin is 500, and the number of time bins in the histogram is 256 (the histogram was merged into a new one with K = 8). The photon count (intensity) of the original data is shown as Fig. 6(a).The averaged lifetime image obtained by FLERIT is shown as Fig. 6(b). The mean (AVE), standard deviation and F-value (normalized precision; F = 1 for the ideal case, and F > 1 or F >> 1 for realistic FLIM algorithms) of the calculated lifetimes for the three areas are listed in Table I. The precision of the average lifetime is similar to the single-exponential case in Section III. Consider the second case of a tri-exponential decay where 1  = 2ns; 2  = 3ns; 3  = 5ns; fD1 = 0.4, 0.33, and 0.2, fD2 = 0.4, 0.33, and 0.2, fD3 = 0.2, 0.34, and 0.6, for the areas A, B, and C, respectively. The photon count (intensity) of the original data is shown as Fig. 7(a).The averaged lifetime image obtained by FLERIT is shown as Fig. 7(b).  Tables I and II confirm that FLERIT offers an interesting feature similar to the previously reported IEM [33] when it deals with multi-exponentials: This is a useful feature, as in some applications such as FRET-FLIM experiments [7], (11) can be used to estimate the FRET efficiency. In many biological applications however, it is desirable to estimate j  and , 1, , Dj f j P  and as previously discussed IEM is only a single exponential approximation.

B. With pixel binning:
Due to limited photons in the histogram, it is challenging to estimate 1 D f accurately using (13). Typically pixel binning is used to increase the photon count by trading off the spatial resolution. Using the synthesized data presented in Fig. 6, we adopted a summation based binning procedure as shown in Fig.  8. The intensity after binning is shown in Fig. 8(a). After binning, 1  , 2  , 1 D f and 2 D f can be estimated, and the averaged lifetime can be calculated as shown in Fig. 9(b) using (14). The performances including the mean, standard deviation and F-value of calculated parameters are listed in Table III. The  table contains more parameters than Table II, as FLERIT resolves all lifetime components and proportional coefficients. The F-value is slightly worse than that in Table I,   TBME-00300-2015.R2 6 Fig. 11(a) shows the 1 D f image after binning. Fig.11 Fig. 6(a). The SNR of AVE  is much better than each individual j  and Dj f . The ratio of the measurement window to 2  (12.5ns/5ns), is only 2.5, much smaller than the recently proposed bi-exponential algorithms [48], indicating that the duty cycle of the laser repetition can be higher giving better photon efficiency.

VI. EXPERIMENTS
To test FLERIT on real data, FLIM experiments were carried out on HeLa cells ubiquitously expressing EGFP using a commercial scanning confocal FLIM system

A. Experimental set up
Data was acquired using a Leica SP5 scanning confocal microscope fitted with a PicoHarp 300 TCSPC module. Excitation was with a tunable white light laser (WLL) operating at 488 nm and 40 MHz. Detection was with a single channel MPD SPAD, collecting the majority of EGFP emission. All images where Nyquist sampled, 512x512 pixels and with predefined total image integration times set to 10s, 60s, 180s and 600s.

B. Sample preparation
HeLa cells were plated onto 25mm glass coverslips previously coated with 50ug/ml poly-D-Lysine hydrobromide (UV irradiated for sterility) and grown for 24 hours at 37 o C and 5% CO2. Cells were cultured in DMEM (Gibco, 31053) supplemented with 100U/ml penicillin, 100ug/ml streptomycin (Gibco, 15140), 10% heat-inactivated foetal bovine serum (Gibco, 10500064), 1X Glutamax (Gibco, 35050) and 1mM sodium pyruvate (Gibco, 11360). Following 24 hour growth, cells were transfected using Turbofect transfection reagent (Thermo Scientific, R0531) with 2ug pEGFP-N1, a discontinued Clontech plasmid encoding enhanced green fluorescent protein, and incubated for a further 24 hours at 37 o C and 5% CO2 to allow expression of the encoded EGFP. Cells were washed with 1X phosphate buffered saline, fixed with 4% paraformaldehyde and blocked with 50mM ammonium chloride prior to mounting on slides with MOWIOL 4-88 and allowing to set overnight before imaging. Fig. 12 shows the epifluorescence image an example cell with ubiquitously expressed EGFP. Fig. 13(a)-(d) show the average fluorescence lifetime image for the acquisition time of 10s, 60s, 180s, and 600s (the maximum photon count are 54, 251, 756, and 1939 respectively). The figures show that the deviations of the lifetime decrease as the acquisition time is increased. Fig. 14(a)-(d) show the histograms of AVE  for the acquisition time of 10s, 60s, 180s, and 600s. Figs. 14 show that the standard deviation can be improved with a longer acquisition and it is inversely proportional to the square root of the acquisition, in agreement with the conclusion given in [28]. Fig. 15(a)-(d) show the histograms of 1  and 2  for different acquisitions. The bi-exponential ingredient is 1.58%, 2.12%, This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TBME.2015.2491364, IEEE Transactions on Biomedical Engineering TBME-00300-2015.R2 7 2.12% and 2.21% respectively for 10s, 60s, 180s, and 600s. The peaks of 1  and 2  histograms are located at around 850ps and 3ns. The average lifetime is about 2.8ns, in accordance with Ref [49].

D. Time consumption of the data analysis
We have run the data analysis using MATLAB® on DELL Optiplex 7010 desktop. For a 512x512 image, it takes 26.7s. If FLERIT is implemented in a hardware similar to Ref. [31], the computational burden can be further decreased by adopting fast multistage Wiener filtering method [50], Lanczos algorithm [51], and propagator method [52]. With more and more hardware multipliers or intellectual property cores embedded in DSP processors and FPGA circuits, the proposed method can be realized in such embedded systems for real-time applications.

VII. CONCLUSIONS
In this paper, we proposed a new method called FLERIT. The derivations of FLERIT have been carried out by introducing a signal model. The performances of FLERIT were demonstrated on both synthesized and experimental data. The new method does not require any prior information and it can be applied to both gated CCDs and TCSPC systems with limited number of timing channels. For more accurate analysis, the model can be extended to include the IRF, but it is an independent work not covered in this paper. Simulations and experiments show that FLERIT can provide single-exponential average fluorescence lifetimes similar to the previously reported IEM method or multiple exponential analysis. FLERIT can extract fluorescence lifetimes by only a few time gates, which can reduce the data throughput between a parallel TCPSC front-end and a data analyzing system. The computation burden of FLERIT is much less than traditional fitting methods making it suitable for implementations in embedded systems.