Estimation of the number of fluorescent end-members for quantitative analysis of multispectral FLIM data

Multispectral fluorescence lifetime imaging (m-FLIM) can potentially allow identifying the endogenous fluorophores present in biological tissue. Quantitative description of such data requires estimating the number of components in the sample, their characteristic fluorescent decays, and their relative contributions or abundances. Unfortunately, this inverse problem usually requires prior knowledge about the data, which is seldom available in biomedical applications. This work presents a new methodology to estimate the number of potential endogenous fluorophores present in biological tissue samples from time-domain m-FLIM data. Furthermore, a completely blind linear unmixing algorithm is proposed. The method was validated using both synthetic and experimental m-FLIM data. The experimental m-FLIM data include in-vivo measurements from healthy and cancerous hamster cheek-pouch epithelial tissue, and ex-vivo measurements from human coronary atherosclerotic plaques. The analysis of m-FLIM data from in-vivo hamster oral mucosa identified healthy from precancerous lesions, based on the relative concentration of their characteristic fluorophores. The algorithm also provided a better description of atherosclerotic plaques in term of their endogenous fluorophores. These results demonstrate the potential of this methodology to provide quantitative description of tissue biochemical composition. © 2014 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) Timeresolved imaging References and links 1. 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, 152–161 (2008). 2. K. S. M. Mycek and N. Nishioka, “Colonic polyp differentiation using time-resolved autofluorescence spectroscopy,” Gastrointest. endosc. 48, 390–394 (1998). 3. L. Marcu, “Fluorescence lifetime in cardiovascular diagnostics,” J. Biomed. Opt. 15, 011106 (2010). #208447 $15.00 USD Received 17 Mar 2014; revised 2 May 2014; accepted 2 May 2014; published 13 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.012255 | OPTICS EXPRESS 12255 4. J. A. Jo, B. Applegate, J. Park, S. Shrestha, P. Pande, I. Gimenez-Conti, and J. Brandon, “In vivo simultaneous morphological and biochemical optical imaging of oral epithelial cancer,” IEEE Trans. Biomed. Eng. 57, 2596– 2599 (2010). 5. J. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” IEEE J. Sel. Topics Appl. Earth Observ. 5, 354–379 (2012). 6. H. Xu and B. W. Rice, “In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique,” J. Biomed. Opt. 14, 064011 (2009). 7. S. S. Vogel, P. S. Blank, S. V. Koushik, and C. Thaler, “Spectral imaging and its use in the measurement of Förster resonance energy transfer in living cells,” in Fret and Flim Techniques. Laboratory Techniques in Biochemistry and Molecular Biology , T. Gadella, ed. 33, 351–394, (Elsevier, 2009). 8. J. G. Gert-Jan Kremers, Erik B. van Munster and J. Theodorus W. J. Gadella, “Quantitative lifetime unmixing of multiexponentially decaying fluorophores using single-frequency fluorescence lifetime imaging microscopy,” Biophys. J. 95, 378–389 (2008). 9. O. Gutierrez-Navarro, D. Campos-Delgado, E. Arce-Santana, M. Mendez, and J. Jo, “A fully constrained optimization method for time-resolved multi-spectral fluorescence lifetime imaging microscopy data unmixing,” IEEE Trans. Biomed. Eng. 60, 1711–1720 (2013). 10. 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, 606–617 (2014). 11. J. Harsanyi, W. Farrand, and C.-I. Chang, “Determining the number and identity of spectral endmembers: An integrated approach using Neyman-Pearson eigenthresholding and iterative constrained rms error minimization,” presented at the 9th Thematic Conference on Geological Remote Sensing, Pasadena, California, USA (1993). 12. C.-I. Chang and Q. Du, “Estimation of number of spectrally distinct signal sources in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens. 42, 608–619 (2004). 13. J. Bioucas-Dias and J. Nascimento, “Hyperspectral subspace identification,” IEEE Trans. Geosci. Remote Sens. 46, 2435–2445 (2008). 14. H. Akaike, “A new look at the statistical model identification,” IEEE Trans. Automat. Contr. 19, 716–723 (1974). 15. J. Rissanen, “Modeling by shortest data description,” Automatica 14, 465–471 (1978). 16. R. Heylen and P. Scheunders, “Hyperspectral intrinsic dimensionality estimation with nearest-neighbor distance ratios,” IEEE J. Sel. Topics Appl. Earth Observ. 6, 570–579 (2013). 17. C. Andreou and V. Karathanassi, “Estimation of the number of endmembers using robust outlier detection method,” IEEE J. Sel. Topics Appl. Earth Observ. 7, 247–256 (2014). 18. A. Ambikapathi, T.-H. Chan, C.-Y. Chi, and K. Keizer, “Hyperspectral data geometry-based estimation of number of endmembers using p -norm-based pure pixel identification algorithm,” IEEE Trans. Geosci. Remote Sens. 51, 2753–2769 (2013). 19. J. R. Lakowicz, Principles of Fluorescence Spectroscopy (Springer, 2006). 20. T. Zimmermann, “Spectral imaging and linear unmixing in light microscopy,” Molecular Biology 95, 245–265 (2005). 21. P. Vallotton, A. Phatak and M. Berman, “Spectral Imaging and Unmixing” in Fluorescence Applications in Biotechnology and Life Sciences, E. M. Goldys ed. ( Wiley-Blackwell, 2009). 22. M. Craig, “Minimum-volume transforms for remotely sensed data,” IEEE Trans. Geosci. Remote Sens. 32, 542– 552 (1994). 23. S. Boyd and L. Vandenberghe, Convex Optimization (Cambridge University Press, New York, 2004). 24. D. Manolakis and G. Shaw, “Detection algorithms for hyperspectral imaging applications,” IEEE Signal Process. Mag. 19, 29–43 (2002). 25. F. W. Young, J. de Leeuw, and Y. Takane, “Regression with qualitative and quantitative variables: An alternating least squares method with optimal features,” Psychometrika 41, 505–526 (1976). 26. B. C. Levy, Principles of Signal Detection and Parameter Estimation (Springer Publishing Company, 2008). 27. 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,” Biophys. J. 81, 1265–1274 (2001). 28. 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, 394–401 (2011). 29. 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, 046012 (2013). 30. Robert M . Clegg and A. Periasamy FLIM Microscopy in Biology and Medicine (Chapman and Hall/CRC, 2009). 31. 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, 172–181 (2011). 32. J. Nocedal and S. J. Wright, Numerical Optimization (Springer, 2000). #208447 $15.00 USD Received 17 Mar 2014; revised 2 May 2014; accepted 2 May 2014; published 13 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.012255 | OPTICS EXPRESS 12256 33. K. Arakawa, K. Isoda, T. Ito, K. Nakajima, T. Shibuya, and F. Ohsuzu, “Fluorescence analysis of biochemical constituents identifies atherosclerotic plaque with a thin fibrous cap,” Arterioscler. Thromb. Vasc. 22, 1002–1007 (2002).


