Anne-Sophie Montcuquet

Lionel Hervé

Fabrice Navarro

Jean-Marc Dinten

Jérôme I Mars

In vivo fluorescence spectra unmixing and autofluorescence removal by sparse Non-negative Matrix Factorization

Keywords: Source separation, fluorescence optical imaging, autofluorescence

Fluorescence imaging locates fluorescent markers that specifically bind to targets, as tumors: markers are injected to a patient, optimally excited with near infrared light, and located thanks to emitted back fluorescence analysis. To investigate thick and diffusive media, as the fluorescence signal decreases with the light travel distance, the autofluorescence of biological tissues comes to be a limiting factor. To remove autofluorescence and isolate specific fluorescence, a spectroscopic approach, based on Non-negative Matrix Factorization (NMF), is explored. To improve results on spatially sparse markers detection, we suggest a new constrained NMF algorithm which takes sparsity constraints into account. A comparative study between both algorithms is proposed on simulated and in vivo data.

I. INTRODUCTION

Fluorescence spectroscopy and Fluorescence diffuse optical tomography (FDOT) are promising non-invasive, user-friendly and cost effective methods to serve medical diagnostic systems. Fluorescent compounds, called markers, are injected to a patient, and specifically bind to targeted tumors [START_REF] Ntziachristos | Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin v-cy5. 5 conjugate[END_REF], [START_REF] Leblond | Pre-clinical whole-body fluorescence imaging: Review of instruments, methods and applications[END_REF], [START_REF] Bremer | Optical-based molecular imaging: contrast agents and potential medical applications[END_REF], [START_REF] Ballou | Fluorescence imaging of tumors in vivo[END_REF]. This method of targeting allows researchers to detect and localize cancer cells in patients. Once injected, markers are excited by near infrared (NIR) light -between 600 and 900 nm -where the tissue absorption is lower. The excitation wavelength is chosen to ease the biological tissue penetration, and to optimally excite the injected markers. The fluorescence signal emitted back by markers is measured, and the cancer cells may then be localized. So far, NIR fluorescence imaging is mainly used on small animals where some markers are available for injection, and where the layer of biological tissues to explore does not exceed a few centimeters. In such case, a biological tissues intrinsic fluorescence -called autofluorescence [START_REF] Hung | Autofluorescence of normal and malignant bronchial tissue[END_REF] -exists, but is insignificant compared to the fluorescent marker signal and does not prevent localization of the cancer cells. For medical diagnostic application on thick media (human prostate, or breast for example, that is around 4 cm), the autofluorescence is worth considering: as the fluorescence signal gets exponentially weak with the light travel distance, the autofluorescence signal remains constant and turns into a limiting factor. The analysis of fluorescence signals impaired by autofluorescence may lead to a wrong localization of the markers and consequently to a wrong localization of cancer cells: the signal needs to be preprocessed in order to remove autofluorescence.

Several methods to unmix fluorescence spectra and to filter autofluorescence have already been developed and tested on small animal examination equipments. Among those methods are the non linear least squares [START_REF] Wood | Refining epifluorescence imaging and analysis with automated multiple-band flat-field correction[END_REF], principal component analysis (PCA), independent component analysis (ICA) and singular value decomposition (SVD) [START_REF] Levenson | Spectral imaging of deep tissue[END_REF], [START_REF] Hoyt | Systems and methods for in-vivo optical imaging and measurement[END_REF], [START_REF] Cronin | Spectral imaging[END_REF], [START_REF] Lathauwer | Blind source separation by higher-order singular value decomposition[END_REF] methods. The ICA method requires sources to be statistically independent, SVD considers orthogonal sources, and many other a priori knowledges about the nature of the sources are taken into account in those methods. But a principal one remains missing: non-negativity [START_REF] Henry | Current factor analysis receptor models are ill-posed[END_REF]. Many real-world data are non-negative, among are the fluorescence data, and the unmixed fluorescence spectra provided by a separation source method have a physical meaning only when non-negative. In light of this observation, some studies began to consider Nonnegative Matrix Factorization (NMF) to separate fluorescence spectra.

