 Research
 Open Access
 Published:
Twohierarchical nonnegative matrix factorization distinguishing the fluorescent targets from autofluorescence for fluorescence imaging
BioMedical Engineering OnLine volume 14, Article number: 116 (2015)
Abstract
Background
Nonnegative matrix factorization (NMF) has been used in blind fluorescence unmixing for multispectral invivo fluorescence imaging, which decomposes a mixed source data into a set of constituent fluorescence spectra and corresponding concentrations. However, most classical NMF algorithms have ill convergence problems and they always fail to unmix multiple fluorescent targets from background autofluorescence for the sparse acquisition of multispectral fluorescence imaging, which introduces incomplete measurements and severe discontinuities in multispectral fluorescence emissions across the multiple spectral bands.
Methods
Observing the spatial distinction between the diffusive autofluorescence and the sparse fluorescent targets, we propose to separate the mixed sparse multispectral data into equality constrained twohierarchical updating within NMF framework by dividing the concentration matrix of entire endmembers into two hierarchies: the fluorescence targets and the background autofluorescence. Specifically, when updating concentrations of multiple fluorescent targets in the twohierarchical NMF, we assume that the concentration of autofluorescence is fixed and known, and vice versa. Furthermore, a sparsity constraint is imposed on the concentration matrix components of fluorescence targets only.
Results
Synthetic data sets, in vivo fluorescence imaging data are employed to demonstrate and validate the performance of our approach. The proposed algorithm can achieve more satisfying results of spectral unmixing and autofluorescence removal compared to other stateoftheart methods, especially for the sparse multispectral fluorescence imaging.
Conclusions
The proposed algorithm can successfully tackle the sparse acquisition and illposed problems in the NMFbased fluorescence unmixing through equality constraint along with partial sparsity constraint during twohierarchical NMF optimization, at which fixing sparsity constrained target fluorescence can make the update of autofluorescence as accurate as possible and vice versa.
Background
Due to Shimomura, Chalfie and Tsien’s shared discovery and development of the green fluorescent protein [1], noninvasive invivo fluorescence imaging has become a powerful technique at labeling and tracking cells in tissues for biological and medical sciences, such as drug discovery and disease diagnosis [2]. However, the overlap of multiple fluorescence spectra introduces the crosstalk between multiple fluorescent targets in multispectral fluorescence imaging. Especially, when excited in the visible to nearinfrared part of the electromagnetic spectrum, some natural fluorescent molecules in tissues and food (such as collagens, lipofuscin, and melanin) can also produce fluorescence photons called autofluorescence (AF) [2–5], which spreads everywhere in tissues with typically containing other nonsparse instrumentbased background noise (such as dark current, residual bulk image in CCD, and leakage light from exciting filters). Therefore, the strong spatial and spectral overlap between different target fluorophores and background AF makes it very difficult to unmix multiple fluorescent targets, complicating the analysis of multiple fluorophores in multispectral fluorescence imaging.
To solve the overlapping problem in fluorescence imaging, spectral unmixing (SUM) is used to decompose the mixed spectra data \(\mathbf {D}\) into a product of pure spectral signatures \(\mathbf {S}\), i.e., endmembers, and corresponding concentrations \(\mathbf {C}\) (or abundances, mixing weights), representing the contribution to the observed fluorescence radiance from the corresponding endmember. Before we dig into some fluorescence unmixing methods, it will help us to have a general understanding of the following observations in fluorescence imaging: first, due to invivo fluorescence imaging measuring the diffuse radiance from the surface of scattering animal tissue after the fluorescent light transport in tissues, a fluorophore inside animal tissue might make a significant contribution to a cloud of neighboring pixels instead of only the pixel geometrically associated with it, such that the quantity \(\mathbf {C}\) for multiple fluorophores represents the signal distribution of these multiple fluorophores instead of the concentrations of fluorophores in the pixel [3]; second, the emitted fluorescent light travel in diffusive tissues can result in the decreasing intensity and spectral distortions. However, the spectrum distortion of signatures \(\mathbf {S}\) is insignificant in front of intensity loss of deep embedded fluorophores [6]. Therefore, the signatures of multiple diffusive fluorescent sources are relatively stable in contributing to the mixed spectra data in fluorescence imaging; third, as the fluorescence signal attenuates exponentially with the distance light travels, a background AF signal remains constant and becomes a serious limiting factor. The analysis of a fluorescence signal polluted by the background AF may lead to a wrong localization of the target fluorophore. In the light of the above mentioned observations, the spectral unmixing is therefore precisely referred to as a blind source separation approach, which is a necessary preprocessing step to remove unwanted AF and unmix different spectra of interest. By returning separated pure fluorescence contribution data, the SUM method not only gives contrastenhanced [6] signal distribution of each target fluorophore in 2D planar imaging to enable comparing fluorescence sources at similar imaging conditions and comparable depth, but also presents the possibility to perform accurate fluorescence tomographic reconstructions confirming the 3D position of marked samples.
According to different endmember estimation methods, SUM has been classified as either supervised [3, 7] or unsupervised. In practice, a userselected region containing a single fluorescent target or a spectral library for multiple known fluorescent targets present in the sample is chosen to provide the endmember information when implementing supervised SUM. However, it is difficult to select an appropriate region solely containing single fluorescence or to make extensive calibration in the spectral library acquisition ensuring the accuracy of the unmixing algorithm. Therefore, various unsupervised SUM methods [8, 9] (also called blind source separation) have been developed to estimate spectra and concentrations simultaneously without exact priori knowledge about constituent spectra.
To blindly apply spectral unmixing for fluorescence imaging, nonnegative matrix factorization (NMF) [10, 11] has been used to find the matrices \(\mathbf {S}\) and \(\mathbf {C}\) by minimizing the cost function F representable as a difference between \(\mathbf {D}\) and \(\mathbf {CS}\) in fluorescence spectroscopy [4, 12–15]. However, due to the cost function of NMF being not convex in the search spaces of \(\mathbf {S}\) and \(\mathbf {C}\) simultaneously, starting from different initial search points results in different values for the elements of \(\mathbf {S}\) and \(\mathbf {C}\). Therefore, the NMF suffers form the ill convergence problem and the final solution is easily caught in the local minimum of the cost function. Especially for the sparse multispectral fluorescence imaging data that are acquired in a relatively small number (3–10) of spectral bands, the ill convergence problem is easily exacerbated due to the sample deficiency and severe discontinuity of multispectral fluorescence emissions, though this sparse imaging can be fast and costeffective in clinical applications. In multispectral fluorescence imaging, the multipletarget fluorophores with their unimodal spectrum emission peaks are sparse sources in both the spatial and the spectral domains. Such source sparsity is widely used as a regularized constraint in fluorescence unmixing studies [4, 12–14, 16], among which a few studies were interested in the sparsity of matrix \(\mathbf {S}\). However, unlike the hyperspectral imaging [17, 18], the sparse acquisition of multispectral fluorescence imaging is inevitable to cause spectrum leakage and fence effect in the spectrum of fluorescence target such that the rows of \(\mathbf {S}\) are not true expressions of spectrally sparse fluorescence emissions. Furthermore, Montcuquet et al. [4] reported that, due to the currently used fluorescent markers traveling in blood and lymphatic canals and always spreading socalled nonspecific signal everywhere in tissues, the spatial sparsity constraints alone may seem inappropriate for the unmixing methods to detect the fluorescent source signal from the nonspecific signal. Moreoverand to the best of our knowledgeNMF methods fail in discriminating the spatial sparsity of the multipletarget fluorescence signal and the spatial nonsparsity of the AF (including the nonspecific signal).
Considering the abovementioned factors, we propose a twohierarchical NMF (thNMF) procedure to distinguish the updating of matrix components for the multiple fluorescent targets from the background AF. Our proposed twohierarchical NMF is first inspired by the equality constraint in alternating least squares procedures [19], where the known columns of \(\mathbf {C}\) are fixed during the update process and the unknown columns of \(\mathbf {C}\) are updated with the subcost function which has already removed the part of known columns of \(\mathbf {C}\) and the corresponding rows of \(\mathbf {S}\). To some degree, equality constrained least squares algorithm can be regarded as separating one known hierarchy from one unknown hierarchy during the update process although the known hierarchy is fixed within hierarchical alternating least squares (HALS) algorithm [20–22].
By exploiting the difference between the spatial distributions of target fluorophores and background AF within the equality constrained HALS framework, we propose thNMF with the following two main contributions: first, rather than updating the whole matrix \(\mathbf {C}\) (and corresponding \(\mathbf {S}\)) in minimizing the cost function F, the thNMF treats the diffusive AF separately from the sparse fluorescent targets and divides the whole minimization process into two different subprocesses, which fix the known column of matrix \(\mathbf {C}\) for background AF (resp. target fluorophores) during the updating unknown columns of \(\mathbf {C}\) for target fluorophores (resp. background AF) with the subcost function having removed the part of known columns of \(\mathbf {C}\) and the corresponding rows of \(\mathbf {S}\), whereby we can prevent thMNF from mixing these two different components in NMF decomposition; second, the twohierarchical scheme can potentially facilitate imposing different appropriate constraints on the target and background components in the two subprocesses to obtain desired unmixing results. In this paper, according to the spatial sparsity difference between the target fluorophores and background AF, we solely impose the partial sparsity constraint on the fluorescent targets to accurately separate the sparse multitarget fluorophores and nonsparse background AF. The remainder of this paper is organized as follows. Section II describes the materials and methods of proposed thNMF algorithm. Section III provides experimental results on synthetic and invivo fluorescence imaging data. A discussion and conclusions are given in section IV.
Methods
Twohierarchical NMF
Most existing spectral unmixing approaches are based on linear mixture model, which can be expressed as
where matrix \(\mathbf {D}\in \mathbb {R}_{+}^{N\times L}\) describes multispectral fluorescence imaging data where each column corresponds to a single fluorescent image of a special spectral band. N is the number of total pixels in a single image and L is the number of acquisition spectral bands. \(\mathbf {C} = [{\mathbf {c}_1},{\mathbf {c}_2}, \ldots ,{\mathbf {c}_k}] \in {\mathbb {R}_+^{N \times K}}\) is called concentration matrix where each column \({\mathbf {c}_i}\) represents the concentration distribution of a special endmember and \(\mathbf {S} = {[{\mathbf {s}_1},{\mathbf {s}_2}, \ldots ,{\mathbf {s}_k}]^T} \in {\mathbb {R}_+^{K \times L}}\) is the endmember spectrum matrix where \(\mathbf {s}_i\) is a vector corresponding to the emission spectrum of the ith endmember. K denotes the number of endmembers. \(\mathbf {E}\) is the related noise matrix.
The aim of NMF methods is to find the nonnegative best solution couple (\(\mathbf {C,S}\)) whose product \(\mathbf {CS}\) best approaches the original data matrix \(\mathbf {D}\). To solve the NMF problem, an alternating nonnegative least squares (ANLS) algorithm [22, 23] regards it as a convex optimization subproblem to find the optimal factor \(\mathbf {C}\) corresponding to a fixed factor \(\mathbf {S}\) and vice versa. When partitioning variables into one single variable at a time in each convex subproblem in the framework of ANLS, we can get a simple univariate quadratic problem, such as \({\mathbf {C}_{ik}} \leftarrow \text {argmin}{_{{\mathbf {C}_{ik}} \ge 0}}\left\ {\mathbf {DCS}} \right\ _2^2\), which admits a closedform solution. Moreover, since the optimal value for a given entry of \(\mathbf {C}\) does not depend on the other entries of the same column, one can optimize alternatively whole columns of \(\mathbf {C}\), with
where \({\mathbf {R}_k} = \mathbf {D}  \sum \nolimits _{i \ne k} {{\mathbf {C}_{:i}}{\mathbf {S}_{i:}}}\) is the kth residual matrix. This new method is referred to as HALS [20–22]. Under the framework of HALS algorithm, it seems convenient to add different constraints to each column of \(\mathbf {C}\) or each row of \(\mathbf {S}\). HALS algorithms with a sparsity constraint [25] and a constraint for maximum spatial dispersion [26] have been applied in the area of hyperspectral unmixing.
Actually, in the HALS computation framework, the multiple fluorescent targets and background AF are amenable to being treated as different hierarchies due to their having fundamentally different characteristics in the spatial and spectral distributions. Firstly, the fluorescent targets are locally accumulated in specific biological tissues and have sparse distributions in the spatial areas, while background AF propagates at all directions and diffuses uniformly over large spatial areas. Secondly, the spectrum of a single fluorescent target is narrow and only has one emission peak in spectral excitation and emission, while the background AF widely spreads at all visible light range and has a slowly varying spectral distribution. However, this spectrum distinction between fluorescence target and background AF is substantially weakened in the sparse multispectral fluorescence imaging, which introduces spectrum leakages between neighboring spectral bands and severe discontinuities of spectral emissions for both fluorescence targets and background AF. Therefore, it is reasonable for NMF to emphasize on the continuous characteristics of the spatial distributions for fluorescent targets and background AF and only introduce appropriate regularization constraint into the concentration matrix \(\mathbf {C}\) in this work.
As for exploiting the different spatial regularizations for the fluorescent targets and background AF, we divide the whole update of concentration matrix \(\mathbf {C}\) into two different but relevant subproblems for the fluorescent targets and background AF. Inspired by the equality constraints on variables during alternating least squares procedures [19] and the HALS, we propose twohierarchical NMF for the sparse multispectral fluorescence imaging. Assuming that the original cost function with an appropriate sparsity constraint \(J(\mathbf {C})\) is \(F\left( \mathbf {C},\mathbf {S} \right) =\left\ \mathbf {D}\mathbf {CS} \right\ _{2}^{2}+\lambda J\left( \mathbf {C} \right)\) , where \(\lambda \in {{\mathbb {R}}^{+}}\) is a scalar that weighs the contribution of the sparsity measure function \(J(\mathbf {C})\) of the matrix \(\mathbf {C}\), which is usually regarded as a regularization term. According to the ANLS algorithm, we can obtain the update rules
Considering the concentration distinction between fluorescent targets and background AF, we let \({\mathbf {C}=[\mathbf {C}_F\ \mathbf {C}_A]}\) and \({\mathbf {S} = [\begin{array}{l} {{\mathbf {S}_F}} \\ {{\mathbf {S}_A}} \\ \end{array}]}\). Then the problem of Eq. (5) can be divided into two subproblems:
where \({\mathbf {C}_F}\) and \({\mathbf {S}_F}\) denote the concentrations and fluorescence spectra of fluorescent targets, \({\mathbf {C}_A}\) and \({\mathbf {S}_A}\) denote the concentration and spectrum of background AF. \({\mathbf {D}_F} = \mathbf {D}{\mathbf {C}_A}{\mathbf {S}_A}\) and \(\mathbf {D_A} = \mathbf {D}{\mathbf {C}_F}{\mathbf {S}_F}\) represent the residual matrix. \(\alpha \in {{\mathbb {R}}^{+}}\) and \(\beta \in {{\mathbb {R}}^{+}}\) are scalars which weigh the contributions of the spatial regularization functions \(J(\mathbf {C}_F)\) of the matrix \(\mathbf {C}_F\) and \(J(\mathbf {C}_A)\) of the matrix \(\mathbf {C}_A\), respectively. In the real experiment, we assume that only one endmember is enough to represent background AF which includes all noise and other background signals, so that it means that \(\mathbf {C}_F \in \mathbb {R}_ + ^{N \times (K  1)}\), \({\mathbf {S}_F} \in \mathbb {R}_ + ^{(K  1) \times L}\), \({\mathbf {C}_A} \in \mathbb {R}_ + ^{N \times 1}\) and \({\mathbf {S}_A} \in \mathbb {R}_ + ^{1 \times L}\).
Because the update rule of HALS is the same as the classical multiplicative rule for NMF although the solution may be different. We can use the multiplicative update rule to optimize the problems of Eqs. (4) and (5) at last. Besides, many constraints can be easily introduced into NMF based on the scheme of multiplicative update rules, such that we can simply apply the appropriate constraints to optimize the solutions of \({\mathbf {C}_F}\) and \({\mathbf {C}_A}\).
Besides the convenience of introducing appropriate constraints to different hierarchies of endmembers, there is another important reason why we propose the twohierarchical structures for the solutions to subproblems of the Eq. (5). Although there are some classical NMF algorithms which can only add appropriate constraints to some specific endmembers, they don’t change the update rules of basic NMF framework but add constraints into some specific endmembers after the update of whole \(\mathbf {C}\). In fact, updating whole \(\mathbf {C}\) cannot prevent the classical NMF from mixing the target fluorophores and background AF in the real experiment. However, the hierarchical updating of Eqs. (6) and (7) for optimizing Eq. (5) is fully different from the updating framework of classical NMF, its advantage can be explained with the equality constraints in the work of Van Benthem et al. [19]. When updating the concentrations of multiple fluorescent targets \({\mathbf {C}_F}\), we assume that the concentration of AF \({\mathbf {C}_A}\) is fixed and known, and vice versa. In an iterative process of updating \(\mathbf {C}\), it seems that this assumption is a special “equality constraint” for the solution to \({\mathbf {C}_F}\) (or \({\mathbf {C}_A}\)). The equality constrained optimization is used to fix some components and to free other ones, such that it not only enables us to avoid unexpected components in the concentration matrix but also helps us gets more precise information on unknown concentrations. With an ongoing iterative process of hierarchical updating, we can obtain the desired concentrations of \({\mathbf {C}_F}\) and \({\mathbf {C}_A}\) according to different appropriate constraints. Furthermore, a more accurate \({\mathbf {C}_F}\) (or \({\mathbf {C}_A}\)) will lead to a more accurate \({\mathbf {C}_A}\) (or \({\mathbf {C}_F}\)) based on the hierarchical updating.
Constraints
Based on the basic framework of thNMF, different appropriate constraints can be imposed separately on each hierarchy of multiple fluorescent targets and background AF in invivo fluorescence imaging. As mentioned in the introduction section, a spatial sparsity constraint is appropriate for concentrations of multiple fluorescent targets with its improving the uniqueness of the decomposition along with enforcing a localbased representation. There are several popular ways to impose sparsity constraints on NMF. In Hoyer’s nonnegative sparse coding [27] framework, Frobenius norm is combined with \({{l}_{1}}\)norm of matrix \(\mathbf {C}\) or \(\mathbf {S}\) (the sparseness penalty term) to construct a cost function. Hoyer’s second method [28] has enforced a desired sparsity degree by means of nonlinear projection at each NMF iteration based on a sparsity measure from the relationship between \({{l}_{1}}\) and \({{l}_{2}}\)norms. With concentration matrix’s kth column being a nonzero N dimensional vector, the sparsity measure of \(\mathbf {C}_k\) is defined as
The first technique is conceptually simpler than the second technique but requires determining an appropriate penalty parameter \(\lambda\) which has a wide range of parameter values, while the second technique is more complicated but allows user to choose a priori sparsity value that ranges from 0 for nonsparse result to 1 for extremely sparse result. Being different from Hoyer’s second method enforcing the solution of matrix \(\mathbf {C}\) to reach a fixed sparsity value, Montcuquet et al. proposed a thresholding step in each updating iteration to help the algorithm converge to a more accurate solution [4].
For the framework of thNMF, it is easy to combine with different sparsity constraints mentioned above. For the simplicity and speed of iteration, we choose the \({{l}_{1}}\)norm as the sparsity measure function in the Eq. (6) to enforce the concentrations’ sparseness for multiple fluorescent targets, i.e., the Eq. (6) for the thNMF algorithm (see Algorithm 1) can be modified as
In the hierarchy of fluorescent targets, The updating rule of \(\mathbf {C}_F\) is similar to the method first proposed by Hoyer [27]:
The \({\alpha }\) is a small regularization parameter that balances the tradeoff between the approximation error and constraints. Compared with the sparsity value \(\phi\), the penalty parameter \(\alpha\) for \({{l}_{1}}\)norm is difficult to control. For the convenience of determining the penalty parameter and quantitative analysis of the sparsity degree among different algorithms, we use the similar method of [24] to adaptively control the \(\alpha\): \(\alpha\) for \(\mathbf {C}_F\) is initialized to 0.01, and after each iteration, \(\alpha\) is increased by 5 percent if average sparsity of \(\mathbf {C}_F\) is below the target sparsity \(\phi\), and is decreased by 5 percent otherwise. The desired sparsity value \(\phi\) is set to 0.65–0.85 for invivo fluorescence imaging. In the real application, it may occur that the \(\alpha\) is fluctuating during the iteration. To deal with this phenomenon, if the parameter \(\alpha\) fluctuates during the adjacent 20 iterations, the adaptive procedure of \(\alpha\) is terminated and \(\alpha\) becomes constant.
For the simplicity and overall performance of spectral unmixing, we only choose the appropriate constraint on the fluorescent targets rather than on the background AF though it is easy to introduce an appropriate constraint on the AF in Eq. (7).
During the iterative process, the entries of \(\mathbf {D}_F\) and \(\mathbf {D}_A\) are set to 0 if they are less than 0. At jth iteration, the value of cost function \(F{(\mathbf {C,S})^j}\) will be computed and the algorithm will be stopped if \(\left {F{{(\mathbf {C,S})}^j}  F{{(\mathbf {C,S})}^{j  1}}} \right < tol\), which is chosen equal to 0.0001.
Initialization of \(\mathbf {C}\) and \(\mathbf {S}\)
Due to NMF suffering from the illposed problem, the final decomposition solution strongly depends on NMF’s initialization. Moreover, for the NMF with different sparsity constraints on the different endmemebers, it is important to have relatively accurate initial estimates because specific endmembers correspond to specific sparsity degrees for the concentrations \(\mathbf {C}\). Some authors initialize \(\mathbf {S}\) using asymmetric Gaussian functions [14] or exponential initialization matrix [12] to simulate the quintessential fluorochrome spectra, while other works [15] exploit the statistical characteristics of the source spectra signals for the initialization of \(\mathbf {S}\). By adopting similar methods in the fluorescence imaging [4] and hyperspectral unmixing [29], we use calibration spectra to initialize the signature matrix \(\mathbf {S}\). The calibration spectra of multiple target fluorescence (fluorescent dyes) are acquired at the multiple emission filters by precalibration exvivo experiments in the same imaging conditions as the following invivo experiments, while the AF spectrum is the average spectrum acquired in some chosen regions of mouse with no fluorescent dyes. The concentration matrix \(\mathbf {C}\) is obtained by keeping spectra \(\mathbf {S}\) fixed during the first 10 iterations with the classical multiplicative update for NMF.
Results
Experimental results with synthetic data
In this section, we define a simulated experiment to quantitatively evaluate the performance of thNMF algorithm compared to other four classic spectral unmixing algorithms: multivariate curve resolutionalternating least squares (MCRALS) [3], HALSbased NMF with \({{l}_{1}}\)norm constraint on \(\mathbf {C}\) (NMF\({{l}_{1}}\)) [24], NMF with sparsity constraints on \(\mathbf {C}\) (NMFsc) [28], sparse NMF (sNMF) [4]. MCRALS has been used in the chemometrics of biological tissues [17] and small animal fluorescence imaging system [3], while NMF\({{l}_{1}}\), NMFsc and sNMF are different NMF algorithms with different sparsity constraints. According to the abovementioned initializations for \(\mathbf {C}\) and \(\mathbf {S}\), we use the same initial calibration spectra [4] and concentration matrix estimates for all NMFbased algorithms in all following experiments.
To test the proposed method’s performance on normal multispectral image data, we build a typical phantom [4] that is composed of two specific fluorescence objects and one strong AF object. The fluorescence spectra of specific fluorescence parts are fitted curves of Alexa Fluor (AF) 488 (see Fig. 1a) and 594 (see Fig. 1b) fluorescence dyes at their spectrum wavelength from 488 to 720 nm with interval of 5 nm, while the fluorescence spectrum of AF is a relatively slowvarying curve in the same range (see Fig. 1c). The concentration vectors of AF488 and AF594 consist of two parts: one part is pure fluorescent dye and other part is mixed with AF488 and AF594 (see Fig. 1a, b). The simulated fluorescent target signal (or AF signal) is a product of the concentration vector and the corresponding fluorescence spectrum. We can see that AF’s intensity in Fig. 1a, c has the same order of magnitude as that for the fluorescent target signal. Finally, a simulated phantom is obtained by adding the two specific fluorescent signals and the AF signal together (see Fig. 1d). At last, we provide a sensitivity study with the simulated data that are added with different zeromean white Gaussian noises with different signaltonoise ratios (SNRs).
We use spectral angle distance (SAD) [30] to evaluate the endmember spectrum estimation, and mean square error (RMSE) [31] to evaluate the concentration estimation. These two criteria are used to measure the similarity between the unmixed result and its standard (or ground truth) value. In general, the smaller value of the SAD (or RMSE) means the better similarity and the better unmixing performance. For the kth endmember, the SAD is defined as
where \(\mathbf {s}\) is the standard spectral vector and \(\hat{\mathbf {s}}\) its estimate. RMSE is defined as
where \(\mathbf {C}_{ki}\) and \({{\hat {\mathbf{C}}}_{ki}}\) are the standard and the estimated concentration of a specific pixel, respectively. At last, the average value of all endmembers’ SAD (or RMSE) is used to evaluate the performance of estimating fluorescence spectra and concentrations.
In this synthetic experiment, for all algorithms the parameter settings are as follows: an unimodality of spectra is introduced in 1th and 2nd endmember for MCRALS; the penalty coefficient of sparsity constraint is 0.01 for NMF\({{l}_{1}}\); the sparsity of concentration distribution is 0.65 at NMFsc; the sparsity value for both the first and second concentrations is 0.73 at sNMF. The number of endmember for all algorithms is 3.
Figure 2 shows the ground truth and unmixed fluorescence spectra of synthetic data when SNR = 30 dB. The unmixed AF spectrum for MCRALS in Fig. 2b is totally wrong and the unmixed spectrum of AF488 is somewhat dissimilar to the ground truth (Fig. 2a). All unmixed fluorescence spectra of NMFsc in Fig. 2d have large artifacts. The fluorescence spectra of AF488 and AF594 for NMF\({{l}_{1}}\) in Fig. 2c and for sNMF in Fig. 2e have abnormal offsets and distortions between 480 and 650 nm. All unmixed endmembers by our thNMF method (Fig. 2f) obtain the best fluorescence spectra that are most similar to the ground truth of spectra.
The normalized unmixing results of all endmembers obtained with all algorithms are displayed in Fig. 3, where the unit of the color bar is normalized from 0 (zero intensity) to 1 (the maximum intensity of the final unmixing results obtained with all methods). The unmixed AF488 and AF594 for MCRALS in Fig. 3b and NMFsc in Fig. 3d are not localized in spatial distributions, such that they are very dissimilar to the ground truth of simulated AF488 and AF594 fluorescences in Fig. 3a. There are obvious holes in the spatial distributions of unmixed AF by NMF\({{l}_{1}}\) in Fig. 3c, NMFsc in Fig. 3d and sNMF in Fig. 3e. Moreover, from the color scale representations of unmixing results, we can see that the concentrations of AF488 and AF594 for MCRALS, NMF\({{l}_{1}}\), NMFsc and sNMF are much larger than the ground truth value of concentrations in Fig. 3a. Though sNMF can get desirable sparse spatial distribution for AF488 and AF594, the unmixed concentrations are larger than the ground truth concentrations of these two fluorescence targets; besides, the unmixed AF have two holes in the locations of fluorescence targets. Compared with other methods, only our thNMF method in Fig. 3f achieves the best unmixing results that are very similar to the ground truth of simulated data.
Figure 4 manifests the mean values of SAD and RMSE as functions of the SNRs ranging from \(\infty\) (without any noise) to 15 dB in steps of 5 dB for three endmembers. We can see that the NMFsc obtains the worst results for both the SAD and the RMSE. Meanwhile, the proposed thNMF has not only the best mean SAD but also the smallest mean RMSE to demonstrate its robustness against noise corruption.
Experimental results with real fluorescent data
In this section, we apply thNMF algorithm to invivo fluorescent data acquired with a BioReal Quick View 3000 imaging system (BioReal Sciences, Austria), which is an epiillumination planar imaging system equipped with electron multiplying charge coupled device (EMCCD) and 150 watt Xenon arc lamp as a light source. The EMCCD (Andor Technology, Ireland) has high quantum efficiency in the range of 400–950 nm which is our interest of fluorescence imaging. In the epiillumination fluorescence imaging system, the EMCCD gets the emission light through the multiple bandpass (30 nm) emission filters covering 525, 542, 579, 624 and 716 nm (Semrock, Rochester, USA).
In the following three invivo experiments, we use one nude and two BALB/c mice with subcutaneous specific fluorescent dyes to confirm the effectivity of thNMF for preclinical applications. This subcutaneous injection of dye is sometimes adopted in clinical fluorescence unmixing researches [2–4] and can easily spread the dye through mouse’s body. The animal experiments were approved by our institutional review board. The nude mouse has weak AF and is often used in fluorescence imaging, while the BALB/c mouse has much stronger AF than nude mouse due to the fur of BALB/c mouse being easy to be excited. Furthermore, two fluorescent dye pairs AF488/AF594 and AF488/AF555 are used as fluorescent markers, in which the AF488/AF555 pair has more spectral overlap between two dyes’s spectra than the AF488/AF594 has (see the following sections for detail). Both of these multipletarget fluorescence emissions have spectral overlap with the spectrum of complex background AF. To reduce the AF in stomach region, the mouse may be fed with alfalfafree diet in preclinical applications. However, the mice are normally fed before invivo experiments in this work. To demonstrate the raw images and unmixing results clearly, all the observations and unmixing results are displayed in pseudocolor mapping based on logarithmic scaling (or linear scaling) and overlaid with the corresponding grayscale photographic image of mouse.
We demonstrate four classic spectral unmixing algorithms (MCRALS, NMF\({{l}_{1}}\), NMFsc and sNMF mentioned in the above section) for comparison purpose. All algorithms assume that there are three endmembers for spectral unmixing, or K=3. The parameter settings for all algorithm are as follows: an unimodality constraint of spectra is introduced in the first and second endmember for MCRALS; the penalty coefficient of sparsity constraint is 0.001 for NMF\({{l}_{1}}\); the sparsity of concentrations is 0.88 in NMFsc; the sparsity value for both the first and second concentrations is 0.87 for sNMF.
Unlike the invitro fluorescence, the photons from the invivo target fluorophores are located inside tissues and may be greatly attenuated before reaching the mouse surface due to light scattering and absorption. This attenuation varies with fluorescent dye’s injection depth, tissue type, and wavelength, it is hence impossible to obtain the true concentrations of all endmembers in invivo multiple target fluorophores. Besides, the true spectrum of AF is unknown in the real invivo experiments. Therefore, we use SAD to evaluate the spectrum estimation of endmember AF488 and AF594. Besides, the visual inspection is still a useful assessment method to evaluate algorithm’s performance. We also introduce a sparsity value \(\phi\) of column vector of concentration matrix \(\mathbf {C}\) [27] to evaluate the concentration estimation. In the real experiments, the sparsity values of fluorescent dyes \(\mathbf {C}_F\) must have relatively large values because of the spatial sparseness of fluorescence targets while the sparsity value of AF \(\mathbf {C}_A\) must have a relatively small value because of the dispersive distribution of AF.
The nude mouse
In the first two experiments, the AF488 and AF594 fluorescent dyes are diluted to 0.1 μg ml\(^{1}\). The raw images are sparsely acquired through five emission glass filters from 525 to 716 nm. The first three images are excited at 474 nm, the last two images are excited at 565 nm. In the first nude mouse experiment, AF488 is injected at the bottom of the body with 20 ng dye while AF596 near the neck with the same quantity, and a mixture of each dye with 10 ng is located at the middle of the body.
Figure 5(a1)–(a5) display the raw spectral images acquired from five emission filters. The pseudocolor mapping is based on logarithmic scaling to clearly display the detailed change of weak AF in all spectral channels. The AF488 is displayed on the channels 525, 542 and 579 nm while the AF594 on the 624, 716 nm. It seems that there is almost no crosstalk between these two fluorescent dyes. However, the AF is weakly presented on the background surface of mouse and overlapped with AF488 and AF594 in all multispectral observations. Therefore, this nude mouse experiment is appropriate to test the unmixing performance in removing weak AF from the multiple fluorescence targets.
Five algorithms’ unmixing results for invivo fluorescence imaging are displayed in Fig. 5(b1)–(b3) for MCRALS, Fig. 5(c1)–(c3) for NMF\({{l}_{1}}\), Fig. 5(d1)–(d3) for NMFsc, Fig. 5(e1)–(e3) for sNMF, Fig. 5(f1)–(f3) for our thNMF. MCRALS and NMF\({{l}_{1}}\) methods still overlap the fluorescent dyes with background AF. Though the unmixed fluorescent dyes obtained with NMFsc and sNMF are separated to some degree, the unmixing is not complete and the fluorescent dyes with their strong intensities are still overlapped with the AF in Fig. 5(d3) and Fig. 5(e3). Only the proposed thNMF can successfully separate multiple fluorescent dyes and AF.
Table 1 quantifies all algorithms’ unmixing performances using the SAD similarity, while Table 2 represents the sparsity values of concentrations for all endmembers. Table 1 shows that the unmixing results of thNMF get the smallest value of SAD for the AF488 endmember and the smallest value of mean SAD. The proposed thNMF also obtains the second smallest value of SAD for the AF594. Though MCRALS gets the smallest value of SAD for the AF594 endmember, the sparsity value 0.8077 of concentration for AF achieved with MCRALS is larger than the sparsity value 0.6138 achieved with thNMF, which is in conflict with the reality that the sparsity value of fluorescent dye is much larger than that of AF. Considering these factors comprehensively, we confirm that the proposed thNMF gets the best unmixing results with fully removing the AF for the nude mouse experiment.
The BALB/C mouse 1
The experimental condition for BALB/c mouse is similar to previous nude mouse experiment. Fig. 6(a1)–(a5) show the raw spectral images acquired from five emission filters. The AF is apparently overlapped with the AF488 at the first three channels and the AF594 at the last two channels, respectively. The BALB/C mouse introduces more AF and more attenuation of fluorescence signal than the nude mouse dose due to BALB/C mouse hair, which further lowers the contrast between the target fluorescence and background AF in BALB/C mouse compared to the nude mouse. The AF distribution is more irregular and large than that of nude mouse. Moreover, BALB/C mouse’s AF intensity is much stronger than nude mouse’s AF intensity. Especially, the AF at the upper right part of BALB/C mouse has an intensity comparable to the AF488 intensity at the channels of 525, 542 and 579 nm.
Five algorithms’ unmixing results are illustrated in Fig. 6(b1)–(b3) for MCRALS, Fig. 6(c1)–(c3) for NMF\({{l}_{1}}\), Fig. 6(d1)–(d3) for NMFsc, Fig. 6(e1)–(e3) for sNMF, Fig. 6(f1)–(f3) for our thNMF. The unmixed fluorescent dyes obtained with MCRALS, NMF\({{l}_{1}}\) and NMFsc are still overlapped with background AF, while the unmixed fluorescent dyes obtained with NMFsc in Fig. 6(d1)–(d3) are still overlapped at different locations. Though the unmixing results of sNMF (Fig. 6(e1)–(e3)) are similar to thNMF’s results (Fig. 6(f1)–(f3)), the concentration distribution of AF is affected by the fluorescent dyes in the middle of body such that there is an obvious hole in AF at the bottom of body in Fig. 6(e3). Among all unmixing results from all five algorithms, the unmixing result of thNMF is best and the crosstalk between AF and the two fluorescent dyes is reduced to a negligible level.
Table 3 quantifies the different fluorescence unmixing performances using the SAD similarity, while Table 4 represents the sparsity values of concentrations for all endmembers. The evaluation results in Table 3 are similar to those of nude mouse experiment. Our thNMF obtains the smallest mean SAD, the smallest SAD value of AF488 and the second smallest SAD value of AF594, while sNMF gets the smallest SAD value of AF594 endmember. All SAD values of AF594 endmember are close to each other for all algorithms except NMF\({{l}_{1}}\) algorithm.
Table 4 displays that the sparsity value of AF for NMF\({{l}_{1}}\) is larger than that of AF488. Besides, the sparsity value of AF for NMFsc is equal to that of AF488 (or AF594). These sparsity results simply do not conform to the reality that the sparsity value of AF is much smaller than that of fluorescent dye. However, thNMF and sNMF can achieve good sparsity value for all concentrations in Table 4. More precisely, thNMF has better sparsity values of AF488 and AF than sNMF dose. Taking all these factors into account, we can find that thNMF can obtain the best unmixed result as compared with other four algorithms in the experiment of BALB/c mouse.
The BALB/C mouse 2
In the last experiment, the AF488 and AF555 used as two target fluorophores are diluted to 0.1 μg ml\(^{1}\). AF488 (50 ng) is subcutaneously injected at the bottom right of the mouse body and AF555 (50 ng) on the left and right sides of the middle part of body. Fig. 7(a1)–(a4) show raw fluorescence images that are sparsely acquired at the 525, 542, 579, 624 nm emission filters. Using linear scaling for the pseudocolor mapping can clearly demonstrate the detailed change of strong AF, which has clear spectral overlap with AF488 and AF555 in all spectral channels. The AF488 is excited at 474 nm and displayed nearly on all the four channels while the AF555 is excited at 500 nm and emitted on the 579 and 624 nm channels. In addition to the strong AF overlapping greatly with the AF488 and AF555, there is also distinct crosstalk between these two target fluorophores.
We compare the proposed method to MCRALS, NMF\({{l}_{1}}\) and sNMF that have good performances in previous experiments. The unmixing results are illustrated in Fig. 7(b1)–(b3) for MCRALS, Fig. 7(c1)–(c3) for NMF\({{l}_{1}}\), Fig. 7(d1)–(d3) for sNMF, Fig. 7(e1)–(e3) for thNMF. Except for the MCRALS, the other three methods, i.e. NMF\({{l}_{1}}\), sNMF and thNMF, can separate the fluorescent targets in Fig. 7(c1)–(c2), (d1)–(d2), (e1)–(e2) from AF, respectively. But there are missing parts in AF component in Fig. 7(c3) for NMF\({{l}_{1}}\) and Fig. 7(d3) for sNMF. Though the sparse acquisition from only four spectral channels exacerbates the illposed problem of spectral unmixing, the proposed thNMF can fully remove the intense AF from the two target fluorophores and obtain best unmixing results among all algorithms. Besides these visual assessment, Table 5 confirms that the proposed algorithm achieves the smallest values of average SAD.
Discussion and conclusions
By spatially distinguishing the target fluorophores from background AF within equality constrained [19] HALS framework, we have proposed thNMF algorithm for invivo multispectral fluorescence imaging to divide the whole update process of concentration matrix \(\mathbf {C}\) into two relevant subprocesses, at which the concentrations of target fluorophores and background AF are separately updated with a spatial sparsity constraint being solely imposed on the target fluorophore. A detailed comparison with other four stateofarts spectral unmixing algorithms is conducted to confirm the excellent performance of thNMF algorithm. In NMFbased spectral unmixing, the initial estimates of the endmembers and AF are important for solving the illposed problem of NMF. As in [4], all the methods for comparison in this work use calibration spectra to initialize the spectral matrix \(\mathbf {S}\). This initialization alone can not guarantee NMF’s convergence to an ideal solution. Furthermore, the very sparse acquisition of multispectral fluorescence imaging has incomplete measurement and severe discontinuity problems over multispectral fluorescence emissions with various complex AF, such that it aggravates the illposed nature of the inverse problem of spectral unmixing. Under such sparse acquisition, the unimodality of spectra constraint (for MCRALS) and other sparsity constraints for some statesofarts NMFbased spectral unmixing algorithms can not facilitate achieving the optimal unmixed solutions. However, by employing equality constraint along with partial sparsity constraint during HALS optimization where fixing background AF (resp. sparsity constrained target fluorophores) can make the update of sparsity constrained target fluorophores (resp. background AF) as accurate as possible, the thNMF algorithm can successfully tackle the illposed and sparse acquisition problems in the NMF framework according to the different spatial distribution of the background AF and the target fluorophores.
Considering the different spatial characteristics of endmembers from different hierarchies in the invivo fluorescence imaging, the proposed thNMF algorithm can easily add different spatial regularization constraints into different hierarchies to get optimal solution. Although we concentrate on sparse multispectral imaging and fluorescence unmixing, twohierarchical updating that distinguishes the target from background signal can be an interesting inspiration for research of multispectral fluorescence imaging and blind source separation. This is because most multimodality signals can be decomposed into two hierarchies: one hierarchy of localized signals of interest with its sparse representation and another hierarchy of background noisy signal with its widely distributed spatial existence.
Abbreviations
 NMF:

