Deconvolution of fluorescence lifetime imaging microscopy by a library of exponentials

Fluorescence lifetime microscopy imaging (FLIM) is an optic technique that allows a quantitative characterization of the fluorescent components of a sample. However, for an accurate interpretation of FLIM, an initial processing step is required to deconvolve the instrument response of the system from the measured fluorescence decays. In this paper, we present a novel strategy for the deconvolution of FLIM data based on a library of exponentials. Our approach searches for the scaling coefficients of the library by non-negative least squares approximations plus Thikonov/l2 or l1 regularization terms. The parameters of the library are given by the lower and upper bounds in the characteristic lifetimes of the exponential functions and the size of the library, where we observe that this last variable is not a limiting factor in the resulting fitting accuracy. We compare our proposal to nonlinear least squares and global non-linear least squares estimations with a multi-exponential model, and also to constrained Laguerre-base expansions, where we visualize an advantage of our proposal based on Thikonov/l2 regularization in terms of estimation accuracy, computational time, and tuning strategy. Our validation strategy considers synthetic datasets subject to both shot and Gaussian noise and samples with different lifetime maps, and experimental FLIM data of ex-vivo atherosclerotic plaques and human breast cancer cells. © 2015 Optical Society of America OCIS codes: (170.1610) Clinical applications; (170.3650) Lifetime-based sensing; (170.3880) Medical and biological imaging; (170.6510) Spectroscopy, tissue diagnostics; (170.6920) Time-resolved imaging. References and links 1. L. Marcu, P. M. W. French, and D. S. Elson, Fluorescence Lifetime Spectroscopy and Imaging: Principles and Applications in Biomedical Diagnostics (CRC, 2014). 2. N. Anthony, K. Berland, and P. Guo, “Principles of fluorescence for quantitative fluorescence microscopy,” in FLIM Microscopy in Biology and Medicine, A. Periasamy and R. M. Clegg, eds. (Chapman and Hall/CRC, 2009), pp. 35–63. 3. A. Periasamy and R. M. Clegg, FLIM Microscopy in Biology and Medicine (Chapman and Hall/CRC, 2009). 4. S. Shrestha, B. Applegate, J. Park, X. Xiao, P. Pande, and J. Jo, “High-speed multispectral fluorescence lifetime imaging implementation for in vivo applications,” Opt. Lett. 35(15), 2558–2560 (2010). 5. M. O’Donnell, E. R. McVeigh, H. W. Strauss, A. Tanaka, B. E. Bouma, G. J. Tearney, M. A. Guttman, and E. V. Garcia, “Multimodality cardiovascular molecular imaging technology,” J. Nucl. Med. 51(1), 38S–50S (2010). #243670 Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 (C) 2015 OSA 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23748 6. J. Park, P. Pande, S. Shrestha, F. Clubb, B. E. Applegate, and J. A. Jo, “Biochemical characterization of atherosclerotic plaques by endogenous multispectral fluorescence lifetime imaging microscopy,” Atherosclerosis 220(2), 394–401 (2012). 7. V. De Giorgi, D. Massi, S. Sestini, R. Cicchi, F.S. Pavone and T. Lotti, “Combined non-linear laser imaging (twophoton excitation fluorescence microscopy, fluorescence lifetime imaging microscopy, multispectral multiphoton microscopy) in cutaneous tumours: first experiences,” J. Eur. Acad. Dermatol. Venereol. 23(3), 314–316 (2009). 8. N. Galletly, J. McGinty, C. Dunsby, F. Teixeira, J. Requejo-Isidro, I. Munro, D. Elson, M. Neil, A. Chu, P. French, and G. Stamp, “Fluorescence lifetime imaging distinguishes basal cell carcinoma from surrounding uninvolved skin,” Br. J. Dermatol. 159(1), 152–161 (2008). 9. J. M. Jabbour, S. Cheng, B. H. Malik, R. Cuenca, J. A. Jo, J. Wright, Y.-S. L. Cheng, and K. C. Maitland, “Fluorescence lifetime imaging and reflectance confocal microscopy for multiscale imaging of oral precancer,” J. Biomed. Opt. 18(4), 046012 (2013). 10. M. Mycek, K. Schomacker, and N. Nishioka, “Colonic polyp differentiation using time-resolved autofluorescence spectroscopy,” Gastrointest. Endosc. 48(4), 390–394 (1998). 11. A. J. Walsh, R. S. Cook, M. E. Sanders, L. Aurisicchio, G. Ciliberto, C. L. Arteaga, and M. C. Skala, “Quantitative optical imaging of primary tumor organoid metabolism predicts drug response in breast cancer,” Cancer Res. 74(8), 5184–5194 (2014). 12. A. J. Walsh, R. S. Cook, H. C. Manning, D. J. Hicks, A. Lafontant, C. L. Arteaga, and M. C. Skala, “Optical metabolic imaging identifies breast cancer glycolytic levels, sub-types, and early treatment response,” Cancer Res. 73(20), 6164–6174 (2013). 13. B. K. C. Lee, J. Siegel, S. E. D. Webb, L. S. Fort, M. J. Cole, R. Jones, K. Dowling, M. J. Lever, and P. M. W. French, “Application of the stretched exponential function to fluorescence lifetime imaging,” J. Biophys. 81(3), 1265–1274 (2001). 14. S. Pelet, M. J. R. Previte, L. H. Laiho, and P. T. C. So, “A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation,” Biophys. J. 87(4), 2807–2817 (2004). 15. S. C. Warren, A. Margineanu, D. Alibhai, D. J. Kelly, C. Talbot, Y. Alexandrov, I. Munro, M. Katan, C. Dunsby, and P. M. W. French, “Rapid global fitting of large fluorescence lifetime imaging microscopy datasets,” PLOS One 8(8), e70687 (2013). 16. J. A. Jo, Q. Fang, T. Papaioannou, and L. Marcu, “Fast model-free deconvolution of fluorescence decay for analysis of biological systems,” J. Biomed. Opt. 9(4), 743–752 (2004). 17. A. S. Dabir, C. A. Trivedi, Y. Ryu, P. Pande and J. A. Jo, “Fully automated deconvolution method for on-line analysis of time-resolved fluorescence spectroscopy data based on an iterative Laguerre expansion technique,” J. Biomed. Opt. 14(2), 024030 (2009). 18. P. Pande, and J. A. Jo,“Automated analysis of fluorescence lifetime imaging microscopy (FLIM) data based on the Laguerre deconvolution method,” IEEE Trans. Biomed. Eng. 58(1), 172–181 (2011). 19. J. Liu, Y. Sun, J. Qi, and L. Marcu, “A novel method for fast and robust estimation of fluorescence decay dynamics using constrained least-squares deconvolution with Laguerre expansion,” Phys. Med. Biol. 57(4), 843–865 (2012). 20. O. Gutierrez-Navarro, D. U. Campos-Delgado, E. R. Arce-Santana, K. C. Maitland, S. Cheng, J. Jabbour, B. Malik, R. Cuenca, and J. A. Jo, “Estimation of the number of fluorescent end-members for quantitative analysis of multispectral FLIM data,” Opt. Express 22(10), 12255–12272 (2014). 21. O. Gutierrez-Navarro, D. Campos-Delgado, E. Arce-Santana, M. Mendez, and J. Jo, “Blind end-member and abundance extraction for multi-spectral fluorescence lifetime imaging microscopy data,” IEEE J. Biomed. Health Inform. 18(2), 606–617 (2014). 22. H. Xu and B. W. Rice “In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique,” J. Biomed. Opt. 14(6), 064011 (2009). 23. F. Fereidouni, G. A. Blab, and H. C. Gerritsen, “Blind unmixing of spectrally resolved lifetime images,” J. Biomed. Opt. 18(8), 086006 (2013). 24. F. Fereidouni, G. A. Blab, and H. C. Gerritsen, “Phasor based analysis of FRET images recorded using spectrally resolved lifetime imaging,” Methods Appl. Fluoresc. 2(3), 035001 (2014). 25. D. U. Campos-Delgado, O. Gutierrez Navarro, E. R. ArceSantana, and J. A. Jo, “Extended output phasor representation of multi-spectral fluorescence lifetime imaging microscopy,” Biomed. Opt. Express 6(6), 2088–2105 (2015). 26. J. Zhang and B. Morini, “Solving regularized linear least-squares problems by the alternating direction method with application to image restoration,” Electron. Trans. Numer. Anal. 40, 356–372 (2013). 27. P. Gong and C. Zhang, “A fast dual projected newton method for l1-regularized least squares,” in Proceedings of the Twenty-Second international joint conference on Artificial Intelligence 2, 1275–1280 (2011). 28. G. Landi, “A modified newton projection method for l1-regularized least squares image deblurring,” J. Math. Imaging Vis. 51(1), 195–208 (2015). 29. M. V. Afonso, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “An augmented Lagrangian approach to the constrained optimization formulation of imaging inverse problems,” IEEE Trans. Image Process. 20(3), 681–695 (2011). #243670 Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 (C) 2015 OSA 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23749 30. K. Vishwanath, and M. A. Mycek, “Do fluorescence decays remitted from tissues accurately reflect intrinsic fluorophore lifetimes?” Opt. Lett. 29(13), 1512–1514 (2004). 31. J. G. Proakis and D. G. Manolakis, Digital Signal Processing, 4th ed. (Prentice Hall, 2006). 32. J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed. (Springer-Verlag, 2006). 33. R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge University, 1985). 34. The MathWorks Inc., Parallel computing toolbox user’s guide R2015a, http://www.mathworks.com/help/ pdf doc/ distcomp/distcomp.pdf. 35. The MathWorks Inc., Optimization toolbox user’s guide R2015a, http://www.mathworks.com/help/pdf doc /optim /optim tb.pdf.