NMF was notably investigated by Paaetero and Tapper, and gained popularity in 2001, through the works of Lee and Seung [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. In 2004, Gobinet used the NMF decomposition on spectroscopic data to unmix several pure fluorescence spectra and get a fluorescence chemical mapping of wheat grain sections [START_REF] Gobinet | Application of nonnegative matrix factorization to fluorescence spectroscopy[END_REF]; in 2007, they used NMF to dewax Raman signals of human skin biopsies [START_REF] Gobinet | Preprocessing and source separation methods for Raman spectra analysis of biomedical samples[END_REF]. NMF was also applied to fluorescence spectroscopy, to unmix several markers by Xu et al. [START_REF] Xu | Improved sensitivity by applying spectral unmixing prior to fluorescent tomography[END_REF], [START_REF] Xu | In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique[END_REF]. Their work is the closest one from our research, but is still applied only to small animal in vivo imaging, while we try to focus on models and tumor/healthy skin ratios close to deep embedded markers problems. For in vivo fluorescence spectroscopy, the unmixing problem is referred to as a blind source separation problem since the spectra may vary according to the fluorescent dye biological environment. Fluorescence spectra to separate are also supposed statistically dependent, which filters out many methods (such as ICA). Finally, the NMF algorithm seems to be, in many ways, more suitable in blind positive spectra separation than other separation methods. Several NMF algorithms exist, based on diverse criteria to minimize and optimization methods: we present in the paper a classical NMF algorithm that deals with multiplicative update rules [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. But as all blind source separation methods, it is impossible to find a unique NMF decomposition [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-way Data Analysis and Blind Source Separation[END_REF]. Without any constraint, there is an infinity of decomposition solutions. By considering only non-negative solutions, which is inherent to the NMF method, we already restrain the solution set. To relax more the non-uniqueness issue and to refine the solution set, we suggest to apply appropriate regularizations and a priori knowledge considerations. Such additional constraints usually deal with sparsity and smoothness constraints [START_REF] Cichocki | New algorithms for nonnegative matrix factorization in applications to blind source separation[END_REF], [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]. Kim and Park got for example interested in sparse NMFs algotihms [START_REF] Kim | Non-negative matrix factorization based on alternating non-negativity constrained least squares and active set method[END_REF] by L 1 -norm constraint term minimization, Cichocki et al. have presented cost functions based not anymore on the Kullback-Leibler divergence but on Csiszár's ϕ-divergence [START_REF] Cichocki | Independent Component Analysis and Blind Signal Separation[END_REF] while other approaches use alternative cost functions formulations [START_REF] Hamza | Reconstruction of reflectance spectra using robust nonnegative matrix factorization[END_REF], [START_REF] Dhillon | Generalized nonnegative matrix approximations with Bregman divergences[END_REF]. For fluorescence imaging application, we propose a regularized NMF algorithm that takes spatial sparsity constraints into account to improve NMF decomposition and tested on spectroscopic simulated and in vivo data.

In this article, we first introduce the NMF method, and one of the most popular associated algorithm, proposed by Lee and Seung in 2001 [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. From that classical algorithm, to specifically detect spatially sparse fluorescent sources, we defined a new constrained algorithm with sparsity constraints. A comparative study is then run in a second part to compare NMF results obtained with or without taking sparsity constraints into account. To illustrate that study, a simulated example on which several tests are run is proposed. Finally, we present experimental in vivo data on the last part. Up to three fluorescent sources with overlapping emission spectra have to be unmixed: two fluorescent markers -ICG loaded into nanoparticules (ICG-LNP) and Alexa 750 -plus the autofluorescence signal. Once more, unmixing results obtained with the NMF algorithm and the sparse NMF algorithm are compared.

II. THEORY

A. Non-Negative Matrix Factorization

For initial non-negative mixed data M , NMF proposes to find a couple of matrices (A, S) with non-negatives coefficients, whose product optimally approaches M . The classical NMF definition says [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]:

Given a non-negative matrix M ∈ R X×Y , find non-negative matrices A ∈ R X×P and S ∈ R P ×Y such that :

M AS (1)
where non-negative matrices are matrices whose all factors are nonnegative and P stands for the number of sources to unmix.

Applied to spectroscopy, matrix A is considered as the weights matrix, and S as the spectra matrix; both matrices A and S contain respectively as much columns and lines P as fluorescent sources to separate.

To find the best matrices A and S that satisfy equation 1, two distinct steps must be considered. First, a criterion that links M with A and S has to be defined. Then, in a second step, the criterion has to be optimized. Optimizing this criterion -generally by minimizing it -under the nonnegativity constraint would lead to the best couple of solutions (A, S). Several criteria (the Euclidean distance, the Kullback-Leibler divergence [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]...) and different optimization methods (Alternating Least Square, multiplicative update rules) may suit to NMF method.

B. Sparsity constraints

Besides the two basic steps -criterion choice and optimization -that lead to the NMF decomposition, additional constraints may be considered. Those constraints may directly modify criterion to optimize, or would imply a new optimization step.

Several studies got interested in sparsity constraints, applied to NMF. Among them we should name Hoyer et al. [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] and Stadlthanner et al. [START_REF] Stadlthanner | Extended sparse nonnegative matrix factorization[END_REF] who defined new sparse NMF algorithms. Both research groups got interested in sparsity of matrix S, while no assumption was made on A. But we can easily extend these algorithms in order to constrain matrix A. A difference remains between both algorithms: as Hoyer constrained all rows of S to have a common sparsity value, Stadlthanner extended Hoyer's algorithm by providing a different sparsity coefficient to each row of matrix S. This change makes the extended sparse NMF (esNMF [START_REF] Stadlthanner | Extended sparse nonnegative matrix factorization[END_REF]) algorithm more suitable for blind source separation (BSS) problems, where sources may have different sparsenesses. We implemented these algorithms in order to test them on our data. In our case, we got interested in sparseness of matrix A, and considered that sparseness of all columns of A could differ from a source to another. Sparse columns of A in spectroscopy is the expression of spatially sparse fluorescence signals. Indeed, for local specific markers distributions, we expect to get peaked and sparse weight columns in matrix A. A contrario, when unmixing is not complete, a residual autofluorescence baseline surrounds the specific marker peak. This property is depicted Figure 1.

Marker unmixing result

Weight matrix A specific markers

Inaccurate unmixing Correct unmixing

Sparse data Non-sparse data An intuitive reasoning to improve unmixing would be to smooth over the unwanted autofluorescence residuals on spe-cific markers weight profiles, by thresholding the smallest values: in other words, we look for sparse A columns for specific markers contributions. We introduce a sparsity value [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]: let us consider a weight matrix A of size X × P , the sparsity of a given column A p (p ∈ (1, P )) of A is:

sparsity (A p ) = √ X - X x=1 |a xp | X x=1 a 2 xp √ X -1 (2)
Sparsity value ranges from 0 for non-sparse results to 1 for extremely sparse results, as depicted Figure 2. In the next section, we propose an NMF algorithm that finds the perfect threshold so that our matrix A sparsity remains as close as possible to a chosen sparsity value.

III. MATERIAL AND METHODS

A. Simulated example

In order to test our unmixing method, we define a simulated example to run the NMF algorithms on. We propose to build a simulated phantom, composed of an autofluorescence part and of a specific fluorescence one (to compare to a tumor pointed out by an injected fluorescent marker). The specific fluorescence part is the product of a weight vector A 1 by a fluorescence spectrum S 1 (see Figure 3 -a). A pari, the autofluorescence part is the product of the weight vector A 2 by the fluorescence spectrum S 2 (see Figure 3 -b). Finally the total simulated phantom is obtained by adding the specific fluorescence and the autofluorescence parts (Figure 3 -c). In fluorescence optical imaging, signal ratio between healthy tissue and tagged tumor depends on biomarkers injected.

From bibliography, and from experience, ratios from 3 to 15 [START_REF] Jin | In vivo optical imaging of integrin alphav-beta3 in mice using multivalent or monovalent crgd targeting vectors[END_REF] (for more specific-to-tumor markers) are usual. New generation of activatable fluorescent markers reach ratios from 24 to 180 [START_REF] Urano | Selective molecular imaging of viable cancer cells with pH-activatable fluorescence probes[END_REF], [START_REF] Weissleder | In vivo imaging of tumors with protease-activated near-infrared fluorescent probes[END_REF], [START_REF] Razkin | Activatable fluorescent probes for tumour-targeting imaging in live mice[END_REF], depending on wavelength range observed, and conditions of experimentation: ex vivo or in vivo, and localization of tumor. For this simulation example, to reproduce deep embedded markers detection, we chose a ratio tumor/healthy tissue with highest intensity value approximately equal to 1.

B. Contrast definition

In order to evaluate the tumor detection, we introduce the contrast C T,N which is measured between a tumorous area T and a normal (or healthy) tissues area N . The value of C T,N characterizes the tumor detection after autofluorescence removal, on simulation and experimental results. Average intensity of fluorescence signal is measured on both concerned ROIs: T and N are respectively the average intensities in photons per pixel of areas T and N on analyzed images (see Figure 3-c).

c T,N = T - N T + N (3)
The closer to one the contrast value gets, the better the detection will be.

C. Classical NMF algorithm

In 2001, NMF popularity increased after Lee and Seung published two new NMF algorithm, based on the use of multiplicative update rules that minimize specific criteria.

The NMF decomposition looks for the best solution couple (A, S) whose product best approaches the initial data V . Classically, to find matrices A and S, a chosen criteria is iteratively minimized.

Different criteria to minimize, or cost functions, can be used: we may cite the square of the euclidean distance between V and AS, and the Kullback-Leibler divergence criterion [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. Here, we define the cost function F to minimize as the square of the Euclidean distance between V and AS [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], lower bounded by 0:

F = Nx x=1 Ny y=1 (v xy - P p=1 a xp s py ) 2 = V -AS 2 2 (4)
Here, . 2 denotes the standard Euclidean norm. In order to get solutions (A, S), the following optimization problem is thus considered:

Problem 1: Find couple (A, S) such as (A, S) = argmin (A,S)≥0 V -AS 2 2
A classical gradient descent method may be used to solve Problem 1. Nevertheless, multiplicative update rules that minimize F have been developed [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. They offer a good compromise between speed and ease of implementation to solve Problem 1:

Theorem 1 The distance V -AS 2 is non increasing under the update rules:

S xp ← S xp (A t V ) xp (A t AS) xp A py ← A py (V S t ) py (ASS t ) py (5)
The proof of this theorem is given in Lee and Seung's publication [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF].

We precisely got interested in those update rules because of their ease of implementation and speed, for which they were initially created. Moreover, for initial non-negative matrices A and S, the constraint of non-negativity is inherent to the method: thanks to multiplicative rules, results remain nonnegative all along the iterations. The cost function definition, followed by the optimization step compose the classical NMF algorithm; an intrinsic initialization step chose non-negative matrices to start the algorithm [START_REF] Montcuquet | Nonnegative matrix factorization: a blind spectra separation method for in vivo fluorescent optical imaging[END_REF]. In our case, since a few information is available on fluorescence spectra, and to restrain the solution set, we chose to initialize matrix S with calibration spectra.

The NMF algorithm thus defined is described in this section:

Algorithm: NMF For i=1 to N iter 1) Initialize A and S with respectively positive constants and positive spectra models

2) Update S: S ← S (A t V ) (A t AS) 3) Update A: A ← A (V S t ) (ASS t )
4) restart steps 2 and 3 until stopping criterion is obtained (for example when r i-1 -r i 2 < , for r i = V -AS 2 and for a chosen ).

D. Constrained NMF

For our fluorescence imaging problem, since columns A p of A that refer to the fluorescent markers distribution in the medium are expected to be sparse, bring an a priori information on matrix A would help restrain the NMF solution set, and attenuate the non-uniqueness ambiguity. We implemented a sparse NMF algorithm that would help promoting sparsity of matrix A thanks to a thresholding step which at each iteration of algorithm after the optimization part that solves Problem 1 would favor sparse updates. Indeed both minimizing F and dealing with sparsity of columns A p of matrix A is usually simultaneously resolved applying a thresholding on data A p that brings a "sparsity promoting" behavior to the algorithm [START_REF] Bariteaux | An efficient numerical method for general regularization in fluorescence molecular tomography[END_REF].

Working out an appropriate threshold value at each iteration comes to find a number which can be comprised between 0 and the maximum of A p values. In order to relax the thresholding value choice, we prefer to select an initial sparsity value for concerned vector A p thus comprised between 0 and 1 as described in previous section; at each iteration our algorithm will find the best threshold value to eliminate background values and make sparsity of vector A p remain close from the initial chosen sparsity value.

The last iteration will end with a classical updating of matrices A and S (without thresholding step) in order to balance the fluorescence intensities between all fluorescence sources.

In this section, taking the classical algorithm as a basis, we propose to describe our sparse NMF algorithm implementation. The first NMF algorithm steps do not change from the classical one, but a sparsity step is added for each iteration in order to remain close to a wanted sparsity value ϕ defined by the user: Algorithm: NMF with sparsity constraints For i=1 to N iter 1) Initialize A and S with respectively positive constants and positive spectra models

2) Update S: S ← S (A t V ) (A t AS) 3) Update A: A ← A (V S t ) (ASS t )
4) Each column A p of A referring to weights of specific markers, for a wanted sparsity value ϕ for coefficients of column

