Deconvolution of fluorescence lifetime imaging microscopy by a library of exponentials ∗

Daniel U. Campos-Delgado,1, O. Gutierrez Navarro,1 E. R. Arce-Santana,1 Alex J. Walsh,2 Melissa C. Skala,2 and Javier A. Jo3 1

Facultad de Ciencias, Universidad Autonoma de San Luis Potosi, SLP, Mexico of Biomedical Engineering, Vanderbilt University, Nashville, Tennessee, USA 3 Department of Biomedical Engineering, Texas A& M University, College Station, Texas, USA 2 Department

∗ [email protected]

Abstract: 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 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 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 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 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.

1.

Introduction

Fluorescence is an important resource for quantitative characterization of samples [1–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–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–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–22], or analyze the datasets in a lower-dimensional domain by the phasor perspective [23–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–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–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–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/l2 or l1 regularization term [26–29], where this process is parallelized to speed up convergence. We present an advantage of our approach based on Thikonov/l2 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–15], and a constrained Laguerre-

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23750

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, Rn n-dimensional real vectors, Rn×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 1n (0n ), and In denotes the identity matrix of order n. For a vector x ∈ Rn , diag(x) ∈ Rn×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. Table 1. List of Acronyms.

DEGNLS DELB DELE DENLS FWHM FLIM FOV NNLSA PSNR SNR

2.

Deconvolution Estimation by Global Nonlinear Least Squares Deconvolution Estimation by Laguerre Basis Deconvolution Estimation by a Library of Exponentials Deconvolution Estimation by Nonlinear Least Squares Full Width at Half Maximum Fluorescence Lifetime Imaging Microscopy Field of View Non-Negative Least Squares Approximation Peak to Signal Noise Ratio Signal to Noise Ratio

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–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 yk [l] = hk [l]  u[l] + vk [l] =

L−1

∑ u[l − j] hk [ j] + vk [l]

∀l ∈ [0, L − 1], k ∈ [0, K − 1]

(1)

j=0