Introduction
Fluorescence lifetime imaging (FLIM) is a powerful tool in biological and clinical sciences for characterization of tissue samples.Despite several studies showing the potential of multispectral FLIM (m-FLIM) as a promising clinical optical imaging modality, its utility has not yet fully demonstrated [1][2][3][4].One of the main reasons for this is the element of subjectivity involved in interpreting multispectral FLIM images, which is based on qualitative comparison of spectral intensity and lifetime values obtained from sample fluorescence time-resolved measurements with those of known fluorophores.Such comparison-based interpretation of FLIM images is satisfactory only when the fluorescence decay recorded at a pixel in a FLIM image can be attributed to a single fluorophore.However, in most practical applications, there are often more than one fluorescing species contributing to the bulk fluorescence signal.In such cases, it becomes imperative to develop a method that can quantitatively identify the number, the type and the relative abundance of fluorophores present in the sample.
The number of components present in a sample, their identification and relative concentration or abundance are three major problems, whose solutions are needed to provide meaningful descriptions for m-FLIM data.These problems are also present in other fields of application, such as in remote sensing where they are known as unmixing problems [5].Different unmixing schemes have been applied with spectral fluorescence data [6,7], frequency resolved FLIM data [8], and more recently with time-resolved m-FLIM data [9,10].The assumption that the number of components is known a priori is a common approach [5], however this estimation is complicated by the presence of outliers and experimental noise.
From the remote sensing literature, three distinct approaches for the estimation of the number of components can be distinguished: methods based on information theoretic criteria [14,15], eigenvalue based algorithms [11][12][13], and more recently geometry based strategies [16][17][18].Methods based on information theoretic criteria usually select a model, or families of probability density functions, that best fits the data according to a cost function.Some of these algorithms were first employed in passive sensor array processing [14,15].Other approaches are based on the eigenvalue decomposition of the observed data.These techniques also reduces the dimensionality of the problem and the computational cost.The method by Harsanyi et al. (1993) [11] introduced the use of Neyman-Pearson (NP) detection theory to make a decision based on a fixed probability of false alarm.The virtual dimensionality (VD) [12] method is based on the difference between the eigenvalues of the correlation and covariance matrix.This technique detects a noise component when the difference between the eigenvalues have similar magnitudes, and a signal source when the difference is significant.Another popular methodology is hyperspectral signal identification by minimum error (HySime) [13], which estimates the subspace of the eigenvectors that best contains the data according to a least squares criteria.
Geometry-based approaches are a recent trend, which is providing important advances in remote sensing.Hyperspectral intrinsic dimensionality estimation with nearest-neighbour distance ratios (HIDENN) [16] measures two geometrical properties from the data to estimate the fractal dimension of a manifold that contains the data.The intrinsic dimensionality of this manifold is then employed to estimate the effective dimensionality of the data.Outlier detection method (ODM) [17] characterizes data as outliers of noise in a principal component space.The noise components are spherically symmetric, while data signals are expected to present higher radius, in a hyper-ellipsoid fashion.These differences are measured by using an outlier detector based on inter-quartile ranges.The methods proposed by Ambikapathi et al. (2013) [18] are based on convex geometry.The constraints imposed on the abundances in a linear mixture model provide information to decide if the measured data should lie inside a convex hull (Gene-CH).In addition, if the end-members are linear independent they also lie inside an affine hull (Gene-AH).In both cases, the end-members are the vertex of the hulls.An NP detector is employed to validate the end-members.
It should be noted that the problem of estimating the number of components in the context of m-FLIM data has never been attacked before.Hence the technical contribution of this work relies on a new method for estimating the number of components in time-domain m-FLIM data by proposing two simultaneous NP tests to evaluate temporal and spatial consistency among the estimated end-members.These two tests are required mainly since the end-members are highly overlapped in m-FLIM data, compared with remote sensing applications.The first NP test evaluates temporal consistency by the linear independence of the resulting end-members (the estimated fluorescence time-decays), without using dimensionality reductions as in [18].The second NP detector checks the spatial consistency by the abundances or fractional contributions of the end-members in the data sample.Only those end-members with significant abundance levels are considered as original components; otherwise, they could be produced by experimental noise or outliers.By employing our previous contributions in [9,10], the new method implements a complete solution for the blind linear unmixing problem for m-FLIM data.That is, without prior information, our proposal thus estimates simultaneously the number of end-members, their time-domain fluorescence profiles and their relative abundances.Synthetic data was employed to measure quantitatively the performance of the algorithm.In addition, the algorithm was further validated using experimental m-FLIM data measured from in vivo samples of hamster pouches with induced oral carcinoma, and ex-vivo human atherosclerotic plaques.
The paper is organized as follows: Section 2 describes the linear unmixing model employed in the decomposition of m-FLIM data.In Section 3, the proposed NP detectors for linear independence and spatial coherence are presented in detail.Also, the BEAE is briefly introduced, and the final algorithm for estimating the number of components, the end-members and their abundances is detailed at the end of this section.The experiments and their results for the synthetic and experimental m-FLIM datasets are shown in Sections 4 and 5, respectively.Section 6 introduces a discussion on the results, and the final remarks are presented in Section 7.