nonnegative matrix factorization
 SUM:

spectral unmixing (SUM)
 thNMF:

twohierarchical NMF
 HALS:

hierarchical alternating least squares
 ANLS:

alternating nonnegative least squares.
References
 1.
Zimmer M. GFP: from jellyfish to the Nobel prize and beyond. Chem Soc Rev. 2009;38:2823–32.
 2.
Leblond F, Davis SC, Valdes PA, Pogue BW. Preclinical wholebody fluorescence imaging: review of instruments, methods and applications. J Photochem Photobiol B Biol. 2010;98(11):77–94.
 3.
Xu H, Rice BW. Invivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique. J Biomed Opt. 2009;14(6):064011.
 4.
Montcuquet AS, Herve L, Navarro F, Dinten JM, Mars JI. In vivo fluorescence spectra unmixing and autofluorescence removal by sparse nonnegative matrix factorization. IEEE Trans Biomed Eng. 2011;58(9):2554–65.
 5.
Huang Z, Zeng H, Hamzavi I, Alajlan A, Tan E, Mclean DI, Liu H. Cutaneous melanin exhibiting fluorescence emission under nearinfrared light excitation. J Biomed Opt. 2006;11(3):034010.
 6.
Montcuquet AS, Herve L, Navarro F, Dinten JM, Mars JI. Nonnegative matrix factorization: a blind spectra separation method for in vivo fluorescent optical imaging. J Biomed Opt. 2010;15(5):056009.
 7.
