Estimation of multiexponential fluorescence decay parameters using compressive sensing

Abstract. Fluorescence lifetime imaging microscopy (FLIM) is a microscopic imaging technique to present an image of fluorophore lifetimes. It circumvents the problems of typical imaging methods such as intensity attenuation from depth since a lifetime is independent of the excitation intensity or fluorophore concentration. The lifetime is estimated from the time sequence of photon counts observed with signal-dependent noise, which has a Poisson distribution. Conventional methods usually estimate single or biexponential decay parameters. However, a lifetime component has a distribution or width, because the lifetime depends on macromolecular conformation or inhomogeneity. We present a novel algorithm based on a sparse representation which can estimate the distribution of lifetime. We verify the enhanced performance through simulations and experiments.


Introduction
Fluorescence lifetime imaging microscopy (FLIM) is an imaging technique observing the decreasing rate of emitted photon counts with the lapse of time, after applying a laser beam to a cell which is labeled with fluorophores. It provides essential information for cell analysis and can be used for intracellular measurements of fluorescence resonance energy transfer (FRET). 1 The FLIM acquisition methods can be categorized into time-domain and frequency-domain methods. Time-domain methods, such as time correlated single-photon counting (TCSPC) techniques, are more robust to scattering and can measure short lifetimes better than frequency-domain methods. 2,3 FLIM can be obtained using analog mean delay, which is intensity weighted mean delay time; however, it cannot estimate multiple lifetimes; it measures the intensity weighted average lifetime. 4 Recently, many researchers have used FLIM with multiphoton confocal microscopy to noninvasively monitor changes in metabolism. Skala et al. have observed cellular nicotinamide adenine dinucleotide and flavin adenine dinucleotide lifetimes at the earliest stages of cancer development. 5 Roberts et al. have investigated coexisting drugs and their metabolites in the skin using FLIM techniques. 6 Deka et al. have also used FLIM to measure the metabolic rate of cells for in vivo wound healing diagnosis. 7 Accurate estimation of fluorescence lifetime is crucial to the FLIM imaging. A simple estimation method is a nonlinear leastsquares (NLLS) fitting. The NLLS algorithm is a conventional method to estimate lifetimes by minimizing the weighted sum of the squared residuals. [8][9][10][11] However, this method is not effective in low counts, since it does not adopt a Poisson noise probability density function of the photon counting process. FRET has donor and acceptor fluorophores; however, NLLS does not accurately estimate multiexponential decays. Pearson and Neyman's χ 2 -based estimators have been proposed to overcome the Gaussian noise assumption. 12 Sasaki and Masuhara proposed a convolved autoregressive model for fitting multiexponential decay curves. 13 Apanasovich and Novikov also proposed a fitting method, which is mathematically equivalent to Sasaki's and Masuhara method. 14 An improved version was presented by Enderlein and Erdmann, which estimates the optimal fit parameters for multiexponential decay curves by examining an error function. 15 Recently, the maximum-likelihood estimator (MLE) based on the Poisson-likelihood function has proven to be superior in terms of unbiasedness and efficiency from theoretical analysis and computer simulations. 16 Although MLE is an excellent estimator to fit the multiexponential decay parameters, it is sensitive to the initial parameter values. Thus, the choice of the numerical optimization algorithm is significantly important.
The number of fluorescent lifetime components is quite sparse; therefore, a compressed sensing algorithm is suitable for the estimation. We observed that a sparse Poisson intensity reconstruction algorithm (SPIRAL), 17 which is a penalized maximum Poisson-likelihood estimator, is effective for FLIM. SPIRAL has been proposed for the estimation of image intensity from limited photon counts such as tomographic reconstruction in nuclear medicine. Using this algorithm, a unique solution to the minimization can be obtained. This paper proposes a novel algorithm for multiexponential fluorescence decay parameters that can estimate distribution of the lifetime in FLIM. To the best of our knowledge, compressed sensing has not been applied to FLIM yet. We evaluate the performance of our algorithm from simulations and TCSPC data from confocal microscopy, which suffer from very low photon counts. Experimental results demonstrate the efficacy of the proposed approach.

Prior Work
Before proposing a new method for estimating lifetime based on a Poisson distribution with a regularizing function, a brief review of prior work for conventional lifetime estimation is provided. We especially discuss two approaches: modified Neyman's nonlinear least squares (MN-NLLS) and fast fitting, which are embedded in the widely used TCSPC software and SPCImage. 18