A p , is changed into Ãp : ∀x ∈ (1, X), ãxp =    a xp if a xp ≥ max(a xp ) β pmin 0 otherwise (6)
with β pmin the threshold value such as:

β pmin = argmin βp sparsity Ãp (β p ) -ϕ 5
) restart steps 2 to 4 until stopping criterion is obtained (for example when r i-1 -r i 2 < , for r i = V -AS 2 and for a chosen ). 6) if i = N iter then run steps 2 and 3 before to end the algorithm.

As the previous algorithm is implemented with Matlab, we use the fminsearch function (Matlab) to solve the argmin problem. Fminsearch attempts to find a minimum of function sparsity Ãp (k) -ϕ where variable k is unknown, starting at an initial estimate of k. The threshold β pmin which is selected is the value from all β p tested that returns the minimum of the tested function.

Let us consider that the vector A p has, at first iteration of the algorithm, a sparsity value which is smaller (for example sparsity(A p )=0.2) than the chosen initial sparsity value ϕ (for example ϕ = 0.6). Then in step (4) of the algorithm, the fminsearch function will select an optimal threshold value β pmin from all β p tested: when this threshold is applied on A p data, all values in vector A p that are less than the threshold β pmin are set to zero, and new vector Ãp will be sparser than before. In that step of the algorithm, the threshold is chosen so that the value sparsity Ãp (k) -ϕ is minimized, but the threshold does not assure that the obtained sparsity value of A p is equal to ϕ.

A contrario if vector A p has a sparsity value greater (for example sparsity(A p )=0.7) than the chosen initial sparsity value ϕ (for example ϕ=0.6), then the obtained threshold β pmin at the end of step (4) of the algorithm will be equal to 0, since A p is already sparser than expected.

Then in both cases, steps 2 and 3 of the algorithm are restarted, and matrices A and S are updated thanks to the multiplicative update rules that minimize this time the cost function V -AS 2 . At that point, the sparsity value of vector A p may have changed since A has been updated, and sparsity value could thus exceed the chosen initial sparsity value ϕ. In that case, threshold obtained at further iteration would be set to 0, as explained previously. Finally by repeating steps 2 to 4 in turns, the cost function is minimized with sparse solutions being favoured, but priority is still given to the cost function minimization.

Sparsity value selection:

The constrained NMF algorithm asks for an initial sparsity value chosen by the user. When fluorescence sources distributions in media are sparse, looking for a sparse column A p is appropriate; this a priori information helps restraining the solutions set. We empirically noticed that a large choice of initial sparsity values, even far from the true expected sparsity value, help improving the classical NMF results. We propose in next section a set of simulation results that underline this idea.

For now no automatic selection of initial sparsity value has been implemented, and we empirically choose this parameter. NMF algorithm being quick to converge (a few seconds to one minute for a whole body mouse study), manual selection of this initial sparsity value does not prevent fast data processing.

IV. NMF VERSUS SPARSE NMF

In this section, we propose a comparative study between the NMF algorithm without spatial sparsity constraints and our new constrained algorithm that looks for spatially sparse solutions for the specific markers fluorescence.

A. Influence of sparsity value choice

A first comparison is made between both algorithms by simply running them on our simulated data (see Figure 4a) where sparsity of simulated vector A 1 is equal to 0.9017. Initialization for matrix S was the same for both tries, Gaussian were chosen, but lightly translated (50 nm) compared to true expected spectra. We test 11 initial sparsity values for the sparse NMF algorithm, from 0 to 1 to examine influence of sparsity initialization on results. Results of matrices A and S obtained after unmixing are presented Figure 4-b: the NMF algorithm without sparsity constraints did not manage to correctly unmix both sources. Taking sparsity constraints into account, here with an initial sparsity value equal to 0.8, led to an accurate result (see Figure 4-c). The imposed sparsity value of 0.8 finally reached 0.9027 at convergence of sparse algorithm when true sparsity value is 0.9017. We can also underline that a wrong initialization prevented the classical NMF algorithm to converge to the true spectra, but did not influence the sparse NMF algorithm. We propose in next part to study more carefully both algorithms' robustness to initialization.

Finally Figure 5 proposes sparsity values of A 1 and contrast between tumoral and healthy tissues obtained after unmixing for all initial sparsity values tested. We note on graph 5-a that thresholding steps of sparse NMF can only improve the classical NMF results. Only a few values of initial sparsity from 0 to 0.2 would not make any change on unmixing results compared to classical NMF performances (sparsity of unmixed vector A 1 is then equal to 0.2381). But for a large range of initial sparsity values, from 0.2 to 1, sparse NMF considerably improved classical NMF unmixing. In that range, initial sparsity values to bind to were at the origin of beneficial thresholding steps that helped the algorithm to converge to a sparse solution closer to the simulated A 1 .