Linear unmixing of m-FLIM data
The linear mixture model [19][20][21] assumes that every measurement is a linear combination of the end-members.In the case of time-domain m-FLIM data [9], the measurements and the end-members are fluorescence intensity decays measured in different spectral bands simultaneously, which can be concatenated as shown in Fig. 1(a).The fluorescence decays are a combination of the endmembers p j ∈ R L for all j ∈ [1, N] with their respective concentrations α k, j , which are usually referred as abundances in remote sensing literature.This linear mixture is given by The mixture model can be rewritten in matrix form to include all the K measurements available.For this purpose, the end-members matrix P N = [p 1 ,..., p N ] ∈ R L×N gathers the time-domain fluorescence profiles of the N components, and the abundances are grouped in vectors α k = [α k,1 , . . . ,α k,N ] at each spatial location k.Using these definitions, the mixture model can be written as where Y ∈ R L×K contains the available m-FLIM measurements, and the abundance matrix This model considers the following two constraints for the abundances at each spatial location k [5], since these parameters represent the proportional concentration of the end-members: where 1 N represents a vector with unitary entries of dimension N, and the inequality ≥ is interpreted component-wise.Furthermore, the linear mixture model in (1) subject to the constrains in ( 3) and ( 4) can have a geometrical interpretation: the set of mixture samples Y = {y 1 ,...,y K } should be a subset of a convex hull [22], where the end-members p n represent its vertices.The corresponding convex hull can be written as Moreover, our formulation requires that the end-members {p 1 , . . ., p N } are linearly independent [23], and as a consequence Ω N p is a simplex of size N − 1.A notable difference between linear unmixing for mFLIM data and remote sensing information is the incorporation of constraints in the estimation of end-members [10].The fluorescence intensity time decays are strictly positive, and they are usually normalized in order to minimize intensity artifact problems [19,24].Therefore, in our application, there are two additional constraints for the end-members with respect to the linear unmixing model p j ≥ 0. (7)