Modified Neyman's Nonlinear Least-Squares Method
NLLS algorithm is used in over-determined systems and this algorithm gives the solution that minimizes the sum of squared error, which is the difference between the observed data y ¼ ½y 1 ; y 2 ; : : : ; y N T and the fitting model g ¼ ½gðt 1 Þ; gðt 2 Þ; : : : ; gðt N Þ T . [8][9][10] The objective function is as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 4 9 4 Φ NLLS ¼ The NLLS algorithm is a widely used method to obtain the distribution of fluorescence lifetime, since it is simple and fast. However, NLLS assumes that the measured data are acquired under additive Gaussian noise with equal variance. It is not suitable for photon counts which follow a Poisson distribution. Thus, it is less accurate for fluorescence lifetime estimation, as demonstrated by Kim and Seok. 16 Neyman's χ 2 function has been introduced using the property of the Poisson process; therefore, the variance is proportional to the intensity. However, the original Neyman's χ 2 function is not defined when y i is zero. 16 As a remedy, a modified Neyman's NLLS method is proposed, which is obtained by substituting the variance in the denominator with maxðy i ; 1Þ. This objective function can be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 2 9 0 Φ MN-NLLS ¼ Observation of multiexponential decay is modeled as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 6 3 ; 2 2 9 gðt i Þ ¼ where h i Ã ½expð−t 1 ∕τ 1 Þ · · · expð−t N ∕τ

Fast Fitting Using a Convolved Autoregressive Model
For estimating the decay parameters of photon emission curves convolved with an instrumental response function, a convolved autoregressive model was proposed. [13][14][15] The discrete number of photons observed in a finite number of time channels ½t j − Δt∕2; t j þ Δt∕2 is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 5 6 6 where τ j and f j are the decay constant and the amplitude of the j'th component decay. y k is the number of observed photon counts, and h k is the instrumental impulse response function for the k'th time channel. The z-transform of Eq. (5) yields E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 4 7 1ỹ where z j ¼ expð−Δt∕τ j Þ andhðzÞ is the z transform of the estimated instrumental impulse response function h k . Equation (6) can be rearranged as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 3 8 7ỹ The products can be expanded into sums by defining new coefficients b j and c j : After rearranging and expanding, Eq. (6) can be rewritten in the following form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 2 2 7ỹ The inverse z-transform gives the desired autoregressive convolution: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 1 5 4 It can be written in a compact matrix form.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 7 2 3 y ¼ Dp; with y ¼ ½y Mþ1 ;y Mþ2 ;:::;y N T ; p ¼ ½b 1 ;b 2 ;:::;b M ;c 0 ;c 1 ;:::;c M−1 T ; where N is the number of time sampled data points and M is the number of lifetime components. The general form of Eq. (11) is as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 6 3 ; 5 6 2 y ¼ Dp þ ε; where ε is a random vector of residuals. The optimal values for the vectors p and y can be estimated.
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 5 0 9ŷ ¼ Dp; (13) whereŷ andp are the estimated versions of y and p, respectively. From this notation, we have ε ¼ y −ŷ. The covariance matrix of the vector ε is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 4 4 5 covðεÞ where hεi denotes the expectation of ε. The MLE solution p is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 6 3 ; 3 9 1 p ¼ ðD which enables the calculation of the distribution of multiexponential components f j from Eq. (8).

Problem Formulation
Data collected by counting photons hitting a detector are inherently noisy due to low count levels. Harmany et al. developed SPIRAL for photon limited imaging in the spatial domain such as emission tomography. 17 The algorithm is reported to outperform the state-of-the-art emission tomography algorithms. We found that SPIRAL can be adopted for multiexponential lifetime estimation in the time domain. Under a Poisson process, the observation is modeled as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 2 0 2 y ∼ PoissonðAf Ã Þ; where f Ã ¼ ½f Ã 1 ; f Ã 2 ; : : : ; f Ã M is the amplitude of fluorescence lifetimes to be estimated, the matrix A is given in Eq. (3), and y is a length N vector of the observed photon counts; f Ã and A represented the image intensity and a projection matrix for measurements, respectively, in the original SPIRAL. 17