Introduction
Fluorescence is an important resource for quantitative characterization of samples [1][2][3].In fluorescence lifetime imaging microscopy (FLIM), the optic-electronic instrumentation captures the time-resolved fluorescence emission to a laser excitation at each pixel of the imaged sample [4].This technique allows early and non-invasive diagnosis for different pathologies, as cardiovascular and dermatology diseases [5][6][7][8], oral pre-cancer conditions [9], colonic dysplasia [10], or measure therapeutic responses of anticancer drugs [11,12].However, a precise evaluation of the FLIM information requires an initial processing stage in order to deconvolve the instrument response of the system, and to extract the average lifetime map of the sample to identify the fluorophores [13][14][15][16][17][18][19].Other strategies for the quantification of the FLIM information are focused on the direct analysis of the fluorescence decays by a linear unmixing approach [20][21][22], or analyze the datasets in a lower-dimensional domain by the phasor perspective [23][24][25].
The deconvolution process is an inverse problem that in general is ill-posed.The standard deconvolution method assumes a multi-exponential model for the fluorescence decay [13][14][15].In this approach, the parameters of the model (characteristic lifetimes and scaling coefficients) are identified jointly at each spatial location of the sample by nonlinear least squares approximations.When the characteristic lifetimes are considered common to all spatial points, a global perspective is adopted, which largely reduced the processing time but could be affected by local minima during the optimization [13][14][15].Alternatively, a method based on a Laguerre expansion considers that the fluorescence impulse response can be expressed as linear combinations of the Laguerre base [16][17][18], where by a subsequent least-squares optimization, the scaling coefficients of the base are computed to approximate the measured fluorescent decay.However, this methodology could result in a fluorescence response with a non-monotonic decreasing pattern.So, in this case, a constrained least squares approximation has to be pursued to enforce the monotonic decreasing time profile [19].Moreover, in [18], the shape parameter in the Laguerre basis is adaptively optimized per pixel in the FLIM dataset to improve the fitting accuracy, although the computational complexity is raised.In fact, the order of the Laguerre basis has to be also carefully selected in an ad hoc manner to achieve the best estimation according to the studied FLIM dataset.
In this context, this paper introduces a novel methodology for the deconvolution problem that relies on a library of exponentials, where the characteristic lifetimes in the library are chosen from pre-defined lower and upper bounds for the studied sample.As our analysis will show, this range of variation can be rather loose, and the number of elements in the library does not need to be large to achieve good accuracy.The scaling coefficients of the library elements are chosen by non-negative least squares approximations (NNLSA) plus a Thikonov/l 2 or l 1 regularization term [26][27][28][29], where this process is parallelized to speed up convergence.We present an advantage of our approach based on Thikonov/l 2 regularization in terms of estimation accuracy, computational time and tuning strategy of its parameters with respect to nonlinear least squares and global non-linear least squares approaches [13][14][15], and a constrained Laguerre-base method [19].We present our validation strategy under synthetic datasets subject to both shot and Gaussian noise, and experimental FLIM data of ex-vivo atherosclerotic plaques [6] and human breast cancer cells [11].
The notation used in this paper is described next.Scalars are denoted by italic letters, vectors and matrices by lower and upper-case bold letters, and sets by italic calligraphic upper-case letters.R represents the real numbers, R n n-dimensional real vectors, R n×m real matrices of dimensions n × m, and card(X ) the cardinality of a set X .For a real vector x, the transpose operation is denoted by x , the Euclidean norm by x = √ x x, and x 0 represents that each component in the vector is positive or zero.An n-dimensional vector filled with ones (zeros) is represented by 1 n (0 n ), and I n denotes the identity matrix of order n.For a vector x ∈ R n , diag(x) ∈ R n×n characterizes a square matrix with diagonal terms given by x.For a random variable x, x ∼ N (0, σ 2 ) represents that x is normally distributed with zero mean and variance σ 2 .The list of the acronyms used throughout the paper is presented in Table 1.