Estimating the number of components
The accurate estimation of end-members P N and their abundances A N rely on a correct knowledge of the number of components N in the m-FLIM data, which by itself represents a difficult task due to the presence of noise and outliers [5].Outliers are generated by measurement errors or components with very low occurrence in the sample.The proposed algorithm is based on our previous methodology for blind unmixing or BEAE [10] that requires a linear independence property of the end-members.BEAE relies on alternating least squares (ALS) [25] and restricted quadratic optimization to compute the end-members and their abundances in the sample [9].For this new challenge, we propose an iterative methodology to increase gradually the estimation of the number of components, if at each step the resulting end-members satisfy a joint coherence test to validate their independence and their spatial interpretation.Once all candidates are good estimations of the fluorescent components present in the sample, the set of end-members will satisfy the following criteria: Linear Independence: The convex hull ( 5) is a simplex and the set of end-members should be affinely independent [18,23].
Spatial Coherence: Each end-member has a significant relative abundance, since it is contributing to the sample fluorescence signal.
These criteria are evaluated by using NP tests [26], as detailed in the following subsections.

Linear independence criterion
In the BEAE methodology [10], we perform the linear unmixing of the m-FLIM samples while estimating a set of candidate end-members: where N > 0 is the estimated number of components.After an iteration, we add one more candidate to the set (8) and increase our estimation of the number of end-members to N = N +1.As long as the number of candidates N is lower or equal to the real number of components N, i.e.N ≤ N, the set of end-members (8) will be affinely independent and Ω N p will be a simplex.To detect if the j-th candidate end-member p j is a linear combination of the rest, we propose the following approach based on the methodology by Ambikapathi et al. (2013) [18].For each candidate p j ∈ P N , its projection into the closed-subspace spanned by {p 1 , . . ., p j−1 , p j+1 ,..., p N } is calculated as p j = P j β j (9) where The optimal projection is estimated by a quadratic criterion as subject to 1 N−1 β = 1 and β ≥ 0, where • 2 denotes the Euclidean norm.A direct solution for ( 10) can be obtained by using the methodology of Gutierrez-Navarro et al. (2014) [9].In our formulation, we consider that the estimation of the end-members is not perfect, and there is some additive uncertainty which is represented by a zero-mean Gaussian vector w j ∈ R L with covariance Σ > 0, i.e.
where p o j denotes the true j-th end-member.In order to detect fake end-members, we define the j-th error vector where 1) .Note that the new uncertainty term v j is still a zero-mean Gaussian vector, but with covariance κΣ > 0, where Therefore, in our evaluation, there are only two possible cases: • The candidate p j is a linear combination {p 1 ,..., p j−1 , p j+1 ,..., p N }, which implies that μ j = 0 since the j-th end-member p o j can be represented exactly by P o j β j , Ω N p is not a simplex and therefore N > N, i.e. the error vector e j is only affected by the additive uncertainty v j : where N (a, B) denotes a Gaussian distribution with mean a and covariance B.
• The set in ( 8) is linearly independent and p j is a valid end-member, and as a consequence μ j = 0, i.e. the error vector e j is affected by a constant term μ j plus the additive uncertainty v j : The noise covariance matrix of the extracted end-members Σ is estimated by using the methodology from Bioucas-Dias and Nascimento (2008) [13].Hence, from our previous analysis, we are facing a binary hypothesis test, where the observation (error vector) in ( 11) presents a Gaussian distribution.Therefore, to simplify the evaluation, a quadratic decision variable r I j is constructed by the weighted Euclidean norm of the error vector: where the inverse of the resulting covariance in the error (κΣ) −1 is employed to normalize r I j .Since the error vector e j is Gaussian, then the new variable r j will have a χ 2 -distribution.In this way, a binary hypothesis testing problem [26] is defined as: where L − 1 are the degrees of freedom in the distribution associated to the length of the error vector e j , χ 2 (L − 1) denotes a central χ 2 distribution, and χ 2 L − 1, μ j the non-central one with noncentrality parameter μ j .The expression in (18) cannot be evaluated since μ j is unknown.However, we can employ the NP approach [26] to overcome this limitation.This method distinguishes between the null H 0 and alternative H 1 hypotheses by fixing the false alarm probability for this test to a predefined value P I FA , i.e.
where η j is a constant that satisfies (19), and f χ (x, L − 1) is the χ 2 -probability distribution with respect to H 0 .Thus the decision can be made based only on a criterion that employs the distribution of the null hypothesis H 0 for j-th end-member, as and we can define the decision rule for our problem as Since r I j is the sum of L squared independent random variables with finite mean and unit variance, we can further simplify our decision rule.By the central limit theorem, for a large value of L, the variable r I j is close to a normal distribution.In the case of m-FLIM data, L is usually several hundreds.Hence, we propose to simplify the calculations by using a linear mapping of r j to normalize and obtain a new decision variable: Our decision problem for the linear independence criteria is finally stated as