Poisson Log-Likelihood
The sensing matrix A represents multiexponential decay in FLIM: each column represents an exponential decay of a single lifetime component, while the row direction is the time axis. Thus, the observed photon count is modeled as a linear superposition of sparse lifetime components f Ã . In addition, we need to include the effect of impulse response function of the microscopy. Therefore, the sensing matrix A is the same as Eq. (3).
In the Poisson model, the probability of the observation vector is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 6 6 8 where e i is the i'th canonical basis unit vector. Therefore, the negative Poisson log-likelihood is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 6 0 2 where 1 is an n-vector of ones. Since y i is the constant with respect to f , we can ignore the logðy i !Þ term during optimization. A small parameter β ¼ 10 −10 is added to avoid the singularity problem of a log function when f is zero. The maximumlikelihood estimation for the measurement modeled in Eq. (16) is determined by the minimizing the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 4 4 3

Objective Function
We use an objective function which is a sum of the negative Poisson log-likelihood function and L1-norm regularization term. By adding the L1-norm penalty term to the objective function, we can obtain a sparse solution f Ã . 17,19-23 Exponential decay of the fluorescence can be modeled as a sum of a finite number of lifetime components; therefore, sparse representation is an excellent choice. Function Fðf Þ is the negative Poisson log-likelihood function of the observed photon counts y from the Poisson process. We adopt a penalized negative Poisson log-likelihood objective function with a non-negativity constraint: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 3 2 6 ; 2 4 2 subject to f ≥ 0; (20) where ω is the weighting factor for the regularization term. A sequential quadratic approximation to Fðf Þ is employed for SPIRAL, and the global convergence and uniqueness of the solution are established under a mild set of assumptions. 17 The SPIRAL solves the penalized Poisson log-likelihood function; therefore, it can find multiple exponential decay time constants and amplitudes that are sparse and non-negative.

Numerical Simulations
To validate the proposed algorithm, a simulation study is conducted in order to estimate multifluorescence lifetimes and to compare these results to existing methods. A biexponential decay model for a simulation is as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 6 3 ; 7 3 4 y i ¼ where a and τ denote the photon count and lifetime, respectively. The simulation uses two fluorescence components with lifetimes of 2 and 6 ns. To adjust the Poisson noise variance, we employed a scaling factor α: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 6 3 ; 6 4 4 x ¼ Poissonfαyg: The variance and mean are proportional to α; therefore, signal power is proportional to α 2 , while the signal-to-noise power ratio is proportional to α. As α decreases, the signal-to-noise ratio decreases with α. The observed time interval is from 0 to 19.1 ns and the sampling interval is 0.1 ns. The number of lifetime components for sensing matrix A is 196, which covers from 32 ps to 32 ns in logscale. Since the conventional MLE is sensitive to initial values, we set the initial parameter values to the true values, while the initial values are set to zero for other methods. Conventional methods need to fix the number of lifetime components; therefore, they find the best fit after increasing the lifetime component from one to three in most cases. On the other hand, we obtain the distribution of lifetime components using the proposed algorithm. From the distribution, the lifetime is estimated as the amplitude weighted mean of the distribution in each contiguous group and the corresponding intensity is the sum of the intensity components of the group. Table 1 shows the simulation results of lifetime estimation using the MN-NLLS method, [8][9][10] fast fitting method, 15 MLE method initialized with ground truth, 16 and our algorithm. Mean and standard deviation are as shown in Table 1 as calculated from 500 realizations of the Poisson random variable. The accuracy of algorithms decreases as the scaling Poisson variance factor α decreases. As shown in Table 1, the MN-NLLS and fast fitting methods are quite sensitive to noise: the estimated mean value significantly deviates from true lifetime, especially when α ¼ 0.2. In comparison, the proposed method is robust to noise: the relative error for the 2 ns component is 33%, 37.5%, 7%, and 7.5% for NM-NLLS, fast fitting, MLE, and the proposed method. The proposed method provides similar results the as MLE method initialized with the ground truth, since the proposed algorithm is based on MLE. However, the MLE method is not practical since the initial solution must be the ground truth. Otherwise, it falls in a local minimum, which is quite different from the desired solution.
Typical execution times of various estimation methods in Table 1 are approximately as follows: MN-NLLS method 0.41 s, fast fitting method 0.72 s, maximum-likelihood method 2.10 s, and the proposed method 1.25 s. The compared algorithms were implemented by MATLAB®. The execution times were obtained with an Intel Core i7-3770 processor running at 3.40 GHz. Lifetime depends on the macromolecular conformation and also heterogeneity; therefore, the lifetime component has a distribution or width depending on the conformation. Since the proposed method can estimate the distribution of the lifetime, it results in a more accurate estimation for cell where the number of lifetime components is not known.