Problem formulation
We assume that the FLIM instrumentation provides discrete-time measurements of the fluorescence decays of a sample over a spatial domain of K points in the dataset, and with a sampling period of T seconds [1][2][3].We also assume that there is available a measurement of the instrument response of the system.By considering a time window of L samples and a causal response [31], the observation model for the l-th time sample at k-th spatial point is described by where y k [l] and h k [l] denote the measured fluorescence decay and impulse response samples that depend both on the spatial location k, respectively, u[l] characterizes the instrument response samples, represents the convolution operator, and v k [l] describes the random noise related to the instrumentation or measurement uncertainty.In fact, the effect of scattering can be modeled as a scaled factor A k > 0 of the instrument response u[l] in the observation model of Eq. ( 1) [1][2][3], i.e. under the effect of scattering we have now where δ [l] denotes the discrete-time delta function.Although this component could potentially introduce a distortion of the estimated impulse response function, that is this a common problem to all deconvolution methods, and has been shown previously that the scattering effect is not significant for endogenous tissue FLIM [30].The observation model in Eq. ( 1) can be written in a vector notation as where the resulting input matrix U has a toeplitz structure [33] that does not depend on the spatial location of the sample.A diagram of the observation model in our studied framework is presented in Fig. 1.We define the set of FLIM measurements as Y = {y 1 ,...,y K }.Hence the deconvolution problem estimates for the impulse response vector ĥk at each spatial point, such that the estimated fluorescence decay vector ŷk = U ĥk approximates the measurement one y k [1][2][3].As any inverse problem, the deconvolution is an ill-posed formulation that can have multiple solutions.In our problem, there are two key assumptions that help to formulate and bound the search space for the feasible solutions: A. The instrument response u[l] is common to all K spatial points in the dataset, and its samples are non-negative and normalized to sum one in time domain, i.e.
B. The fluorescence impulse response h k [l] at each spatial sample k and time instant l can be represented by the conical combination of a library of N pre-defined exponential functions {e −l/τ 1 , e −l/τ 2 ,...,e −l/τ N } plus a constant offset term, see Fig. 2(a).In this way, the library of exponential functions has known and bounded characteristic lifetimes {τ 1 , τ 2 ,..., τ N }, i.e. τ min ≤ τ n ≤ τ max and 0 < τ min < τ max for all index n ∈ [1, N].As a result, the fluorescence impulse response h k [l] is modeled as: where the non-negative positive coefficient c k,n ≥ 0 represents the contribution of the n-th exponential mode into the fluorescence response at k-th spatial point, and c k,N+1 the offset component.Thus, by the observation model in Eq. ( 3) and putting aside the random or uncertain term, any measured fluorescence decay vector y k in the dataset is assumed to belong to the convex cone where In this way, the L-dimensional vectors {Ub n } N n=1 will represent a library of fluorescence decays that will be used to approximate each measurement in the FLIM dataset, see Fig. 2