This result is comforted by detection results presented on graph 5-b. Again for a range of initial sparsity values from 0.2 to 1, contrast C T,N between tumoral and healthy tissue areas is much closer to 1 than the contrast equal to 0.68 obtained with the classical NMF algorithm. The detection of the tumor is thus sensibly improved by the sparse NMF algorithm processing.

B. Comparison with other sparse NMF algorithms

As mentioned previously, several sparse NMF algorithms already exist; in this section we propose to test the algorithm developed by Hoyer [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF] on the simulated data already presented in the previous paragraph. Hoyers algorithm proposes a projected gradient descent algorithm for NMF with sparseness constraints that uses a projection operator to enforce a desired degree of sparseness. The main difference with our algorithm lies in that last property: as we propose a thresholding step to help the algorithm to converge to a more accurate solution, Hoyers method forces the solution of matrix A (or S in his article) to reach a fixed sparsity value. The sparsity value choice is thus crucial if we want to get a pertinent solution while the initial sparsity value we chose in our algorithm does not necessary has to be reached.

We implemented Hoyers algorithm with the possibility to impose different sparsity values for distinct columns of A as suggested by Stadlthanner in [START_REF] Stadlthanner | Extended sparse nonnegative matrix factorization[END_REF]. We ran it on previous simulated data (cf. Figure 4-a). We first studied sparsity of A 1 obtained after NMF, depending on the initial sparsity value chosen (see Figure 6-a) and compared it to our results: as expected, sparsity of vector A 1 is equal to the desired sparseness at convergence. Thresholding steps are more adapted to our problem since a larger range of initial sparsity values lead to pertinent results. This is especially useful when the size of the tumor can not be estimated and that the initial sparsity value has to be empirically chosen. On Figure 6-b we then compared detection results obtained with Hoyers algorithm or ours on that same example by plotting contrast C T,N obtained between tumoral and healthy areas after unmixing. Contrast obtained on that example (except for a unique initialization of the sparsity value at 0.1) is always greater than contrast found with Hoyers algorithm.

We compared both methods on simulated and in vivo data, and the thresholding method chosen in our algorithm is more robust to a larger choice of initial sparsity values than methods implemented by Hoyer and Stadlthanner.

C. Robustness to initialization

We want to test the robustness to initialization of our sparse NMF algorithm and to compare it with the classical NMF one. By making the initialized matrix S vary, the algorithm may lead to different solutions for the NMF decomposition. The solutions are indeed very sensitive to the initialization [START_REF] Montcuquet | Nonnegative matrix factorization: a blind spectra separation method for in vivo fluorescent optical imaging[END_REF]. On our simulation example, for a matrix S initialized exactly on the simulated spectra, the solutions obtained by both algorithms were the exact ones. We propose in this section to test the robustness to the initialization of both algorithms we compare, with or without sparsity constraints, on our simulated data. To create different initialization spectra for matrix S, we simply translate spectra we used for simulation (see Figure 7a) on a chosen wavelength range, from around 100 nm lower from true spectra emission peaks, to 100 nm upper: obtained range of N initialization spectra for matrix S is presented on Figure 7-b.

We then run both NMF algorithms, with or without sparsity constraints, N times for each initialization spectra tested. Obtained solution sets are presented Figure 8. Solution set obtained with the non-sparse algorithm presents various solutions depending on translated initializations (see Figure 8 -a) while the algorithm that takes sparsity constraints into account leads to a more robust solution set, even with a wide range of initializations.

The sparse NMF algorithm appears to be much more robust to the choice of the initialization than the classical NMF algorithm: the range of solutions for matrix S obtained when taking sparsity constraints into account gives the right positions for the fluorescence peaks when the classical algorithm presents wrong solutions. Only the initializations the closest from the true expected spectra allow the classical NMF algorithm to converge to right spectra. A contrario, the sparse NMF algorithm is robust enough to give results all almost superimposed to the true spectra. In next part, a more detailed study is made on solution set for both algorithms.

D. Range of admissible solutions to the NMF equation

As explained before, the NMF decomposition is not unique: without constraints, an infinity of solutions is admissible. Impose non-negativity to NMF results restrains the solution set.

To compare once more the sparse NMF algorithm and the classical one (without sparsity constraints), we propose to adapt a uniqueness study conducted by Moussaoui et al. in 2005 [32] to our data. This study defines -from a particular solution of NMF decomposition -the range of admissible solutions that lead to a same result. We briefly present the theory proposed by Moussaoui in the next section. We then propose to compare the ranges of admissible solutions on our simulated data obtained after NMF decompositions, with or without sparsity constraints.

1) Problem statement: Let us assume a factorization of V by the product of matrices A and S exists. If we now consider an invertible matrix T of size P ×P , then a new couple ( Ã, S) of solutions is easily found:

X = (AT -1 ) Ã (T S) S (7)
There is an infinity of factorizations of matrix V . The solution range may nevertheless be restricted by the non-negative constraint imposed to à and S.

2) Case for two sources: If two fluorescence sources are considered, for an initial set of data V of size m × n, matrices A and S respectively of size m × 2 and 2 × n are:

A = a 1 a 2 et S = s 1 s 2 (8)
Let us consider an invertible matrix T , of size 2 × 2:

T = 1 -α α β 1 -β (9)
To avoid indeterminacy problems between lines of S (and by default columns of A), the condition α + β < 1 is imposed. Invert matrix of T is given by:

T -1 = 1 1 -α -β 1 -β -α -β 1 -α (10)
Matrices à = A × T -1 and S = T × S may thus be expressed as:

∀k, S = (1 -α) × s 1k + α × s 2k β × s 1k + (1 -β) × s 2k ≥ 0 (11) ∀ , Ã = 1 1 -α -β × (1 -β) × a 1 -β × a 2 -α × a 1 + (1 -α) × a 2 ≥ 0 (12)
a) Admissible solutions range: Solve the preceding inequations for each term of matrices à and S leads to the admissible intervals in which α and β have to be included [START_REF] Moussaoui | Non-negative source separation: range of admissible solutions and conditions for the uniqueness of the solution[END_REF]; with B 1 = {k; s 1k < s 2k } and B 2 = {k; s 1k > s 2k } :

max k∈B1 s 1k s 1k -s 2k ≤ α ≤ min a 2 a 1 + a 2 (13) max k∈B2 -s 2k s 1k -s 2k ≤ β ≤ min a 1 a 1 + a 2 (14)
The ranges of admissible solutions for the NMF equation are constrained by the values of parameters α and β. b) Application to simulated data: We run both NMF algorithms (with or without sparsity constraints) on our simulated data. The initialization of S was perfectly chosen (that means equal to the expected spectra). We obtain a couple of solutions (A, S) for both algorithms. From that couple, ranges for parameters α and β and matrices T are calculated, and ranges of solutions for S are plotted. We compare ranges of admissible solutions obtained, taking or not sparsity constraints into account; results are presented Figure 9, initial solution S is plotted in black, and true solution S in dotted black. Sparsity constraints lead to a more restricted range of solutions, close to true solution expected, and gives once more better results than classical NMF algorithm. The signal of interest emitted by markers becomes quickly indistinguishable from autofluorescence (from around 20 mm deep). We run the classical NMF algorithm and our sparse NMF algorithm on the simulated data and calculate at each depth the contrast C T,N between tumoral and healthy areas after unmixing. The initial sparsity value for our sparse NMF algorithm was chosen equal to 0.7, away from the true simulated sparsity of vector A 1 (equal to 0.9017) on purpose. Results are presented Figure 11: when a contrast equal to 0.4 is reached for a tumor at 6 mm deep in tissues before unmixing (Figure 11-a), classical NMF processing allows to get a similar contrast for a tumor placed 12 mm deep (Figure 11-b). Finally even with no accurate sparsity initialization, contrast obtained with our sparse NMF algorithm is considerably improved compared to both prior cases (see Figure 11-c): tumor is detected with a contrast tumor/healthy tissue equal to 0.4 at 27 mm. Sparsity constraints and thresholding steps helped pushing back the detection limits.