where yk [l] and hk [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 vk [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 Ak > 0 of the instrument response u[l] in the observation model of Eq. (1) [1– 3], i.e. under the effect of scattering we have now yk [l] = (hk [l]+Ak δ [l])u[l]+vk [l], where δ [l] denotes the discrete-time delta function. Although this component could potentially introduce a distortion of the estimated impulse response function, that is (hk [l] + Ak δ [l]) instead of hk [l], 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

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23751

Eq. (1) can be written in a vector notation as ⎤ ⎡ ⎡ u[0] 0 ... yk [0] ⎢ yk [1] ⎥ ⎢ u[1] u[0] ... ⎥ ⎢ ⎢ ⎥=⎢ ⎢ .. .. .. .. ⎦ ⎣ ⎣ . . . . 

yk [L − 1]





u[L − 1] u[L − 2]

yk ∈RL

...

0 0 .. . u[0]

⎤⎡ ⎥⎢ ⎥⎢ ⎥⎢ ⎦⎣ 

hk [0] hk [1] .. .

⎤

⎥ ⎢ ⎥ ⎢ ⎥+⎢ ⎦ ⎣

⇒ yk = U hk + vk

vk [0] vk [1] .. .

⎤ ⎥ ⎥ ⎥, ⎦

(2)

hk [L − 1] vk [L − 1]

 

 hk ∈RL

U∈RL×L

⎡

vk ∈RL

∀k ∈ [0, K − 1],

(3)

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 = {y1 , . . . , yK }. Hence the deconvolution problem estimates for the impulse response vector hˆ k at each spatial point, such that the estimated fluorescence decay vector yˆ k = U hˆ k approximates the measurement one yk [1–3]. As any inverse problem, the deconvolution is an ill-posed formulation that can have

Fig. 1. FLIM observation model.

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. L−1

∑ u[l] = 1

&

u[l] ≥ 0

∀l ∈ [0, L − 1],

(4)

l=0

B. The fluorescence impulse response hk [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 hk [l] is modeled as: hk [l] = ck,1 e−l/τ1 + ck,2 e−l/τ2 + . . . + ck,N e−l/τN + ck,N+1 #243670 (C) 2015 OSA

∀l ∈ [0, L − 1],

(5)

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23752

where the non-negative positive coefficient ck,n ≥ 0 represents the contribution of the n-th exponential mode into the fluorescence response at k-th spatial point, and ck,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 yk in the dataset is assumed to belong to the convex cone Ω = {y ∈ RL | y =

N+1

∑ cn U bn

∀cn ∈ R & cn ≥ 0}

(6)

n=1

where bn = [1 e−1/τn . . . e−(L−1)/τn ] ∈ RL for n ∈ [1, N] and bN+1 = 1L , i.e. Y ⊂ Ω. In this way, the L-dimensional vectors {Ubn }Nn=1 will represent a library of fluorescence decays that will be used to approximate each measurement in the FLIM dataset, see Fig. 2(b).

Library of Exponential Functions

Assumption A is important to avoid numerical scaling problems in the estimation of the coefficients {ck,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 {ck,n }Nn=1 , so we will like to induce some regularization into the solution. As a result, we propose to use a Thikonov/l2 or l1 regularization method for this purpose [26–29]. We see (a) 1 0.5 0 15 10 5 τn (nsec)

10

5

0

15

20

25

30

35

40

25

30

35

40

Library of Fluorescence Decays

Time (nsec) (b) 1

0.5

0 15 10 5 τ (nsec) n

0

5

10

15

20

Time (nsec)

Fig. 2. Resulting libraries for Ts = 250 ps, τmin = 0.25 ns, τmax = 15 ns and N = 25: (a) exponential functions and (b) fluorescence decays.

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– 15], our strategy is linear in the unknown parameters {ck,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–19], our methodology guarantees a smooth fluorescent impulse response with a monotonic decreasing pattern, and also, the order of #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23753

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 }Nn=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 yk = UBck + vk

∀k ∈ [0, K − 1].

(7)

where B = [b1 . . . bN+1 ] ∈ RL×(N+1)

 ck = ck,1 . . . ck,N+1 ∈ RN+1 ,

(8) (9)

i.e. matrix B contains in its columns the information of the library of exponentials, and vector ck its unknown scaling coefficients. Note that a direct least squares approximation from Eq. (7) to compute ck could be really fast but not feasible, since some of the elements in ck can be negative, which will not represent a meaningful solution for the deconvolution problem. 3.

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 ck are estimated for each spatial point k ∈ [0, K − 1] by a non-negative least-squares approximation (NNLSA) [32] with a Thikonov/l2 regularization term as follows  ⎡ ⎤  2  U B     1 y k  , ⎣ √ ⎦ ck − 1N ∑ (ck,n )2 = cmin   0 α diag 0 2 N+1 k   n=1 0 (10) where α > 0 is the regularization weight. An important property of Eq. (10) is that the coefficients vector ck 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 l1 regularized least squares approximation for the deconvolution process [26–28]: 1 α min yk − UBck 2 + ck 0 2 2

N

1 α min yk − UBck 2 + ck 0 2 2

N

1

α

yk − UBck 2 + [1 ∑ ck,n = cmin 2 N k 0 2

0]ck .

(11)

n=1

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 l1 regularization term in Eq. (11).

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23754

4.

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–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 l1 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–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]. 4.1.

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]: 2,3 or 4

hk [l] =

∑

−l τT

δk,i e

i

+ δo

∀k ∈ [0, K − 1], l ∈ [0, L − 1],

(12)

i=1

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 #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23755

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 yok [l] are computed by applying the convolution operator in Eq. (1), i.e. yok [l] = u[l]  hk [l]. 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–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]:  (13) yk [l] = yok [l] + yok [l] · ωk [l] + vk [l] ∀l ∈ [0, L − 1] where ωk [l] ∼ N (0, σk2 ) and vk [l] ∼ N (0, ςk2 ) represent normal random variables related to the shot and Gaussian noise terms, respectively, where the variances σk2 and ςk2 are selected with respect to desired peak signal-to-noise (PSNR) and signal-to-noise ratios (SNR): PSNR = 10 log10 SNR = 10 log10

maxl∈[0,L−1] yok [l]

σk2

∀k ∈ [0, K − 1],

o 2 ∑L−1 l=0 (yk [l]) ςk2

(14)