Shi C, Wang L. Incorporating spatial information in spectral unmixing: a review. Remote Sens Environ. 2014;149:70–87.
 8.
GutierrezNavarro O, CamposDelgado DU, ArceSantana ER, Mendez MO, Jo JA. Blind endmember and abundance extraction for multispectral fluorescence lifetime imaging microscopy data. IEEE J Biomed Health Inform. 2013;18(2):606–17.
 9.
Fereidouni F, Blab GA, Gerritsen HC. Blind unmixing of spectrally resolved lifetime images. J Biomed Opt. 2013;18(8):086006.
 10.
Paatero P, Tapper U. Positive matrix factorization: a nonnegative factor model with optimal utilization of error estimates of data values. Environmetrics. 1994;5:111–26.
 11.
Lee D, Seung H. Algorithms for nonnegative matrix factorization, Proc Adv Neural Inf Process Syst 2000 Conf. 13; 2001. pp. 556–62.
 12.
Woolfe F, Gerdes M, Bello M, Tao X, Can A. Autofluorescence removal by nonnegative matrix factorization. IEEE Trans Image Proc. 2011;20(4):085–1093.
 13.
Pengo T, MunozBarrutia A, Zudaire I, de Solorzano CO. Efficient blind spectral unmixing of fluorescently labeled samples using multilayer nonnegative matrix factorization. PLOS One. 2013;8(11):e78504.
 14.