2) Multiple sites tumors detection: Tumors may also grow at multiple sites due to the metastasis, and the robustness of our NMF algorithm for multiple tumors detection should be examined. We propose on our previous example to simulate a group of 5 tumors; different intensity levels are chosen in order to mimic fluorescence signals emitted by tumors at several depths in tissues (see Figure 12-a and 12-b). Each tumor is associated with the same fluorescence emission spectra S 1 as in previous examples. Depending of depth of tumor, spectra lightly vary (cf. [START_REF] Montcuquet | Nonnegative matrix factorization: a blind spectra separation method for in vivo fluorescent optical imaging[END_REF]). We run the sparse NMF algorithm with several initial sparsity values from 0 to 1: results of sparsity of vector A 1 obtained after unmixing are presented Figure 13. Once more, for a large range of sparsity values, from 0.25 to 1, unmixing results are considerably improved compared to classical NMF results (see green line on Figure 13) even for multiple sites tumors. For smaller initial sparsity values (0 to 0.25), obtained results are similar to the classical NMF ones, and the thresholding step had no effect on data. Sparsity constraints could only improve NMF results but would not return less accurate unmixing results.

We propose in next section to test the sparse NMF algorithm on mice data.

V. in vivo UNMIXING RESULTS

To test our algorithm this time on real data, an in vivo experiment is performed on a mouse. In vivo experiments imply that an autofluorescence signal is necessarily measured. We want to test the NMF algorithm on three sources unmixing: we thus chose two specific markers to separate from each other, and from autofluorescence signal.

A. Feasibility experiment

The animal procedure was in compliance with the guidelines of the European Union (regulation n • 86/609), taken in the French law (decree 87/848) regulating animal experimentation. All efforts were made to minimize animal suffering. The animal manipulation was performed with sterile techniques and approved by the Grenoble Animal Care and Use committee (France) (registration number 20 iRTSV Léti-FNG-02). An adult female nude mouse (Janvier, Le Genest saint-isle, France) was used throughout the experiments. It was housed in approved facilities, at 21±1 • C under diurnal lighting conditions. The mouse arrived at the animal facility two weeks before the experiments start and had free access to food and water.

The animal is placed on a translation stage with N y positions of travel range. To acquire spectrally resolved measurements and a whole scanning, the mouse is illuminated with a laser line at 690 nm, and the stage is translated N y times. For each position of the stage, the emitted back fluorescence signal is collected along a line of N x points by an imaging spectrometer coupled with a charge-coupled device camera (Andor Technologies): a N x × N λ acquisition is measured (see Figure 14). At the end of the scan, N y acquisitions are obtained, each of size N x × N λ . For this experiment, N x was equal to 255 and N λ to 1024, which corresponds to a wavelength range around 590 to 973 nm. Before to run the NMF algorithm, mixed data of size (N y × N x × N λ ) are reordered as a 2-D array of size (N y × N x , N λ ).

To be able to quantify our results, especially the unmixed autofluorescence signal, a first acquisition is performed on the mouse, with no specific fluorescent markers: only the autofluorescence of the mouse tissues is detected. The resulting acquisition is presented Figure 15-a. Besides are acquired ex vivo spectra of both specific markers we use: Alexa 750 and Indocyanine Green loaded into Lipid Nanoparticules (ICG-LNP) [START_REF] Navarro | A novel indocyanine green nanoparticle probe for non invasive fluorescence imaging in vivo[END_REF] (see Figure 16).

To simulate marked tumors after circulating tumor-specific markers have been eliminated everywhere except in the tumor, two glass capillary filled with specific markers are inserted subcutaneously (see Figure 15 autofluorescence, ICG-LNP and Alexa 750 -have overlapping emission spectra.

Once the capillary tubes placed, the fluorescence acquisition of the animal is run. Obtained result is presented Figure 15-c.

B. Results

We ran both algorithms we want to compare on our experimental data. For the sparse NMF algorithm, a sparsity value equal to 0.8 has been empirically chosen. The different unmixing results are presented Figure 17. Both algorithms succeeded to separate the three fluorescence sources: in the three fluorescence images obtained for each algorithm, the highest intensities correspond to true position of markers (see unmixed fluorescence contributions on Figure 17).

Nevertheless, sparse NMF algorithm gave more accurate results. This can first be noticed by comparison with the fluorescence contributions obtained with the classical algorithm, particularly on ICG-LNP and autofluorescence ones: ICG-LNP contribution obtained with classical algorithm (see image 2 on 17 -a) contains background parts that were not accurately separated. Those background parts appear as missing on the autofluorescence contribution (see image 3 on Figure 17 a, comparison can be made with the initial mixed data). Autofluorescence contribution can also directly be compared to the test acquisition without specific fluorescence (Figure 15-a). Finally, all the drawbacks underlined before, as wrong background remaining and not accurate separation, do not appear -or lesser -in the results obtained with the sparse NMF algorithm (see Figure 17-b, images 1,2 and 3).

Resulting A matrices translate all the observations previously made on images 1, 2 and 3. Indeed, the fluorescent markers weights (columns of A) are expected to be spatially sparse to agree with the capillary tubes positions. Such result is reached with the sparse NMF algorithm (see weights of matrix A on Figure 17 -b)), while an incorrect background remains in columns of A for specific markers obtained with the classical NMF algorithm (Figure 17 -a)), as already noticed on images 1, 2 and 3 of Figure 17-a.

Moreover, we can look more in details at spectra obtained in matrix S with both algorithms: compared to emission spectra from our ex vivo and autofluorescence measurements (see Figure 16), specific fluorescence spectra obtained after classical NMF separation are moving further away from expected fluorescence spectra (see Figure 17-a, matrix S) than those obtained with the sparse NMF algorithm (see Figure 17-b, matrix S), that much closer to reality.

VI. CONCLUSION

The specific fluorescence signal of specific markers used in optical imaging allows to detect targets, as tumors. Even if we use infra-red light, the autofluorescence of biological tissues is also detected in this wavelength range: it needs to be removed to get accurate detection results. We introduced the Non-negative Matrix Factorization as blind source separation method to unmix fluorescence spectra and eradicate autofluorescence. In order to improve detection of spatially sparse fluorescent markers, we suggested a new NMF algorithm with sparsity constrained. This algorithm asks for a wished sparsity value to reach: we propose a simplified way to approximate sparseness from the size of tumors we are looking for.

A study has been conducted on simulated data to compare sparse and non-sparse algorithms unmixing performances. Our sparse NMF algorithm appeared to be more robust to the random choice of initialization spectra, and returned results closer to the expected spectra than the non-sparse algorithm. NMF method without constraints offering infinity of solutions, we got interested in the range of admissible solutions for both algorithms: once more, sparse NMF algorithm gives better results than the algorithm without sparsity constraints.