(b).
Assumption A is important to avoid numerical scaling problems in the estimation of the coefficients {c k,n } N+1 n=1 for each location k.Meanwhile, assumption B allows to formulate the estimation problem as a non-negative least-squares approximation, which can be solve efficiently by numerical optimization algorithms [32].Furthermore, for any spatial location k, we want to have just a small subset of the terms in the library of exponentials with significant weights {c k,n } N n=1 , so we will like to induce some regularization into the solution.As a result, we propose to use a Thikonov/l 2 or l 1 regularization method for this purpose [26][27][28][29].We see  two important properties of our approach based on a library of exponential functions in regard to previous schemes from the literature: I.With respect to the standard deconvolution method based on nonlinear least squares [13][14][15], our strategy is linear in the unknown parameters {c k,n } N+1 n=1 , and the user does not need to choose a pre-defined order for the approximation model.II.With respect to the Laguerre base method [16][17][18][19], our methodology guarantees a smooth fluorescent impulse response with a monotonic decreasing pattern, and also, the order of the approximation does not have to be selected a priori and there are no tuning parameters per spatial point.
Thus, in our approach, just the minimum and maximum time constants (τ min , τ max ) of the library of exponentials in Eq. ( 5) have to be defined, which is a more easy to tune condition based on the expected lifetimes in the studied sample.Furthermore, the characteristic lifetimes {τ n } N n=1 in Eq. ( 5) could be easily generated based on regularly spaced grid in the interval [τ min , τ max ].Now, by considering the library of exponential functions, the observation model in Eq. ( 1) at k-th spatial sample can be expressed compactly as where i.e. matrix B contains in its columns the information of the library of exponentials, and vector c k its unknown scaling coefficients.Note that a direct least squares approximation from Eq. ( 7) to compute c k could be really fast but not feasible, since some of the elements in c k can be negative, which will not represent a meaningful solution for the deconvolution problem.

Optimal approximations
First, the input matrix U can be computed from Eq. ( 2) by the samples {u[l]} L−1 l=0 .With this departing information, the scaling coefficients vector c k are estimated for each spatial point k ∈ [0, K − 1] by a non-negative least-squares approximation (NNLSA) [32] with a Thikonov/l 2 regularization term as follows min where α > 0 is the regularization weight.An important property of Eq. ( 10) is that the coefficients vector c k for each location k can be solved in parallel over the whole dataset, since there is no inter dependency or homogeneity spatial condition in our deconvolution formulation.Meanwhile, another alternative is to consider an l 1 regularized least squares approximation for the deconvolution process [26][27][28]: Hence in both optimization problems in Eqs. ( 10) and ( 11), there are two important parameters to set: the number of elements in the library of exponentials N, and the regularization weight α.Next, our evaluation will show that N does not have a profound impact in the estimation accuracy, but affects greatly the computational time.In the meantime, α can affect the accuracy if this parameter is large, as well as the processing time.In our evaluation, we will denote as DELE-L2 the solution of the deconvolution estimation by the library of exponentials with the Thikonov regularization in Eq. ( 10), and DELE-L1 with the l 1 regularization term in Eq. ( 11).