In this way, for our model  in Eq. (13), the standard deviation of the shot noise component is time-varying and equal to yok [l] · σk , which is adjusted by σk to achieve the desired PSNR at k-th pixel. By recalling that hk and yk denote the true vectors for the impulse response and fluorescent decays, respectively, and hˆ k and yˆ 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–3, 20]:   hk − hˆ k  1 K−1 GoFH = ∑ 1 − hk − h ¯ K k=0   yk − yˆ k  1 K−1 GoFY = ∑ 1 − yk − y¯  K k=0 EH =

ˆ 2 ∑K−1 k=0 hk − hk  2 ∑K−1 k=0 hk 

EY =

ˆ k 2 ∑K−1 k=0 yk − y K−1 ∑k=0 yk 2

(15)

¯ = (1/K) ∑K−1 where h¯ = (1/K) ∑K−1 k=0 hk and y k=0 yk are the mean impulse response and mean fluorescent decay in the dataset, respectively. Hence the indexes GoFH and GoFY evaluate the fitting accuracy with respect to the variability in the dataset of the impulse responses and fluorescence decays, respectively; meanwhile EH and EY 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

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23756

responses [20] at each spatial point k:

τk =

t hk 1 L hk

& τˆk =

t hˆ k , ˆ 1 L hk

(16)

where t = [0 Ts 2Ts . . . (L − 1)Ts ] represents a vector of sampling instants. Hence we define the goodness of fit with respect to the percentage average lifetime estimation error as GoFLT = 100 ×

1 K−1 |τk − τˆk | ∑ τk . K k=0

(17)

Table 2. Characteristic lifetimes in synthetic datasets (A), (B) and (C) of 2nd, 3rd and 4th order, respectively.

Dataset (A) (B) (C)

Characteristic Lifetimes τ1 τ2 τ3 τ4 3.75 ns 8.75 ns – – 8.75 ns 3.75 ns 1.25 ns – 1.25 ns 3.75 ns 6.25 ns 8.75 ns

δk,2

δk,1

Dataset 20 40 (A) 60 80 100

20 40 60 80 100 0

0.5

20 40 60 80 100

1

20 40 60 80 100

0

20 40 60 80 100 0

0.5

1

0

20

20 40 60 80 100 0.5

0.5

1

20 40 60 80 100 0

1

20 40 60 80 100 20 40 60 80 100

20 40 60 80 100 0

0.5

0.5

1 δk,4

δk,3

δk,2

Dataset 40 (C) 60

0

δk,3

20 40 60 80 100 20 40 60 80 100

20 40 60 80 100

δk,1

80 100

1

k,2

Dataset 20 (B) 40 60 80 100

0.5 δ

δk,1

1

20 40 60 80 100 0

0.5

1

20 40 60 80 100

20 40 60 80 100 0

0.5

1

Fig. 3. Abundance maps for synthetic datasets (A), (B) and (C) of 2nd, 3rd, and 4th order, respectively.

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 (GoFY , GoFH , EY , EH , GoFLT ), as well as their computational times: #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23757

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 II : the effect of the regularization weight α is now studied by varying this term, i.e. α ∈ {0.0001, 0.0005, 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1.0}. In this test, we set the number of elements in the library of exponentials as N = 25 with PSNR=20 dB and SNR=50 dB. Case III : the performance of DELE-L2, DELE-L1, DELB, DENLS and DEGNLS is evaluated for the three synthetic datasets at the following parameters PSNR=20 dB, SNR=50 dB, N = 25 and α = 0.001. 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. Case V : the effect of photon economy in the synthetic FLIM datasets is lastly evaluated for DELE-L2, DELE-L1, DELB, DENLS and DEGNLS by taking into account different value pairs (PSNR,SNR) ∈ {(15, 45), (20, 50), (25, 55), (30, 60)} dB in the three synthetic datasets, with the following parameters N = 25 and α = 0.001. 4.1.1.

Case I

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 (GoFY , GoFH ) closer to 1.0, and more accurate estimations are highlighted by smaller values in (EY , EH , GoFLT ), we observe in Fig. 4 that as the noise levels decrease (i.e. PSNR and SNR increase), (GoFY , GoFH ) increase and (EY , EH , GoFLT ) 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 GoFLT 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 GoFLT , 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 GoFY and GoFH , and larger errors in EH and EY . 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.

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23758

1

1

GoF

GoF

Y

H

0.98 0.96 0.94