Finally, both algorithms were run on in vivo acquisitions, and successfully unmixed up to three different fluorescent sources. Acquisitions had been obtained on a mouse with subcutaneous specific fluorescent markers. Once more, results are closer to reality with the sparse NMF algorithm than with the non-constrained algorithm, which comforts simulation outcomes.

Nevertheless, even if fluorescent markers are expected to be specific enough to only accumulate around the tumors, reality of in vivo experiments is for now different. Indeed, once IV-injected, currently used markers travel in blood and lymphatic canals and spread everywhere in tissues: one may talk of non-specific signal. Only a few supplementary amount of markers accumulate around the tumor: the signal measured there is called specific signal. Even if more markers will accumulate around the tumor, a small amount of non-specific signal will always be detected everywhere in tissues. In such case, sparsity constraints may seem inappropriate for the unmixing problem. However, injected specific markers are developed to be more and more specific to tumors. We thus may expect in the near future that markers will be specific enough to avoid a non-specific signal detection, and use sparsity constraints for the unmixing step.

As optical imaging tries to detect deeper and deeper embedded tumors, NMF -as a preprocessing tool to remove autofluorescence signal and isolate specific fluorescence contributions, sharpened by several constraints -is an helpful tool. Sparse NMF algorithm we proposed greatly improves the preprocessing results in spatially sparse makers localization.

  Marker-related column (a) Inaccurate unmixing (b) Correct unmixing

Fig. 1 .

 1 Fig. 1.Link between quality of unmixing result and weight matrix A obtained.

Fig. 2 .

 2 Fig. 2. Sparsity value ranges from 0 to 1, from non-sparse to highly sparse signals: example on weights vectors.

Fig. 3 .

 3 Fig. 3. By summing of a specific fluorescence signal a) with an associate fluorescence spectrum and an autofluorescence signal b) with a different fluorescence spectrum, we obtain simulated mixed data c).

Fig. 4 .Fig. 5 .

 45 Fig. 4. First comparison of sparse and non-sparse NMF algorithms: a) simulated data, b) Results obtained with the classical NMF algorithm, c) Results with our sparse NMF algorithm (for an initial sparsity value equal to 0.8). Taking sparsity constraints into account leads to more accurate results.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Comparison between our sparse NMF algorithm and Hoyers. a) Sparsity values of A 1 , and b) contrast C T,N between tumoral and healthy tissues obtained after unmixing for 11 initial sparsity values for A 1 tested between 0 and 1 for both algorithms: Hoyers and ours

Fig. 9 .

 9 Fig. 9. (a) Range of admissible solutions with the classical NMF algorithm, without sparsity constraints (b) Range of admissible solutions, with the constrained NMF algorithm (sparsity constraints).

  On a specific example similar to previous ones, we simulate the intensity decrease of the fluorescence signal emitted by markers when they are moved in depth from 1 mm to 40 mm in tissues: results are presented Figure10.

Fig. 10 .

 10 Fig. 10. A simulated tumor binded by fluorescent markers are moved in depth in breast tissues, from 1 mm to 40 mm. Specific signal emitted by tumor is quickly indistinguishable from autofluorescence signal (results of mixed signals from 20 mm to 40 mm are the same).

Fig. 11 .

 11 Fig. 11. Deep tumors detection results: a) contrast obtained before unmixing (contrast equal to 0.4 at 6 mm), b) contrast after classical NMF processing (contrast equal to 0.4 at 12 mm), c) contrast obtained with our sparse NMF algorithm (contrast equal to 0.4 at 27 mm).

Fig. 12 .

 12 Fig. 12. Simulation of multiple sites tumors. a) Mixed data, b) 5 tumors to detect.

Fig. 13 .

 13 Fig. 13. Sparsity values of A 1 obtained after sparse NMF unmixing, and associate tumors weights results (blue line). Comparison with true simulated sparsity of A 1 (red line) and classical NMF results (green line).

Fig. 14 .Fig. 15 Fig. 16 .

 141516 Fig. 14. Experimental set-up: acquisition on the animal and data processing.Mixed data

Fig. 17 .

 17 Fig. 17. Results of unmixing obtained with both algorithms: separated fluorescence contributions, and matrices A and S at convergence: a) with classical NMF algorithm, without sparsity constraints and b) with the regularized NMF algorithm, when taking sparsity constraints into account.

CEA, LETI, MINATEC, 17 rue des martyrs, 38054 Grenoble cedex 9, France.

GIPSA-Lab, Rue de la houille blanche BP46,

Saint Martin d'Hères, France.


Longyu Jiang
email: long-yu.jiang@gipsa-lab.inpg.fr

Florian Aulanier

Grégoire Le Touzé

Barbara Nicolas

Jérôme I Mars

Raypath Separation with High Resolution Processing

In this paper, a high resolution spatio-temporal analysis algorithm is presented and applied to separate the raypaths recorded with sensors for ocean acoustic tomography. It is a extended version of classical MUSIC, we propose it to an active large band processing dedicated for ocean acoustic tomography.We use spatial and frequency smoothing method to built projector. The experiment is done on simulated datas.

with classical technology using in ocean acoustic tomography, Beamforming, MUSICAL obtains more precise separation results. Besides, the results obtained by MUSICAL are with less artifacts.

I. INTRODUCTION

In the ocean, it is possible to get images of the sound speed variation thanks to acoustic tomography. A common way to proceed is to take advantage of the multi-path properties of the wavefield [START_REF] Munk | Cambridge monographs on mechanics[END_REF]. Each path provides information on the variation of sound speed distribution. Recently, a tomography method has been developed using vertical arrays in emission and reception and double beamforming to separate the different paths and extract more observables. Thanks to 2D processing, the method hugely improve the conventional beamforming performances [START_REF] Iturbe | Ocean acoustic tomography using a double-beamforming algorithm[END_REF]. However, it is still confronted to the main beamforming drawback : the low resolution performances. As a first step of improvement, we apply a high resolution processing which combine MUSICAL algorithm [START_REF] Gounon | High resolution spatio-temporal analysis by an active array[END_REF], spatial and frequency smoothing [START_REF] Paulus | New multicomponent filters for geophysical data processing[END_REF]. We compare its resolution performances to those of conventional beamforming processing on simulated data.

II. HIGH RESOLUTION SPATIO-TEMPORAL ANALYSIS

High resolution methods are based on an accurate model of interference wave fields (the received signals). Then certain parameters are estimated by maximizing or minimizing a function [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]. Based on [START_REF] Gounon | High resolution spatio-temporal analysis by an active array[END_REF], we present a high resolution spatiotemporal analysis method to estimate jointly drection and arrival time, by taking into account the a priori information of the signal in MUSIC with a active antenna. In the following part, we give details.

A. Signal model

The model is built on an acoustic field composed of P raypaths on a active vertical antenna of M sensors.The temporal signal received on the m th sensor is modelized as:

x m (t) = P p=1 a p e p (t -τ m,p ) + b m (t) (1)
In frequency domain, ( 1) is written as:

x m (ν) = P p=1 a p e p (ν) exp(-j2πτ m,p ) + b m (ν) (2)
The arrival time τ m,p can be expressed as follows:

τ m,p = T p + t m (θ p ) (3)
T p represents the delay between the p th source and the reference sensor, t m (θ p ) is the delay between the reference sensor and the m th sensor. t m (θ p ) is function of θ p , which is the direction of raypaths on the antenna.