Spatial coherence criterion
In the BEAE methodology [10], if the number of components is higher than the real value, the extra end-members will be a linear combination of the real ones.However, low signal-to-noise ratios and the presence of outliers may lead to evaluate an end-member that could deceive the linear independence criteria.On the other hand, end-members with high abundance throughout the m-FLIM sample are responsible for most of the fluorescence intensity.Even components with a localized abundance could present high concentration levels, which could be the case for some dyes or molecules of interest.Hence, end-members that represent real components in the sample data should not only be linearly independent but also have a significant spatial contribution, which is quantified by their abundances.Under this hypothesis, we need to discard end-members with a low relative abundance throughout the m-FLIM data, since they will not have spatial coherence.The abundance information of the end-members at each spatial location of the sample is provided by the BEAE algorithm in matrix A N .
To evaluate the spatial coherence, for j-th end-member, its abundance is extracted from matrix A N by the j-th row, and an histogram h j ∈ R B is constructed with B bins or levels, defined by the vector b ∈ R B .Therefore, the i-th element h j,i in h j represents the number of occurrences where the abundance belonged to the range (b i−1 , b i ] with b −1 = 0 for j-th end-member.As a result, each histogram satisfies 1 B h j = K for all j ∈ [1, N], and we define an spatial coherence descriptor s j = [s j,1 , . . ., s j,B ] with binary components: where θ is the minimum number of occurrences for an abundance bin to be considered significant, and h j = [h j,1 , . . . ,h j,B ] .In this work, we fixed the threshold θ = 0.005K, i.e. 0.5% of the total spatial sample.End-members with significant abundance values throughout the m-FLIM sample will present mostly spatial coherence descriptors s j in (24) with all unitary entries.On the contrary, the descriptor obtained from a bad candidate will be zero in entries related to low abundances.Therefore, to evaluate the presence of bad end-member candidates, which present only low-valued abundances, we can use the information provided by the spatial coherence descriptor s j .With this information at hand, we propose a new decision variable r II j ∈ R B to evaluate the spatial coherence of j-th end-member As a result, r II j is an indicator of the mean value of the abundances with significant contributions throughout the m-FLIM data sample.Once again, a binary hypothesis test [26] will be faced: The decision variable r II j has a low value since it represents an end-member with a lowvalued abundance, with the following observation model: where U (0, ω 1 ) represents a uniform distribution in the interval [0, ω 1 ], and ω 1 ∈ (0, 1].
H II 1 : If the j-th end-member does have a significant contribution in the m-FLIM sample then where ω 2 ∈ (0, 1] and ω 2 > ω 1 . Since it is not possible to characterize the decision variable r II j for H II 1 defined by the upper bound ω 2 , as in the linear independence criterion, we propose another NP test based on the spatial coherence descriptor: (28) where f [0,ω 1 ] U (•) denotes the uniform probability distribution in the interval [0, ω 1 ], and P II FA the probability of false alarm for the spatial coherence test.As a result, the key parameter in the spatial coherence criterion is ω 1 .

Blind end-member and abundance extraction algorithm
Our proposal is based on an ALS procedure [10] which minimizes the quadratic cost function min subject to constraints (4), (3), ( 6) and (7), where • F denotes the Frobenius norm for matrices, and ξ > 0 is a weight parameter related to the regularization term.The complete algorithm employed in this paper is described in Algorithm 1.
Algorithm 1: Estimation of the Number of Components and Linear Unmixing of m-FLIM Data.
input : An m-FLIM dataset Y of size L × K. Parameters: minimum abundance θ , upper bound on the distribution of the spatial descriptor ω 1 , the probabilities of false alarm P I FA and P II FA , the maximum number of components N MAX , the number of bins B and a tolerance parameter ε. output: Estimation of the number of components N, the end-members P N and the abundance matrix A N . 1 Preprocess: Set negative elements in Y to 0, and normalize every measurement to satisfy