Datasets (B) & (C), N=25,...,125

0.92 45 15 −3 x 10 4

Dataset (A), N=25,...,125

50

55

60

0.9 45

20

25

30

15

50

55

60

20

25

30

0.01 Dataset (A), N=25,...,125

Dataset (C), N=25,...,125 EH

EY

3

0.95

2

0.005

1

GoFLT (%)

15 15

50

55

60

20

25

30

Dataset (B), N=25,...,125 10 5 0 45

50 SNR (dB) 55

60

15

20

30

PSNR (dB)

25

0 45 Computational Time (sec)

0 45

55

60

15 20 25 60 Datasets (A) & (B), N=125

30

40

50

Datasets (A) & (B), N=25

20 0 45

50 SNR (dB) 55

60

15

20 PSNR (dB) 25

30

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.

Table 3. Case I: Estimation performance in GoFLT for PSNR=20 dB and SNR=50 dB in the synthetic datasets (A), (B) and (C) of 2nd, 3rd and 4th order, respectively.

Dataset (A) (B) (C)

4.1.2.

N = 25 3.4663 3.9569 6.0763

N = 50 3.4660 3.9557 6.0757

GoFLT (%) N = 75 N = 100 3.4652 3.4652 3.9554 3.9553 6.0769 6.0771

N = 125 3.4650 3.9552 6.0777

Case II

We analyze now the effect of the regularization weight α ∈ [0.0001, 1.0] in Eq. (10) on the estimation indexes (GoFY , GoFH , EY , EH , GoFLT ) and the resulting computational time. As in Case I, we consider the results just for DELE-L2. Figure 5 presents the achieved performance with fixed values for PSNR=20 dB, SNR=50 dB and N = 25. For the indexes (GoFY , GoFH , EY , EH ), the best performance is achieved at α = 0.1 but at the cost of a high computational time. Actually, the results in GoFLT show that lower values of α could provide good estimations of the average lifetime, at the price of a modest reduced accuracy illustrated in (GoFY , GoFH , EY , EH ). Therefore looking to balance the tradeoff between estimation accuracy and computational time,

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23759

we chose this weight as either α = 0.001 or 0.01 for the next evaluation. 0.98

0.96

0.96

GoF

Y

GoFH

0.98

0.94 0.92

0.94 0.92

0.9 −4 10

−2

10

0

10

0.9 −4 10

−3

x 10

2 0 −4 10

−2

10

4

LT

6 4 −2

10 α

0 −4 10

0

10

8 GoF

x 10

2

10

2 −4 10

0

10

6

Dataset (A) Dataset (B) Dataset (C)

EH

EY

4

8

0

10

Computational Time (sec)

6

−2

10 −3

−2

10

0

10

6 5 4

Dataset (A) Dataset (B) Dataset (C)

3 2 −4 10

−2

10 α

0

10

Fig. 5. Case II: Validation test of weight parameter α ∈ [0.0001, 1.0] for PSNR=20 dB, SNR=50 dB, and N = 25 with synthetic datasets (A) (solid line), (B) (dashed line) and (C) (cross and dotted line) of 2nd, 3rd and 4th order, respectively.

4.1.3.

Case III

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 (GoFY , GoFH , EY , EH , GoFLT ) 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

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23760

DELB, where DELE-L2 has the advantage of a more intuitive tuning of its parameters. Figure 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 GoFLT in Table 4. Table 4. Case III: Performance Quantification of Synthetic Datasets for DELE-L2, DELEL1, DELB, DENLS and DEGNLS at PSNR=20 dB, SNR=50 dB, and N = 25.

Dataset

α = 0.01 DELE-L2 DELE-L1

GoFY (A) 0.9669 0.9684 (B) 0.9626 0.9615 (C) 0.9626 0.9619 GoFH (A) 0.9583 0.9659 (B) 0.9469 0.9551 (C) 0.9464 0.9544 EY (A) 0.0008 0.0008 (B) 0.0010 0.0010 (C) 0.0011 0.0011 EH (A) 0.0016 0.0010 (B) 0.0024 0.0016 (C) 0.0027 0.0018 GoFLT (%) (A) 3.4216 3.5947 (B) 3.9294 5.0459 (C) 6.0569 6.8342 Computational Time (sec) (A) 3.1501 10.4177 (B) 3.5383 10.8118 (C) 3.0604 10.9490