Neher RA, Mitkovski M, Kirchhoff F, Neher E, Theis FJ. Blind source separation techniques for the decomposition of multiply labeled fluorescence images. Biophys J. 2009;96:3791–800.
 15.
Yu S, Zhang Y, Liu W, Zhao N, Xiao X, Yin G. A novel initialization method for nonnegative matrix factorization and its application in component recognition with threedimensional fluorescence spectra. Spectrochim Acta A Mol Biomol Spectrosc. 2012;86:315–9.
 16.
Kim H, Park H. Sparse nonnegative matrix factorizations via alternating nonnegativityconstrained least squares for microarray data analysis. Bioinformatics. 2007;23(12):1495–502.
 17.
Pigueras S, Krafft C, Beleites C, Egodage K, von Eggeling F, GuntinasLichius O, Popp J, Tauler R, de Juan A. Combining multiset resolution and segmentation for hyperspectral image analysis of biological tissues. Anal Chim Acta. 2015;881:24–36.
 18.
Li Q, He X, Wang Y, Liu H, Xu D, Guo F. Review of spectral imaging technology in biomedical engineering: achievements and challenges. J Biomed Opt. 2013;18(10):100901.
 19.
Benthem MHV, Keenan MR, Haaland DM. Application of equality constraints on variables during alternating least squares procedures. J Chemom. 2002;16(12):613–22.
 20.