Synthetic and experimental validation
In this section, we present the validation of the DELE-L2 and DELE-L1 algorithms by considering two cases: synthetic and experimental FLIM datasets, where all the processing was carried out in Matlab.For both cases, the estimation errors on the measured fluorescence decays will be evaluated, and in the synthetic experiments, we also quantify the error on the resulting fluorescence impulse responses and average lifetimes [1][2][3].Furthermore, our synthetic validation scheme takes into account different scenarios for Gaussian and shot noise in the fluorescence decays [20].In the synthetic and experimental evaluations, we compare our approach against three well-known deconvolution approaches from the literature: nonlinear least squares (DENLS) and global nonlinear least squares (DEGNLS) with a multi-exponential model plus an offset component [1], [15], and a Laguerre-base (DELB) method with constrained third order time derivative [19].
The average computational time per FLIM dataset was computed in the Matlab implementations for the overall lifetime estimations by using a Pentium Intel Core i7-4770 3.5 GHz quad-core processor and 32 GB RAM computer.For this purpose, we employed the commands tic and toc from Matlab in order to measure accurately the execution time.Furthermore, during the evaluation, Matlab was the only active process in the computer.In addition, we execute "for loops" in parallel during the numerical implementations by the Parallel Computing Toolbox of Matlab [34], and we use the command lsqnonneg to solve the NNLSA in Eq. ( 10) from the Optimization Toolbox [35].Meanwhile, we implemented directly the solution described in [28] for the l 1 regularized least squares approximation in Eq. (11).In DENLS and DEGNLS, the characteristic lifetimes were estimated based on constrained minimization of the estimation error of the fluorescent decays, and the scaling coefficients on standard least squares.For this purpose, we employ the function fmincon from the Optimization Toolbox in Matlab [35] by including the gradient information of the error cost function, and adopting the parallel configuration of the algorithm provided by the toolbox.Note that DEGNLS is a fast strategy, since the characteristic lifetimes of the multi-exponential model are assumed to be common over the whole dataset [1][2][3]15].Finally, we implement the methodology in [19] of DELB that is based on a dual formulation of a constrained quadratic optimization problem.This method relies also on NNLSA, which is solved by lsqnonneg at each spatial point [35].