4.1.4.

Deconvolution Technique α = 0.001 DELE-L2 DELE-L1 DELB

DENLS

DEGNLS

0.9662 0.9615 0.9613

0.9679 0.9629 0.9630

0.9624 0.9591 0.9594

0.9469 0.9087 0.9422

0.9363 0.9055 0.9288

0.9560 0.9433 0.9422

0.9641 0.9535 0.9516

0.9555 0.9493 0.9492

0.9442 0.8836 0.9277

0.9326 0.8796 0.9096

0.0009 0.0010 0.0012

0.0008 0.0010 0.0011

0.0010 0.0012 0.0013

0.0019 0.0047 0.0021

0.0027 0.0050 0.0030

0.0018 0.0029 0.0033

0.0011 0.0017 0.0020

0.0017 0.0021 0.0023

0.0024 0.0098 0.0041

0.0034 0.0104 0.0064

3.4910 3.9443 5.9802

3.6027 4.6718 6.4858

3.5371 4.0819 6.5618

6.8582 6.5580 10.1308

8.1836 6.7262 12.8169

2.6402 2.5629 2.3637

10.3692 12.9724 12.6456

1.3861 1.3699 1.4645

27.8680 28.8734 30.6140

0.5934 0.6505 0.7004

Case IV

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 #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23761

Ground−Truth

Dataset (A)

20 40 60 80 100

DELE−L2

6

8

(ns) 4

Ground−Truth 20 40 60 80 100

6

8

DELE−L2

(ns)

Ground−Truth 20 40 60 80 100

DENLS

DEGNLS

4

6

8

(ns) 4

DELE−L1

6

8

DELB

4

6

8

DENLS

4

6

8

(ns)

DEGNLS

20 20 20 20 20 40 40 40 40 40 60 60 60 60 60 80 80 80 80 80 100 100 100 100 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 4 6 8

Dataset (C)

DELB

20 20 20 20 20 40 40 40 40 40 60 60 60 60 60 80 80 80 80 80 100 100 100 100 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 4

Dataset (B)

DELE−L1

4 6 8

4 6 8

DELE−L2

DELE−L1

(ns)

4 6 8

4 6 8

4 6 8

DELB

DENLS

DEGNLS

(ns)

20 20 20 20 20 40 40 40 40 40 60 60 60 60 60 80 80 80 80 80 100 100 100 100 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 20 60 100 2 4 6

(ns) 2 4 6

2 4 6

(ns) 2 4 6

2 4 6

2 4 6

(ns)

Fig. 6. Case III: Validation test of different deconvolution techniques by computing average lifetime maps at PSNR=20 dB, SNR=50 dB, N = 25 and α = 0.001 with synthetic datasets (A), (B) and (C) of 2nd, 3rd, and 4th order, respectively.

Fig. 7, where we observe a more precise characterization with respect to the ground-truth by DELE-L2, DELE-L1 and DELB, compared to DENLS and DEGNLS. In fact, this visual appreciation was corroborated by the index GoFLT which was 8.99 %, 8.58 %, 9.78 %, 12.88 % and 13.67 % by DELE-L2, DELE-L1, DELB, DENLS and DEGNLS, respectively; with the following computational times 2.36 s (DELE-L2), 11.20 s (DELE-L1), 1.64 s (DELB), 28.26 s (DENLS) and 0.62 s (DEGNLS). The rest of the performance indexes (GoFY , GoFH , EY , EH ) showed the same trend described by GoFLT . Hence although all the deconvolution strategies were able to estimate the short average lifetimes in the sample with errors lower than 14 %, we see that DELE-L2 presents the best compromise between accuracy and computational time of all schemes. 4.1.5.

Case V

Since photon economy is strictly related to the signal SNR, we vary the uncertainty in the measured photon counts by different pairs of (PSNR,SNR) ∈ {(15, 45), (20, 50), (25, 55), (30, 60)} dB for the three synthetic datasets (A),(B) and (C) in Fig. 3, and compute the deconvolution process by DELE-L2, DELE-L1, DELB, DENLS and DEGNLS. We set the following parameters in these algorithms: N = 25, α = 0.001, order of Laguerre expansion 8th and its shape parameter 0.85, and interval of variations of the characteristic lifetimes as [0.25, 15] ns. For brevity, we just show the results of the performance indexes GoFY and GoFLT , which are defined in Eq. (15). By recalling that a value of GoFY close to 1.0 indicates better fitting to the fluorescent

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23762