Synthetic experiments
The performance of Algorithm 1 was first evaluated by using computer simulations.Synthetic m-FLIM data was created by considering the five synthetic end-members P 5 depicted in Fig. 1(a).The end-members were generated using the stretched exponential model [27], and designed to be linearly independent.In multispectral time-domain fluorescence techniques [28,29], the statistical noise in the photodetector follows a Poisson distribution.However, when there is a large photon count, this stochastic term approaches a Normal distribution with zero mean [30], whose standard deviation is proportional to the square root of the fluorescence intensity.As a result, the noise components have a variable variance that depends on the actual fluorescence intensity.Hence, the synthetic data was produced according to the following model: where E ∈ R L×K is a noise matrix defined as E = [e 1 ,...,e K ].The noise components e k = [e k,1 , . . ., e k,L ] ∈ R L at position k are proportional to the square root of the intensity at that location and time instant, e k,i = δ k,i √ y k,i for all k ∈ [1, K] and i ∈ [1, L], where δ k,i ∼ N 0, σ 2 .In this way, the standard deviation of the noise component at k location and i time instant will be equal to σ √ y k,i .Since the equivalent signal-to-noise ratios will be time-varying at each spatial sample, the values of σ were selected to satisfy a Peak Signal-to-Noise (PSNR) ratio of 15, 20 and 25 dB according to for all spatial locations k in the sample.The matrices of synthetic abundances A 5 were selected to assign purity levels [18] ρ equal to 1.0, 0.9 and 0.8, see Fig. 1(b).The purity level ρ is a measure of the presence of pure samples in the data.The value of ρ for the j-th end-member is defined by its spatial abundance as A value ρ = 1 indicates the presence of pure samples of the j-th component in the synthetic data.In fact, a ρ value lower than 1 makes the linear unmixing decomposition more difficult [5,22].We analysed the influence of P FA and the upper bound on the distribution of the descriptor ω 1 of the spatial coherence test in the NP evaluations.The probabilities of false alarm were set to the same value for both tests P I FA = P II FA = P FA , and the number of bins B was fixed to 30 bins.This value of B balanced a compromise between accuracy and complexity in our evaluations.The size of each m-FLIM dataset was 75 × 75 × 510, where the first two dimensions correspond to the spatial plane of the sample (the same dimension of the abundance maps in Fig. 1(b), and the third dimension is time, as shown in Fig. 1(a).A total of 50 synthetic data sets were employed for each PSNR level.

In-vivo experimentation with hamster pouches
We employed in-vivo m-FLIM data taken from the cheek pouches of a golden Syrian hamster.The data was obtained from Jabbour et al. (2013) [29], where the protocol employed to treat the hamster and the tissue extraction are detailed.One dataset corresponds to the left pouch, which was left untreated for control purposes.The right pouch was treated with 7, 12dimethylbenz[α]anthracene (DMBA) to induce oral carcinogenesis.Both tissue samples have an approximate area of 16 × 16 cm 2 , and they were imaged using a 400 × 400 × 1200 m-FLIM dataset.The first two dimensions correspond to the X-Y plane of the tissue sample.The third dimension corresponds to time-domain fluorescence intensity decays measured with a temporal resolution of 160 ps in three wavelength bands: 390 nm, 452 nm, and >500 nm.
The m-FLIM samples were processed together to estimate their common end-members and abundances.Measurements over the edges of both samples were discarded to reduce border effects.For the estimation of the end-members and due to computational limitations, the size of the initial m-FLIM dataset was reduced in the X-Y plane to 300 × 300 × 1200, but afterwards the complete m-FLIM dataset was processed to estimate the abundances in the original X-Y plane 400 × 400 measurements.Since in the experimental data a high PSNR is considered, the parameters employed in the NP tests were ω 1 = 0.2, B = 30 bins and the threshold for the minimum number of significant measurements θ was set to 225 positions.

Ex-vivo experimentation with atherosclerotic plaques
The proposed methodology was also tested with m-FLIM data from ex-vivo human coronary arteries.The samples were obtained from Park et al. (2011) [28].The temporal resolution of the measurements was 250 ps, and they were recorded in the wavelength bands: 390 ± 20 nm, 452 ± 22.5 nm and 550 ± 20 nm.The datasets were classified by an expert [28] using histopathology analysis.Two samples were labeled as Low-Collagen/Lipids (LCL), two datasets were identified as having High-Collagen (HC), and another two as High-Lipids (HL), while the last three samples were classified as Mixtures (Mix), since none component was dominant.Nine individual datasets with a field of view 2 × 2 mm 2 (60 × 60 × 510) were analysed jointly with Algorithm 1 to estimate the number of components and estimate their quantitative descriptions.The same configuration of parameters employed in section 4.2 were used for these experiments.

Oral carcinogenesis in hamster pouches
Three end-members were detected ( N = 3) when solving the hamster m-FLIM datasets simultaneously.By recalling that j-th end-members p j , represents the concatenated response at three wavelengths (p λ 1 j , p λ 2 j , p j ) ] , we can evaluate the accuracy and interpretation of the estimated end-members by the resulting average lifetime and normalized intensity at each frequency band [31].For this purpose, we compute in discrete-time the average lifetime for each j-th end-member at wavelenght band λ i : where pλ i j represents the end-member p λ i j after performing the time-deconvolution [31] with the instrument response, and t is the time vector associated to the periodic sampling during the experiment.The normalized intensity for each wavelength band λ i and j-th end-member is expressed as The average lifetimes and normalized intensities are shown in Table 2. Based on this information, the end-members were identified as Collagen, Porphyrin and NADH/FAD.The endmembers calculated are depicted in Fig. 2. The abundance maps obtained from 400 × 400 measurements of the m-FLIM dataset are shown in Fig. 3.