Synthetic evaluation
Three synthetic datasets were generated by using the measured instrument response in [6] with a sampling interval T = 250 ps, a length of 186 samples (L = 186), and a multi-exponential model for the impulse responses with 2, 3 and 4 fluorophores.Each dataset has a spatial dimension of 100 × 100 pixels, i.e. there are K = 10, 000 spatial points per dataset.The instrument response signal {u[l]} L−1 l=0 has a sharp rising time and a subsequent exponential decay with a full-width at half-maximum (FWHM) of 1.525 ns, and its normalized to sum one to avoid numerical problems, i.e. ∑ l u[l] = 1.The fluorescence impulse response model in the synthetic datasets is expressed as a sum of 2, 3 or 4 exponential terms at each spatial point k [1-3]: where the abundance at spatial location k of the i-th characteristic mode τ i is expressed as δ k,i ∈ [0, 1] and it was chosen according to the maps in Fig. 3, and δ o = 0.001 defines an offset term.We denote the 2nd order dataset as (A), the 3rd order as (B), and the 4th order as (C) in Fig. 3.In all the datasets, the first component δ k,1 in dominant in the center of the sample (see first column in Fig. 3.Meanwhile, in dataset (A), the second component δ k,2 is more abundant over the boundary of the sample; in dataset (B), the second component δ k,2 is dominant in the right upper corner with some presence in the lower part of the sample, and the third one δ k,3 is more abundant in the left side but it has also dominance in the lower right corner; and finally in dataset (C), the highest abundances of the second δ k,2 , third δ k,3 and fourth δ k,4 components are alternated in the upper right, upper left and lower right corners, respectively.Hence, the three datasets (A), (B) and (C) describe rich heterogeneous samples for 2nd, 3rd and 4th order models with spatially dependent average lifetime maps.The characteristic lifetimes τ i in Eq. ( 12) are defined in Table 1.Next, the synthetic uncorrupted fluorescence decays y o k [l] are computed by applying the convolution operator in Eq. ( 1), i.e.
In the electronic/optic instrumentation used in FLIM, there is measurement uncertainty due to low photon counts which is characterized by shot or Poisson noise [1][2][3], plus the effect of electrical noise which is modelled by Gaussian noise.Hence, to reproduce a practical scenario in our synthetic tests, we included both shot and Gaussian noise to the synthetic datasets according to the following model [20]: where k ) represent normal random variables related to the shot and Gaussian noise terms, respectively, where the variances σ 2 k and ς 2 k are selected with respect to desired peak signal-to-noise (PSNR) and signal-to-noise ratios (SNR): In this way, for our model in Eq. ( 13), the standard deviation of the shot noise component is time-varying and equal to y o k [l] • σ k , which is adjusted by σ k to achieve the desired PSNR at k-th pixel.By recalling that h k and y k denote the true vectors for the impulse response and fluorescent decays, respectively, and ĥk and ŷk their estimations by the deconvolution algorithms, we compute four performance indexes to evaluate the goodness of fit and accuracy of the studied methodologies over the whole dataset [1][2][3]20]: where h = (1/K) ∑ K−1 k=0 h k and ȳ = (1/K) ∑ K−1 k=0 y k are the mean impulse response and mean fluorescent decay in the dataset, respectively.Hence the indexes GoF H and GoF Y evaluate the fitting accuracy with respect to the variability in the dataset of the impulse responses and fluorescence decays, respectively; meanwhile E H and E Y quantify the estimation error energy with respect to the dataset energy of the same signals.In addition, we propose a new performance metric in regard to the measured τ k and estimated τk average lifetimes of the impulse responses [20] at each spatial point k: where t = [0 T s 2T s . . .(L − 1)T s ] represents a vector of sampling instants.Hence we define the goodness of fit with respect to the percentage average lifetime estimation error as  For DELE-L2 and DELE-L1, we consider τ min = 0.25 ns and τ max = 15 ns for the library of exponentials in order to cover the characteristic lifetimes of the synthetic datasets in Table 1.For these datasets, we propose five evaluation scenarios in order to identify the effects of the parameters in DELE-L2 and DELE-L1, and compare their performance with respect to the metrics (GoF Y , GoF H , E Y , E H , GoF LT ), as well as their computational times: Case I : the effect of the number of elements in the library of exponentials N will be analyzed by varying this value, i.e.N ∈ {25, 50, 75, 100, 125}.In this evaluation, the regularization weight is keep constant at α = 0.001.Furthermore, we make the evaluation under different conditions of both shot and Gaussian noise by considering the following value pairs (PSNR,SNR) ∈ { (15,45), (20, 50), (25, 55), (30, 60)} dB.
Case IV : the performance of DELE-L2, DELE-L1, DELB, DENLS and DEGNLS is next evaluated with a synthetic dataset with short and heterogeneous average lifetimes in the sample at the following parameters PSNR=20 dB, SNR=50 dB, N = 25 and α = 0.001.In this first evaluation scenario, the purpose is to study the performance of DELE-L2 and DELE-L1 by varying the size of the library of exponentials at different pair levels of shot and Gaussian noise.Since the performance of DELE-L2 and DELE-L1 was similar, we present just the results for DELE-L2.Figure 4 shows the results for the three synthetic datasets.Keeping in mind that better fitting with respect to the variability of the information will be represented by values of (GoF Y , GoF H ) closer to 1.0, and more accurate estimations are highlighted by smaller values in (E Y , E H , GoF LT ), we observe in Fig. 4 that as the noise levels decrease (i.e.PSNR and SNR increase), (GoF Y , GoF H ) increase and (E Y , E H , GoF LT ) decrease for all datasets.Consequently, as expected, better signal strength with respect to the uncertainty sources in the dataset will imply more precise estimations by our deconvolution methodology.However, the number of elements in the library of exponentials N does not affect significantly the estimation performance.To look in more detail at this behavior, Table 3 shows the results of GoF LT as a function of N for PSNR=20 dB and SNR=50 dB.This table illustrates that in general an increment in N shows just a slight decrease in GoF LT , since the resolution in the library of exponential is increased.Nonetheless, the parameter N largely affects the optimization complexity, where as N increases, the resulting computational time will raise.In addition, the computational time is also raised as PSNR and SNR increase.This last tendency is not contradictory, since we are dealing with an inverse problem in the deconvolution process, so the noise in the dataset helps to avoid flat surfaces in the estimation error cost function to improve convergence.To explore further the reduction in N with respect to the estimation accuracy, we compute the evaluation for N = 10, but in this case, we indeed observe a difference with lower values in GoF Y and GoF H , and larger errors in E H and E Y .Hence we did not find a modification in the estimation performance for N ≥ 25; although the computational time is largely incremented for N > 25.
For the subsequent evaluations, we define as N = 25 the number of elements in the library of exponentials, to balance the computational load and accuracy.3. Case I: Estimation performance in GoF LT for PSNR=20 dB and SNR=50 dB in the synthetic datasets (A), (B) and (C) of 2nd, 3rd and 4th order, respectively.