Ground−Truth

DEEL−L2

DEEL−L1

DELB

DENLS

DEGNLS

20

20

20

20

20

20

40

40

40

40

40

40

60

60

60

60

60

60

80

80

80

80

80

80

100

20

60

100

100

20

0.7 (ns)

0.35

60

0.35

100

100

20

0.7 (ns)

60

0.35

100

100

20

0.7 (ns)

60

0.35

100

100

20

0.7 (ns)

60

0.35

100

100

20

0.7 (ns)

60

0.35

100 0.7 (ns)

Fig. 7. Case IV: Validation test of different deconvolution techniques by computing average lifetime maps in a sample with short characteristic lifetimes at PSNR=20 dB, SNR=50 dB, N = 25 and α = 0.001.

decays in the dataset, and GoFLT close to 0 better estimation of the average lifetimes, Fig. 8 illustrates that our proposed deconvolution schemes DELE-L2 and DELE-L1 achieve the best responses for the three datasets and tested pairs SNR/PNSR. Meanwhile, DENLS and DEGNLS show the less accurate fitting. In fact, the four strategies DELE-L2, DELE-L1, DENLS and DELB present a consistent pattern, where as the uncertainty is lowered (i.e. the pair SNR/PSNR increases), the fitting performance improves. However, DEGNLS could exhibit some inconsistencies, which are inherited by the nonlinear nature of the optimization process involved in its solution. Dataset (A)

Dataset (C)

Dataset (B)

0.99

1

0.98

0.98

0.97

0.96

0.99 0.98

GoFY

GoFY

0.97

0.96 0.95

0.96 0.95

0.94

DELE−L2

0.94

0.92

DELE−L1

0.93 0.94 0.93 45

0.9

50

55

60

0.88 45

50

20

DENL

55

0.91 45

60

DEGNL

50

55

60

20 25 PSNR (dB)

30

SNR (dB)

SNR (dB) 15

DELB

0.92

25 PSNR (dB)

30

15

9

9

8

8

7

7

6

6

SNR (dB)

20 25 PSNR (dB)

30

15 20

DELE−L2 DELE−L1 DELB

15

5

LT

5

4

4

3

3

2

2

1 45

50

55

60

1 45

10

5

50

55

60

0 45

30

15

SNR (dB)

SNR (dB) 15

DEGNL

GoF

GoFLT

GoFLT

DENL

20 25 PSNR (dB)

30

15

20 25 PSNR (dB)

55

60

20 25 PSNR (dB)

50

30

SNR (dB)

Fig. 8. Case V: Validation test for photon economy in the synthetic FLIM datasets by using different pairs of PSNRs and SNRs at N = 25 and α = 0.001.

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23763

4.2. 4.2.1.

Experimental evaluation Ex-vivo atherosclerotic plaques datasets

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 DELEL2, 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.

GoFY

0.9 0.85 0.8 0.75 DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

0.04

EY

0.03 0.02

Computational Time (sec)

0.01

10

5

0

Fig. 9. Estimation performance in ex-vivo atherosclerotic plaques datasets by different deconvolution techniques.

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 GoFY and EY in Eq. (15), and also the computational time. Figure 9 shows the boxplot graphs for indexes GoFY and EY , and the analysis of the processing time over the 60 datasets for DELE-L2, DELE-L1, DELB, DENLS and DEGNLS algorithms. As

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23764

4

(ns)

20 40 60 pixel 5 6 7 8 9

20 40 60

(ns)

20 40 60 pixel 5

10

(ns)

20 40 60

8

20 40 60

20 40 60 pixel 5 6 7 8 9

DELE−L1

DELB

5

10

4

20 40 60

10

20 40 60

8

20 40 60 pixel 4

20 40 60

(ns)

8

20 40 60

5 6 7 8 9

DENLS

DEGNLS

10

(ns)

20 40 60 pixel

5 6 7 8 9

20 40 60 pixel 5

6

DEGNLS

20 40 60 pixel

(ns)