Component j
Wavelength λ Lifetime Normalized Intensity 390 nm 452 nm >500 nm 390 nm 452 nm >500 nm τ λ 1 The end-members were identified as Collagen, Porphyrin and NADH/FAD according to the lifetimes calculated in Table 2 .

Atherosclerotic plaques
M-FLIM samples from ex-vivo atherosclerotic plaques were analysed using the methodology in Algorithm 1.Three wavelengths were also considered in the m-FLIM measurements: 390 nm, 452 nm, and 550 nm.Once more, three were detected ( N = 3) and their deconvoluted lifetimes were also calculated to eq. (34).The average lifetimes are in Table where now the end-members were as Collagen, Elastin and Low Density Lipoproteins (LDL).The end-members and their abundances are shown in Figs. 4 and 5, respectively.

Discussion
Blind linear unmixing of m-FLIM data allows to obtain a quantitative description of the biochemical composition of the biological samples, which is easy to interpret and validate.Usually, the interpretation and validation is difficult without prior knowledge of the number of components in the sample, and their characteristic signatures.This task is even more difficult when unreliable or partial information is available, as is usually the case when working with biomed- Table 3.Average Lifetime and Normalized Intensities of the end-members extracted from m-FLIM datasets of atherosclerotic plaques.data.In the literature, there are some approaches targeted for remote sensing applications [5], however, there are intrinsic differences.First, satellite imagery usually contains a high number of end-members (> 30), while the number of endogenous fluorophores present in an m-FLIM sample is usually low.This property could be seen as an advantage, but in m-FLIM data, the end-members are highly overlapped which complicates the linear unmixing problem.Furthermore, the BEAE method [10] takes advantage of the non-negativity nature of fluorescence, which is translated into a normalization of the end-members that helps to constraint the solution space [23,32], and reduce the signal variability in the data [10].

h h h h h h h h h h h h h h h
In the evaluation process, experimentation with synthetic data is the only way to quantitatively measure the performance of the proposal, since it is difficult to estimate the exact and localized abundance of biological samples without destroying them.These experiments are also useful to determine the best parameter configuration in the proposed methodology.According to the results obtained in Table 1, the proposed algorithm presented good performance in estimating the number of components in the synthetic m-FLIM data at PSNR levels of at least 20 dB, independently of P FA and ρ.For PSNR = 20 dB, ω 1 = 0.1 performed the best, while for PSNR = 25 dB either ω 1 = 0.1 or ω 1 = 0.2 performed similarly well.
The processing of the m-FLIM data of hamster pouches determined the presence of three components.Using their estimated lifetime values and relative normalized intensities for each measured emission band, we could identify the first one as Collagen, whose fluorescence has a peak emission at ≈ 400 nm and lifetime values in the order of 4-6 ns.The second end-member was identified as Porphyrin, whose fluorescence has not emission below 500 nm and lifetime values longer than 7-8 ns [29].The third end-member presented lifetimes that did not match any expected end-member found in hamsters oral mucosa.The lifetime of this third end-member in the 390 nm wavelength band alone matched the description of NADH, however, the lifetime is too long in the wavelengths bands 450 nm and > 500 nm.These lifetimes actually seem to match the description of FAD, which is another of the expected primary endogenous fluorophores [29].This is the reason why the third end-member was labelled as a mixture of NADH and FAD.We believe that m-FLIM data presented no pure samples of NADH and FAD, which is why they could not be separated.The third end-member profile was characterized by a peak emission at the 450 nm channel, and lifetime values between ≈1.5-3.5 ns.These fluorescence characteristics can be associated to NADH, whose fluorescence emission also peaks at ≈450 nm with lifetime values between ≈0.5-3 ns, depending whether it is in a free or bound state.On the other hand, FAD has a peak emission ≈510 nm and similar lifetime values.Since the third end-member profile also has significant contribution in the 500 nm band, it could be possible that this end-member reflects the contribution of both NADH and FAD.The fact that these two molecules are usually co-localized within the cell cytoplasm further support this hypothesis.Thus, the third end-member was attributed to a combined contribution of both NADH and FAD.
The healthy pouch, depicted in Fig. 3(a) presented a high relative concentration of Collagen, with some NADH/FAD contribution in the lower left section, and practically no Porphyrin through the imaged area.On the contrary, the right pouch presents a tumour which is visible in the lower right section of the abundance maps, as shown on Fig. 3(b).The resulting abundances reveal that the tissue surrounding the tumour contains high levels of Collagen and lower levels of Porphyrin, in a similar fashion to the healthy tissue in Fig. 3(a).The tumour area presents high abundance of Porphyrin, reaching relative concentration values of 100% in some areas.The relative concentration of NADH/FAD was also high in the tumor region, and has relative high presence throughout the entire imaged are of the DMBA treated cheek pouch.The quantitative results obtained in the hamster m-FLIM samples matches the qualitative description provided by Jabbour et al. (2013) [29].The detection of Porphyrin and NADH/FAD could play