Dataset
GoF LT (%)  In this new scenario of the synthetic validation, we consider a comparison among DELE-L2, DELE-L1, DELB, DENLS and DEGNLS at PSNR=20 dB, SNR=50 dB.Based on the previous two cases, for DELE-L2 and DELE-L1, we define their parameters as N = 25 and evaluate two conditions α = 0.001 and 0.01.Meanwhile, for DELB and based on trial-and-error, we choose an 8th order approximation in the Laguerre expansion and a shape parameter of 0.85 [19].Finally, for DENLS and DEGNLS, we set the order of the multi-exponential model as 2nd, 3rd and 4th for the three synthetic datasets, and the interval of variations of the characteristic lifetimes as [0.25, 15] ns.The performance results are presented in Table 4 for the three synthetic datasets.As mentioned earlier, DEGNLS has the fastest computational time of all schemes, but at the cost of lower estimation accuracy.Hence consistently DELE-L2 and DELE-L1 outperformed DENLS and DEGNLS in the indexes (GoF Y , GoF H , E Y , E H , GoF LT ) for all datasets; and with respect to DELB, DELE-L2 and DELE-L1 were more accurate in most of the tested conditions.For the two cases α = 0.001 and 0.01, DELE-L2 did not show noticeable differences in the accuracy, although there was a considerable reduction in computational time with α = 0.001.However, for DELE-L1, the performance was roughly the same for either value of α, but with a more higher processing time than DELE-L2 and DELB.Consequently, in term of accuracy and computational time, we see a comparable performance between DELE-L2 and  6 shows the ground-truth for the average lifetime maps in the three synthetic datasets, as well as the estimations by the tested methodologies DELE-L2 and DELE-L1 for α = 0.001, DELB, DENLS and DEGNLS.In these plots, we observe that DELE-L2, DELE-L1 and DELB have a good agreement with the ground-truth, and DENLS and DEGNLS present just some underestimation at some regions, which is consistent with the reported index GoF LT in Table 4.In this next evaluation with synthetic data, we generate a sample with short and heterogeneous average lifetimes by considering the 3rd order model in dataset (B) with characteristic lifetimes τ 1 = 225 ps, τ 2 = 450 ps, and τ 3 = 675 ps.The sampling time in this case has to be reduced with respect to Cases I to III to address this scenario, and we set it to T = 45 ps.The parameters of the mixed noise model in Eq. ( 13) are PSNR=20 dB and SNR=50 dB.Since the expected average lifetimes are lower compared to Cases I to III, we specify their range of variations as [45, 2500] ps in DELE-L2, DELE-L1, DENLS and DEGNLS.In DELB, we maintain an 8-th order Laguerre expansion but we reduced its shape parameter to 0.8 by trial-and-error to improve the fitting performance.For DELE-L2 and DELE-L1, we take the same parameters as in Case III: N = 25 and α = 0.001.The resulting average lifetime maps are illustrated in The first experimental evaluation considers the analysis of multi-spectral FLIM (m-FLIM) datasets of ex-vivo atherosclerotic plaques [6], where human coronary artery segments were obtained from 8 autopsy cases within 48 h of the time of death, according to a protocol approved by the Texas A&M University Institutional Review Board.Each arterial segment was longitudinally opened and imaged from the lumen side with a temporal resolution of 250 ps.The imaged field of view (FOV) had an area of 2 mm × 2 mm at 60 × 60 pixels.All the measurements include three wavelength bands: 390 ± 20 nm, 452 ± 22.5 nm and 550 ± 20 nm.Each effective time trace has a total of L = 167 samples per wavelength.Sixty individual datasets were analysed by DELE-L2, DELE-L1, DELB, DENLS and DEGNLS algorithms to estimate the fluorescence decays and impulse responses.As in the synthetic case, for DELE-L2, DELE-L1, DENLS and DEGNLS, we consider τ min = 0.25 ns and τ max = 15 ns; and for DELEs methodologies, we set N = 25.Meanwhile, the parameters of DELB were also selected as in Case III of the synthetic test.The m-FLIM measurements from atherosclerothic plaques contain mostly only 3 fluorophores: collagen, elastin and low density lipoproteins.The first two present a similar response in the 390 nm and 452 nm bands.Therefore, we analyze only the third wavelength which provides distinct information for classification purposes, i.e. a total of K = 3, 600 spatial points for deconvolution.For all datasets, the instrument response has the same time-domain characteristics with a FWHM of 1.525 ns.Since we do not have a ground-truth to compare the estimated fluorescence impulse responses, the evaluation is only carried out with respect to the accuracy in the estimated fluorescence decays by indexes GoF Y and E Y in Eq. ( 15), and also the computational time.Figure 9 shows the boxplot graphs for indexes GoF Y and E Y , and the analysis of the processing time over the 60 datasets for DELE-L2, DELE-L1, DELB, DENLS and DEGNLS algorithms.As