Cichocki A, Zdunek R, Amari S. Hierarchical ALS algorithms for nonnegative matrix and 3D tensor factorization. In: Davies ME, James CJ, Abdallah SA, Plumbley MD, editors. Independent component analysis and signal separation. Lecture notes in computer science, vol. 4666. Berlin: Springer; 2007. pp. 169–176.
 21.
Gillis N, Glineur F. Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization. Neural Comput. 2011;24:1085–105.
 22.
Ho ND. Nonnegative matrix factorizationalgorithms and applications. Ph.D. thesis, Universite Catholique de Louvain, LouvainlaNeuve, Belgium; 2008.
 23.
Van Beek J, Tits L, Somers B, Deckers T, Janssens P, Coppin P. Reducing background effects in orchards through spectral vegegation index correction. Int J Appl Earth Observ Geoinform. 2015;34:167–77.
 24.
Gillis N. Nonnegative matrix factorization: Complexity, algorithms and applications. Ph.D. thesis, Univ. Catholique de Louvain, LouvainlaNeuve, Belgium; 2011.
 25.
Jia S, Qian Y, Li J, Li Y, Ming Z. Hierarchical alternating least squares algorithm with sparsity constraint for hyperspectral unmixing. In: Proceedings of 2010 IEEE 17th Int. Conf. Image Processing; 2010. pp. 2305–2308.
 26.