In [START_REF] Iturbe | Ocean acoustic tomography using a double-beamforming algorithm[END_REF], the term e p (ν) is deterministic. It is composed of a modulus and a phase expressing a priori information about the signal. And a p 's are considered random and decorrelated. In ocean acoustic tomography, each source corresponds to a path, e(t) is the emitted signal and a p corresponds to the complex amplitude of each path.

B. Estimation of interspectral matrix

Spatial and frequency smoothing can be introduced respectively to estimate the interspectral matrix. Making a strong smoothing either distance or frequency may introduce significant bias in the estimation of the matrix. Thus, Paulus [START_REF] Paulus | New multicomponent filters for geophysical data processing[END_REF]combines the two types of smoothing methods together. ( K s K f respectively denotes the spatial and frequency smoothing factor).

From the single available observation x, it is possible to generate a set of (2K s + 1) spatially offset recurrences x ks . These 2K s + 1 recurrences may themselves frequentially be shifted to obtain K = (2K s + 1)(2K f + 1)recurrences x ks,k f . Then we can then estimate the wideband interspectralmatrix by the following formula:

Γ = (2K s + 1)(2K f + 1) -1 2Ks+1 ks=1 2K f +1 k f =1 x ks,k f x * ks,k f (4)

C. Principle of the Algorithm

By eigendecomposition, the signal subspace is spanned by the first P eigenvectors of Γ, and its complementary, the noise subspace is spanned by the orthogonal M F -P last eigenvectors. The orthogonal projection onto the noise subspace is estimated as:

Π B = M F k=P +1 µ k µ k * (5)
Finally, the HR algorithm consists in maximizing the following fuction:

F (θ, T ) = (a(θ, T ) * Π B a(θ, T )) -1 (6)
With: a(θ, ): wide-band steering vector and is written as: a(θ, T ) = [e(ν 1 )e -2iπν1T d + (ν 1 , θ), . . . , e(ν F )e -2iπνF T d + (ν F , θ)] + d(ν i , θ) is the conventinal steering vector used in narrow band analysis [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF].

III. SIMULATIONS

We illustrate the method presented previously on a set of simulated data from a software propagation parabolic equations. In this experiment, we use one source and a vertical antenna of 61 receivers. The source is fixed at 40.5 m under the ocean. The 61 receivers are regularly spaced in the water column between 30 m and 60 m. We choose the 31th sensors as the reference one. The distance between the source and the reference sensor is 2 km.

Separation results are demonstrated by Fig. 1 and Fig. 2. Globally, for each arrival raypaths, estimation values obtained by MUSICAL vary in a smaller scale and are with less artifact. Specially, shown by Fig. 3 and Fig. 4, differring to Beamforming, the first two arrival raypaths are separated totally.

Fig. 1 .Fig. 3 .Fig. 4 .

 134 Fig. 1. Separation results with Beamforming


Anne-Sophie Montcuquet

Lionel Hervé

Fabrice Navarro

Jean-Marc Dinten

Jérôme I Mars

NON-NEGATIVE MATRIX FACTORIZATION UNDER SPARSITY CONSTRAINTS TO UNMIX IN VIVO SPECTRALLY RESOLVED ACQUISITIONS

Keywords: Fluorescence imaging, positive source separation, fluorescence in vivo spectroscopy

Fluorescence imaging in diffusive media is an emerging imaging modality for medical applications which uses injected fluorescent markers (several ones may be simultaneously injected) that bind to specific targets, as tumors. The region of interest is illuminated with near infrared light and the emitted back fluorescence is analyzed to localize the fluorescence sources. To investigate thick medium, as the fluorescence signal decreases with the light travel distance, any disturbing signal, such as biological tissues intrinsic fluorescence -called autofluorescence -, is a limiting factor. To remove autofluorescence and isolate each specific fluorescent signal from the others, a spectroscopic approach, based on Non-negative Matrix Factorization, is explored. We ran an NMF algorithm with sparsity constraints on experimental data, and successfully obtained separated in vivo fluorescence spectra.

INTRODUCTION

In fluorescence imaging, specific markers are injected to a patient, and bind specifically to targeted compounds, like tumors [START_REF] Ntziachristos | Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin v-cy5. 5 conjugate[END_REF]. Several specific markers may be injected at once, and bind to different compounds or organs. The region of interest is illuminated with near infrared (NIR) light; an optimal wavelength range may be defined -between 600 and 900 nm -where the tissue absorption is lower and the marker optimally excited. Finally, the emitted back fluorescence signal is measured and the fluorescent source is localized. To investigate thick media for medical diagnostic application ( 4 cm), as the fluorescence signal gets exponentially weak with the light travel distance, the autofluorescence of tissues becomes a limiting factor. The analysis of a fluorescence signal impaired by autofluorescence may lead to a wrong localization of the markers: the signal needs to be preprocessed in order to remove autofluorescence. For in vivo fluorescence spectroscopy, the unmixing problem is referred to as a blind source separation problem since the spectra may vary according to the fluorescent dye biological environment. Fluores-cence spectra to separate are also supposed statistically dependent, which filters out many methods (such as ICA). Nonnegative Matrix factorization (NMF) differs from classical source separation methods (SVD, PCA) in that it forces the matrices factors to be non-negative, and thus suits to spectra separation. We propose to test this method on spectroscopic data, and to define a new regularized NMF algorithm constrained by sparsity.

NON-NEGATIVE MATRIX FACTORIZATION

For initial non-negative mixed data M , NMF proposes to find a couple of matrices (A, S) with non-negatives coefficients, whose product optimally approaches M . The classical NMF definition says [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]:

Given a non-negative matrix M ∈ R X×Y , find non-negative matrices A ∈ R X×P and S ∈ R P ×Y such that :

M AS (1)
where non-negative matrices are matrices whose all factors are nonnegative and P stands for the number of sources to unmix.

Applied to spectroscopy, matrix A is considered as the weight vector, and S as the spectra matrix, both containing as much elements P as fluorescent sources to separate.

To find the particular matrices A and S that satisfy equation 1, two distinct steps must be considered. First, we need to define a criterion that link M with A and S. Then, in a second step, optimizing this criterion under the non-negativity constraint would lead to the best solution possible. Several criteria (the Euclidean distance, the Kullback-Leibler divergence [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]...) and different optimization methods (Alternating Least Square, update rules) may be chosen to obtain an NMF decomposition.

SPARSITY CONSTRAINTS

For local specific markers, we expect to get peaked weight columns in matrix A. A contrario, when unmixing is not complete, a residual autofluorescence baseline surrounds the specific marker peak. This property is depicted figure 1. An intuitive reasoning to improve unmixing would be to smooth over the unwanted autofluorescence residuals on specific markers weight profiles, by thresholding the smallest values: in other words, we look for sparse A columns for specific markers contributions. We introduce a sparsity value [START_REF] Hoyer | Non-negative matrix factorization with sparseness constraints[END_REF]: let us consider a weight matrix A of size X × P , the sparsity of a given column A p (p ∈ (1, P )) of A is:

sparsity (A p ) = √ X - X x=1 |a xp | X x=1 a 2 xp √ X -1 (2)
Sparsity value ranges from 0 for non-sparse results to 1 for extremely sparse results, as depicted figure 2.