Fig. 1 .
Fig. 1.The 5 synthetic end-members (a) and their abundance maps (b), with purity level ρ = 1, employed in the synthetic evaluation.The end-members of 510 samples represent the concatenated fluorescence intensity time-decays at three different wavelength bands.The x and y positions (75 × 75) in (b) correspond to spatial coordinates in the sample, while the color depicts the quantitative abundance values.

Fig. 2 .
Fig.2.End-members extracted from the m-FLIM datasets of hamster pouches.The second end-member presents practically no emission fluorescence signal on the first two bands.The end-members were identified as Collagen, Porphyrin and NADH/FAD according to the lifetimes calculated in Table2

#Fig. 3 .
Fig. 3. Abundance maps (400 × 400) corresponding to the end-members in Fig. 2 detected in the m-FLIM datasets of hamster pouches.The abundance maps shown in row (a) correspond to the quantitative description obtained from the healthy hamster pouch.Meanwhile, the diseased tissue produced the abundance maps shown in row (b).

Fig. 4 .
Fig. 4. End-members extracted from the m-FLIM datasets of atherosclerotic plaques, which were identified as Collagen, Elastin and LDL.

Fig. 5 .
Fig. 5. Abundance maps (60 × 60) for the m-FLIM datasets of atherosclerotic plaques, where each row correspond to a quantitative description.The m-FLIM datasets were identified through histopathology analysis as LCL, HC, HL and Mix, as indicated in the picture.

N t+1 25 END 26 end 27 end
the initial candidates as P 2 = [p 1 , p 2 ] where [10], there is only one component present in the sample.9Sett=0andP N t = P N and calculate A N t using BEAE in[10]13Set t = t + 1 14 until L(P N t+1 , A N t+1 ) − L(P N t , A N t ) < ε;15 A N = A N t+1 and P N = P N t+1 16 if N > 2 then // We have at least 3 candidates 17 for j ∈ [1, N] do // Evaluate Candidates 18 For j-th end-member, evaluate the decision variables r I j and r II j in (16) and (25), respectively 19 Test the linear independence criterion in (23) 20 Test the spatial coherence criterion in (28) 22 Set N = N − 1 23 Perform the ALS procedure one last time (Lines 8 to 13) 24 Set A N = A N t+1 and P N = P // No fake candidates detected, test more end-members 28 Choose the next candidate as the measurement with the largest estimation error, i.e. the measurement with the poorest estimation by the N end-members p N+1 y ko where k o = arg max 29 Set P N+1 = [P N p N+1 ] and N = N + 1 #208447 -$15.00USD Received 17 Mar 2014; revised 2 May 2014; accepted 2 May 2014; published 13 May 2014 (C) 2014 OSA

Table 2 .
Average lifetimes and normalized intensities of the end-members extracted from the m-FLIM datasets of hamster pouches.The 2nd end-member presented no response in the first and second wavelength band.
The method proposed in this paper performs the linear decomposition of time-domain m-FLIM data, while it evaluates the accuracy of the solution by the linear independence and spatial coherence NP tests.As a result, our proposal takes a decision based on the quality of the end-members and their abundances obtained in an iterative fashion.The iteration stops once the current solution is not able to satisfy the linear independence criterion among the end-members, and the spatial coherence related to their abundances.To the best of the authors' knowledge, this is the first completely blind algorithm tailored for linear unmixing of time-domain m-FLIM #208447 -$15.00USD Received 17 Mar 2014; revised 2 May 2014; accepted 2 May 2014; published 13 May 2014 (C) 2014 OSA