Chen W, Guillaume M. HALSbased NMF with flexible constraints for hyperspectral unmixing. EURASIP J Adv Signal Process. 2012;1:1–14.
 27.
Hoyer P. Nonnegative sparse coding. In: Proceedings of the 12th IEEE Workshop on Neural Networks for Signal Processing; 2002. pp. 557–565.
 28.
Hoyer P. Nonnegative matrix factorization with sparseness constraints. J Mach Learn Res. 2004;5(9):1457–69.
 29.
Tang W, Shi Z, An Z. Nonnegative matrix factorization for hyperspectral unmixing using prior knowledge of spectral signatures. Opt Eng. 2012;51(8):087001.
 30.
van der Meer F. The effectiveness of spectral similarity measures for the analysis of hyperspectral imagery. Int J Appl Earth Obs Geoinf. 2006;8(1):3–17.
 31.
Plaza A, Martinez P, Perez R, Plaza J. A quantitative and comparative analysis of endmember extraction algorithms from hyperspectral data. IEEE Trans Geosci Remote Sens. 2004;42(3):650–63.
Authors’ contributions
SSH, YZ & BJQ participated in literature search, image analysis, manuscript writing and editing. All the authors read and approved the final manuscript.
Acknowledgements
This work was supported in part by the National Natural Science Foundation of China (61271320 and 60872102) and the Small Animal Imaging Project (No. 06545). The authors would like to thank Dr. Yonggang Zhao for inspiring discussion, Liming Tan for contributing to the project, and Dr. Jinliang Pen for supporting the experiments in fluorescence imaging. We are thankful to the anonymous reviewers for their valuable comments that greatly helped to improve this paper.
Competing interests
The authors declare that they have no competing interests.
Author information
Additional information
Shaosen Huang and Yong Zhao are contributed equally to this work
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Huang, S., Zhao, Y. & Qin, B. Twohierarchical nonnegative matrix factorization distinguishing the fluorescent targets from autofluorescence for fluorescence imaging. BioMed Eng OnLine 14, 116 (2015). https://doi.org/10.1186/s1293801501074
Received:
Accepted:
Published:
Keywords
 Multispectral fluorescence imaging
 Equality constraint
 Spectral unmixing
 Autofluorescence
 Nonnegative matrix factorization
 Twohierarchical