Human breast cancer cell datasets
The last experimental evaluation consists on processing fluorescence lifetime images of 13 human breast cancer cell samples used to measure metabolic inhibition by cyanide treatment [12].All cell lines were acquired from the American Type Culture Collection.The cancerous mammary epithelium cell lines (MCF10A) were grown in DMEM (Invitrogen) with 10% FBS and 1% penicillin: streptomycin.For fluorescence imaging, cells were plated at a density of 10 6 cells per 35 mm glass-bottom imaging dish (MatTek Corp.) 48 hours before imaging.Media of the MCF10A dishes was removed and replaced with cyanide supplemented MCF10A growth media (4 mmol/L NaCN, Sigma).The cells were allowed 5 minutes for the cyanide to react, and post-cyanide fluorescence images were acquired.The fluorescence images were captured by a custom-built multiphoton microscope (Prairie Technologies).Since in the samples, the primary target is to quantify the fluorescence response of coenzymes NADH and FAD, there are two optic excitations: one tuned to 750 nm for NADH excitation, and another one to 890 nm for FAD excitation.To isolate the response to NADH and FAD, there are two bandpass filters: 440 ± 40 nm (NADH) and 550 ± 50 nm (FAD).Each pixel has a dwell time of 4.8 μs to acquire images of dimension 256 × 256.The fluorescence lifetime image for each optic excitation was collected using time-correlated single-photon counting electronics, where the instrument response FWHM was 260 picoseconds for both laser excitations.In this testbench, the temporal resolution of the measurements is 49 ps, and the length of the time responses is 175 and 190 samples for the 1st and 2nd channels, respectively.The expected characteristic lifetimes are smaller in this application compared to the atherosclerotic plaques datasets, and the measured fluorescent decays have a strong noise component.Hence for DELE-L2, DELE-L1, DENLS and DEGNLS, we consider now τ min = 0.25 ns and τ max = 7.5 ns; and for DELE-L2 and DELE-L1, we set again N = 25.Once more and based on trial-and-error, an 8-th order is used in the Laguerre expansion for DELB with a shape parameter of 0.85.We evaluate both channels for the 13 datasets.The performance results are once more just evaluated with respect to the estimation of the fluorescent decays by GoF Y and E Y , and also by the computational time.Figure 11 illustrates the same tendency of the initial test with atherosclerotic plaques datasets, where DENLS is clearly improved by DELE-L2, DELE-L1 and DELB in terms of better fitting and less computational time, and DEGNLS is overcome in estimation accuracy.The mean values for GoF Y are 0.7741, 0.7143 and 0.7691 for DELE-L2, DELE-L1 and DELB, respectively, and for E Y are 0.0218, 0.0342 and 0.0235.Meanwhile, the processing times are 4.0349 s, 19.5209 s and 4.2709 s for DELE-L2, DELE-L1 and DELB, respectively.Consequently, now DELE-L2 presents a slight better fitting performance and less computational time with respect to DELB, but a much lower processing time compared to DELE-L1.Finally, the resulting average lifetime mappings are presented in Fig. 12 for datasets 4 (channel 2) and 8 (channel 1), where in this scenario, we observe an overestimation by DENLS and DEGNLS compared to DELE-L2, DELE-L1 and DELB.Consequently, this evaluation corroborates that DELE-L2 shows a good compromise between estimation performance and processing time over the other deconvolution strategies.

Conclusions
In this paper, we have proposed two new deconvolution strategies based on a library of exponentials and regularized non-negative least squares approximations (DELE-L2 and DELE-L1).
The range of variations of the characteristic lifetimes in the library is selected based on a predefined knowledge of the studied sample.For the regularization terms, we suggest Thikonov or l 1 weights in the optimization strategies.We extensively tested our proposal under synthetic datasets by considering both shot and Gaussian noise and samples with different lifetime maps, and in the experimental evaluation, we employed two types of datasets: ex-vivo atherosclerotic plaques and human breast cancer cells.We concluded that our approach with the Thikonov regularization weight (DELE-L2) exhibits a good compromise between estimation accuracy, computational time and tuning strategy compared to the previous methodologies from the literature.As a part of the future work, we will pursue to develop a blind deconvolution strategy based on the library of exponentials perspective that could simultaneously estimate the instrument and fluorescence impulse responses.We will also derive a toolbox and graphical user interface in Matlab to simplify the use of our proposed algorithms.

Fig. 4 .
Fig. 4. Case I: Validation test of the number of elements in the library of exponential functions N ∈ {25, 50, 75, 100, 125} at different pairs of PSNRs and SNRs in the synthetic datasets (A) (solid lines), (B) (dashed lines) and (C) (cross and dotted lines) of 2nd, 3rd and 4th order, respectively.

Fig. 10 .
Fig. 10.Estimated average lifetime maps in ex-vivo atherosclerotic plaques datasets by different deconvolution techniques.

Fig. 11 .
Fig. 11.Estimation performance in human breast cancer cell datasets by different deconvolution techniques.

Fig. 12 .
Fig. 12.Estimated average lifetime maps in human breast cancer cell datasets by different deconvolution techniques.

Table 1 .
List of Acronyms.