Sparsity= 0.3 Sparsity= 0.5 Sparsity= 0.8 Sparsity= 0.9 Fig. 2. Sparsity value ranges from 0 to 1, from non-sparse to highly sparse signals.

In the next section, we propose an NMF algorithm that finds the perfect threshold so that our matrix A reaches a chosen sparsity value.

SPARSITY-CONSTRAINED NMF ALGORITHM

Classical algorithm

We first present a classical NMF algorithm without any constraint. In 2001, NMF popularity increased after Lee and Seung published two new NMF algorithm, based on the use of multiplicative update rules that minimize specific criteria.

Cost function definition

Different criteria to minimize, or cost functions, can be used [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]: we may cite the square of the euclidean distance between M and AS, and the Kullback-Leibler divergence criterion. Here, we get interested in the cost function F , chosen as the square of the Euclidean distance between M and AS [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], lower bounded by 0, and defined as:

F = Nx x=1 Ny y=1 (m xy - P p=1 a xp s py ) 2 = M -AS 2 (3)

Optimization

The following optimization problem is thus considered:

Problem 1 Minimize F with respect to A and S subject to the constraints A, S ≥ 0.

Lee and Seung proposed multiplicative update rules to solve Problem 1: those iterative rules offer a good compromise between speed and ease of implementation to solve the optimization problem (see problem 1); they are presented in the form of the following theorem:

Theorem 1 The distance M -AS 2 is non increasing under the update rules:

S xp ← S xp (A t M ) xp (A t AS) xp A py ← A py (M S t ) py (ASS t ) py (4)
The proof of this theorem is given in Lee and Seung's publication [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF].

We precisely got interested in those update rules because of their ease of implementation and speed, for which they were initially created.

Constrained NMF

In this section, taking the classical algorithm as a basis, we propose a new algorithm that takes into account sparsity constraints. The first NMF algorithm steps do not change from the classical one, but several sparsity steps are added for each iteration in order to reach a wanted sparsity value:

Algorithm: NMF with sparsity constraints 1. Initialize A and S with respectively positive constants and positive spectra models 2. Update S: S ← S (A t M )

(A t AS) 3. Update A: A ← A (MS t ) (ASS t )
4. Each column A p of A referring to weights of specific markers, for a wanted sparsity value spa for coefficients of column A p , is changed into Ãp :

∀x ∈ (1, X), ãxp =      a xp if a xp ≥ max(a xp ) β p 0 otherwise (5)
with β p such as:

β p = argmin βp sparsity Ãp (β p ) -spa
5. restart steps 2 to 4 until stopping criterion is obtained (for example when M -AS 2 = for a chosen ).

RESULTS

For in vivo experiments, an autofluorescence signal is necessarily measured. Then several specific markers may be used to simulate marked targets, such as tumors. In this section, we test NMF to unmix three overlapping different fluorescence sources, including the autofluorescence.

Methods

We perform feasibility experiments on a mouse. The animal procedure was in compliance with the guidelines of the European Union (regulation n86/609), taken in the French law (decree 87/848) regulating animal experimentation. All efforts were made to minimize animal suffering. The animal manipulation was performed with sterile techniques and approved by the Grenoble Animal Care and Use committee (France) (registration number 20 iRTSV Léti-FNG-02).

An adult female nude mouse (Janvier, Le Genest saint-isle, France) was used throughout the experiments. It was housed in approved facilities, at 21±1C under diurnal lighting conditions. The mouse arrived at the animal facility two weeks before the experiments start and had free access to food and water.

To acquire spectrally resolved measurements, the animal is illuminated with a planar laser at 690 nm. The emitted back fluorescence signal is collected along a line of N x points by a spectrometer coupled with a charge-coupled device camera (Andor Technologies): a N x × N λ acquisition is measured (see Figure 3). A translation stage, covering N y steps, is then used to get a scanning of the whole animal: N y acquisitions are obtained, each of size N x × N λ . Before to run the NMF algorithm, our data of size (N y × N x × N λ ) are reordered as a 2-D array of size (N y × N x , N λ ).

Fig. 3. Experimental set-up: acquisition on the animal and data processing.

To draw a parallel between acquisitions with or without specific fluorescence, a first acquisition of the animal is run as a reference, without any specific markers inserted (Figure 4-a). Then, several specific markers are added to simulate marked targets. First, a unique glass capillary tube filled with a specific marker is inserted subcutaneously on the animal, to classically simulate a marked target. The chosen marker is here 5 µl of Indocyanine Green loaded into Lipid Nanoparticules (ICG-LNP) [START_REF] Navarro | A novel indocyanine green nanoparticle probe for non invasive fluorescence imaging in vivo[END_REF] at 0.35 µmol/l. We will use this first example to illustrate the impact of sparsity constraints, running the algorithm with different wanted sparsity values. Subsequently, a second capillary tube filled with 5 µl of Alexa 750 (an other specific marker) at 0.1 µmol/l is inserted. Three distinct fluorescent sources -autofluorescence, ICG-LNP and Alexa 750 -whose emission spectra are overlapping have now to be unmixed. The resulting intensity scanning are presented Figure 4-b and Figure 4-c. In this precise example, we chose to constrain the sparsity level to 0.9 for the specific markers -ICG-LNP and Alexa 750 -weight matrices.

Results

We run the NMF algorithm on our data. For the first experiment, with a unique capillary tube (see data Figure 4-b), we analyze the sparsity constraint effect on results, by increasing the chosen value spa(see 4.2). The number of components equal to zero in the column corresponding to ICG-LNP in matrix A increases with sparsity value. Moreover, the spectra of matrix S get closer to initializations with sparsity value increasing: S was here initialized on purpose with mean spectra of considered data. Spectra of matrix S getting closer to initialization may mean that results improve. For now, it is just an empirical observation to confirm. Finally, as expected, a strong sparsity value allows to remove the autofluorescence residuals and improves results (see Figure 5). Figure 6 presents unmixing results of the second experiment (with two capillary tubes, initial data presented on Figure 4-c), for an high value of sparsity (spa = 0.9). NMF under sparsity constraints successfully unmixed the three components. Nevertheless, a few errors may be noticed on the autofluorescence contributions: holes are present on ICG-LNP and Alexa 750 signals locations.

CONCLUSION AND FUTURE WORK

Beyond the specific fluorescence signal of specific markers used in optical imaging, the autofluorescence of biological tissues is also detected in the wavelength range we use, and needs to be removed to get accurate detection results. To un- mix fluorescence spectra, a blind source separation method was chosen: the Non-negative Matrix Factorization. We presented a classical NFM algorithm, before to introduce sparsity constraints. The constrained NMF algorithm was run on in vivo acquisitions, and successfully unmixed up to three different fluorescent sources. Empirically, results seem to improve with sparsity value imposed to specific markers weights profiles, but we still have to prove the convergence of the constrained algorithm. Moreover, sparsity value is by now empirically chosen: future work will explore methods to predefine the best sparsity value to use, and improve the existent algorithm. As optical imaging tries to detect deeper and deeper embedded tumors, NMF -as a preprocessing tool to remove autofluorescence signal and isolate specific fluorescence contributions, sharpened by several constraints -is an helpful tool.

Fig. 1 .

 1 Fig. 1. Link between the quality of unmixing process and resulting weight matrix A.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Autofluorescence acquisition (b) ICG-LNP tube inserted (c) ICG-LNP and Alexa 750 tubes inserted.

Fig. 6 .

 6 Fig. 6. Two capillary tubes experiment: unmixing results