Experimental Results
We compare the estimated lifetime of MN-NLLS and fast fitting algorithms to that of the proposed algorithm with FLIM data. A Nikon TE2000-S microscope and Picoquant LDH-P-C-440M laser were used for FLIM imaging. The laser operates at a wavelength of 442 nm and a repetition rate of 20 MHz. A lens for the microscope has a numerical aperture 1.4, 60×. The FLIM image scanner is atomic force microscopy controller by PSIA and the fluorometric detector is the R3809U-07 model by Hamamatsu Photonics. A SPC-830 board by Becker-Hickl is used for data processing. A 473-nm long path cut-off filter (Semrock) was used in front of the detector to block the excitation beam. Photon counts were continuously monitored by a photon counter (EG&G Ortec). The thin film containing dihexyloxacarbocyanine iodide (DOCI) and Coumarin 343 dyes were prepared on a cover glass by the spin-coating method. FITC-Aβ11-25 and FITC-Aβ11-25-DABMI to the HeLa cells were incubated for 20 h.
First, we compare the resolving capabilities with a double exponential decay: we coated a mixture of DOCI and Coumarin 343, whose measured lifetimes are 1.59 and 2.82 ns respectively, on polyacrylic acid polymer. As a measure of quality, we use Neyman's χ 2 , defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 6 3 ; 4 0 3 for y i larger than a threshold. In our experiments, the threshold is set to 5. The numbers of the fit parameters and data points are p and N. The lifetime estimation is accurate when χ 2 is close to one.
To experiment with the biexponential decay, we mixed DOCI and Coumarin 343, which have different lifetimes. TCSPC collected data were analyzed with MN-NLLS, the fast fitting method, and the proposed algorithm, and results are summarized in Table 2. We selected six image positions for brevity: the rest of the image is quite similar to the shown data. The proposed method can more accurately estimate double exponential fluorescence lifetimes than others. The MN-NLLS method estimated double exponential components at two positions among the six samples and the fast fitting method was able to resolve the double exponential fluorescence lifetimes at only one position. The performance of our algorithm is satisfactory in terms of χ 2 , which is close to one. The size of the acquired FLIM image is 248 × 254 pixels and the total time sampling is 256 points with time interval 85.9 ps. The proposed method could resolve the biexponential decay for 81% of them, whereas the MN-NLLS and fast fitting methods detected the biexponential for 37% and 17%, respectively. Even the proposed SPIRAL-  based algorithm could not detect the biexponential decay for the all data points due to low photon counts, but it is still much better than the conventional algorithms. Figure 1 shows the lifetime distribution at the (100, 100) pixel position obtained from the proposed method. We can observe two distinct peaks in the distribution. Each estimated lifetime is the amplitude weighted mean of the components and the intensity is the sum of the constituents. Figure 2 presents the observed photon counts at the (100, 50) pixel position of the model acquired from TCSPC and the exponential decay model from estimated lifetime components over time. The black solid line represents the observed photon counts. The red line with a circle, the blue line with a plus, and the green line with a diamond come from the calculation of the biexponential decay estimated from the proposed method, MN-NLLS, and fast fitting method, respectively. The proposed method accurately estimates double exponents; thus, the fitted exponential decay correctly matches the observed data. While the curves from fast fitting and MN-NLLS closely follow the proposed method for the initial stage for the time interval between 0 and 4 ns, these results deviate from the measured photon counts and the proposed method after roughly 4 ns. The fast fitting results diminish faster than those of the MN-NLLS method. Now we demonstrate the performance of the SPIRAL-based algorithm on FLIM-based FRET, which is a powerful tool that provides information about the vicinity between donor-and acceptor-labeled molecules. 24 FRET occurs when donor and acceptor molecules are placed within a distance that is close (<10 nm) enough to transfer the excited state energy from an initially excited donor to an acceptor. 25 After the energy transfer between a donor and an acceptor, the lifetime of the fluorophore is shortened. We adopted a fluorescence amyloid (Aβ) probe, in which fluoresceine isothiocyanate (FITC) is a donor at the end of the N-terminus of Aβ11-25 and 4-dimethylaminophenylazophenyl-4'-maleimide (DABMI) is an acceptor at the end of the C-terminus. FRET provides us not only the conformational information of Aβ, but also the localization of the peptide in the cellular environment. We observed HeLa cells that were grown in Dulbecco's modified Eagles medium supplemented with 10% (v/v) fetal bovine serum and 1% penicillin. For the confocal and fluorescence lifetime imaging, the cells were plated in a 35-mm glass bottom dish (MatTek Corporation) at low density and grown to 80% confluency on culture dishes at 37°C in a humidified 5% CO 2 atmosphere incubator. At first, we analyze a FLIM image without the DABMI acceptor to observe the data without FRET. The fluorescence of the FITC is known to follow a double exponential decay. 26,27 As shown in Table 3, the proposed algorithm estimates double exponential lifetimes better than MN-NLLS and fast fitting algorithms. Table 3 shows the estimation results at six typical pixel data: MN-NLLS and fast fitting methods cannot resolve biexponential lifetimes for half of the pixels, while the presented method estimates biexponential decay at five positions. The size of the acquired FLIM image is 256 × 256 pixels. The proposed method could resolve the biexponential decay for 97% of them, whereas the MN-NLLS and fast fitting methods detected biexponential for 56% and 68%, respectively. Please note that our method does not estimate just two lifetime components; it estimates the distribution of a lifetime, which is a more accurate model for the fluorescence decay. Figure 3 shows the FLIM images and distributions of lifetimes without and with an acceptor in a HeLa cell estimated by SPC software and the proposed method.  Fig. 3(g) is from the proposed method. To generate these FLIM images with multiple lifetimes, we used the intensity weighted mean lifetime at each pixel and it is coded with pseudo color. The left column shows FLIM images and the right column corresponds to lifetime distributions. As stated above, the pseudo color of the FLIM image in Figs. 3(e) and 3(g), which show shortened lifetimes due to FRET, is more yellowish than Figs. 3(a) and 3(c), respectively. SPC software shows a single lifetime parameter: it shows a peak at 2.5 ns in Fig. 3(b) and the peak position changes to 1.8 ns in Fig. 3(f) after FRET. The same TCSPC data are analyzed by SPIRAL: it estimates two distinctive peaks at 1.8 and 3.4 ns as shown in Fig. 3(d). The amplitude of the 1.8 ns component increases after FRET, while the amplitude of the 3.4 ns component decreases, as observed in Fig. 3(h). The distribution of the lifetime after FRET shows that energy transfer dominantly occurs at 1.8 ns. Since a HeLa cell has biological structure, it cannot be uniformly stained. Moreover, the reaction of FITC differs depending on the local environment such as pH distribution or strain. We believe that the spatial nonuniformity that resulted from the proposed method shows the spatial variation of the structure and the local environment of the cell, while the conventional SPC method cannot resolve the fine details. Thus, the ability to show the fine structure of the HeLa cell reflects an improvement in the quality of the proposed imaging technique. Also please note that the SPIRAL analysis of Figs. 3(c) and 3(g) delineates cell boundary much better than the corresponding images of Figs. 3(a) and 3(e) from the SPC software. Although Fig. 3(c) is the image composed by the data obtained before FRET occurs, the data still have double lifetimes due to the characteristics of FITC and the lifetime in each local region varies significantly depending on pH concentration. 1 Therefore, local homogeneous regions in Fig. 3(c) imply similar pH concentrations in the region. Since Fig. 3(g) is an image after FRET occurs, Fig. 3(g) shows the homogeneous regions with blue or yellow pseudo colors. The bluish region is where FRET does not occur, while the yellowish region represents where FRET occurs, since the region with FRET has shorter lifetimes than the region without FRET. Reference 1 supports this phenomenon. The experimental results demonstrate the applicability of the proposed method to the multiexponential decay of fluorophores from low count cell data, where the maximum count is around 40 photons.

Conclusion
We presented a powerful method to estimate the distribution of fluorescence lifetime using Poisson noise modeling and L1-norm regularization. We use the sensing matrix A based on the multiexponential decay of fluorophores and the instrument response function, together with L1-norm regularization for compressive sensing.
We verified that the presented algorithm shows a more robust estimation of lifetime than the maximum-likelihood estimator based on a Poisson noise from computer simulations and a more accurate distribution with real cell with low photon count images with FRET. SPIRAL has the merit of showing the width or distribution of lifetime components. We need extensive work for the theoretical analysis and practical evaluation of SPIRAL. For future work, we can incorporate a spatial prior in order to improve the quality of FLIM images.