20 40 60 pixel 5

6

20 40 60

DENLS

5 6 7 8 9

20 40 60 pixel

(ns)

pixel

pixel 6

20 40 60 pixel

DELB

pixel

pixel

pixel

20 40 60

4

20 40 60 pixel

DELE−L2

Dataset 54

8

pixel

20 40 60

pixel

pixel

Dataset 24

6

DELE−L1

20 40 60

DEGNLS

pixel

8

pixel

6

DELE−L2

20 40 60 pixel

pixel

4

20 40 60 pixel

20 40 60

DENLS

pixel

20 40 60 pixel

20 40 60

DELB pixel

20 40 60

DELE−L1 pixel

Dataset 18

pixel

DELE−L2

20 40 60

(ns)

20 40 60 pixel 5

10

(ns)

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

expected, there is a consistency in the performance by GoFY and EY , where DELE-L2, DELEL1 and DELB clearly improve the estimation accuracy of DENLS and DEGNLS. Moreover, DENLS provided the largest computational time of the studied methodologies. In the overall, with respect to accuracy and computational time, the performance of DELE-L2 and DELB is similar, and DELE-L1 just had a larger processing time. The mean values for GoFY are 0.8752, 0.8689 and 0.8797 for DELE-L2, DELE-L1 and DELB, respectively, and for EY are 0.0106, 0.0116 and 0.0099. Meanwhile, the processing times are 0.9506 s, 4.0363 s and 0.6295 s for DELE-L2, DELE-L1 and DELB, respectively. Finally, Fig. 10 presents the resulting average lifetime for datasets 18, 24 and 54 with each deconvolution methodology. As expected, these plots illustrate good agreement in the resulting lifetime maps for DELE-L2, DELE-L1 and DELB. Meanwhile, DENLS and DEGNLS present some underestimation in some regions of the maps. Hence, as in the synthetic evaluation (Case III), a similar performance is visualized for DELE-L2 and DELB, although the tuning strategy for the parameters of DELE-L2 is more intuitive. 4.2.2.

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 106 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, #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23765

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.

GoFY

0.8 0.75 0.7 0.65 0.6

DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

DELE−L2

DELE−L1

DELB

DENLS

DEGNLS

0.08 EY

0.06 0.04

Computational Time (sec)

0.02

80 60 40 20 0

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

The performance results are once more just evaluated with respect to the estimation of the fluorescent decays by GoFY and EY , 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 GoFY are 0.7741, 0.7143 and 0.7691 for DELE-L2, DELE-L1 and DELB, respectively, and for EY 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 #243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23766

DEEL−L1

DEGNLS

100

100

100

100

100

150

150 200

200

250 50

150

50

1

150

50

(ns)

0

1

DEEL−L2

150

250

250

50

pixel

2

(ns)

0

1

DEEL−L1

150 200

250

250

pixel

2

150 200

250

250

pixel

0

150

pixel

50

pixel

50

pixel

50

200

150

250

50

pixel

2

(ns)

0

1

DELB

2

(ns)

0

100

100

100

100

150

200

200

200

200

250

250

250

250

50

150

250

50

pixel

0

1

150

250

50

pixel

2

(ns)

0

1

150

pixel

100

pixel

50

pixel

50

150

(ns)

0

1

150

250 50

150

250

50

pixel

2

(ns)

(ns)

2

200

250

pixel

2

250

DEGNLS

50

150

1

DENLS

50

150

150

pixel

50

pixel

pixel

DENLS

50

250

Dataset 8 Channel 1

DELB

50

pixel

pixel

DEEL−L2

Dataset 4 Channel 2

0

1

150

250

pixel

2

(ns)

0

1

2

(ns)

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

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. 5.

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 l1 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. Acknowledgments This study was supported by the American Heart Association (Beginning Grant-in-Aid Grant 0765102Y), the National Institutes of Health (NIH) (1R03CA191860, 1R01HL111361, 1R01CA138653), and a Fulbright-Garcia Robles grant.

#243670 (C) 2015 OSA

Received 25 Jun 2015; revised 8 Aug 2015; accepted 22 Aug 2015; published 2 Sep 2015 7 Sep 2015 | Vol. 23, No. 18 | DOI:10.1364/OE.23.023748 | OPTICS EXPRESS 23767