In Vitro and In Vivo Multispectral Photoacoustic Imaging for the Evaluation of Chromophore Concentration

Abstract

Multispectral photoacoustic imaging is a powerful noninvasive medical imaging technique that provides access to functional information. In this study, a set of methods is proposed and validated, with experimental multispectral photoacoustic images used to estimate the concentration of chromophores. The unmixing techniques used in this paper consist of two steps: (1) automatic extraction of the reference spectrum of each pure chromophore; and (2) abundance calculation of each pure chromophore from the estimated reference spectra. The compared strategies bring positivity and sum-to-one constraints, from the hyperspectral remote sensing field to multispectral photoacoustic, to evaluate chromophore concentration. Particularly, the study extracts the endmembers and compares the algorithms from the hyperspectral remote sensing domain and a dedicated algorithm for segmentation of multispectral photoacoustic data to this end. First, these strategies are tested with dilution and mixing of chromophores on colored 4% agar phantom data. Then, some preliminary in vivo experiments are performed. These consist of estimations of the oxygen saturation rate (sO${}_2$) in mouse tumors. This article proposes then a proof-of-concept of the interest to bring hyperspectral remote sensing algorithms to multispectral photoacoustic imaging for the estimation of chromophore concentration.

1. Introduction

Photoacoustic imaging is a hybrid medical imaging technique that provides access to functional information. Nanosecond laser pulse illumination can provide optical properties of biological tissues that are related to molecular composition. The transmitted optical energy is absorbed by the optical absorbers in the tissues, which then undergo thermal expansion. At each laser pulse, this phenomenon generates ultrasonic waves that can be detected at the tissue surface by a conventional ultrasound system. The use of several laser wavelengths provides the multispectral absorption characteristics that can be used to distinguish the imaged tissues from each other [1,2]. In addition, as the absorption by a chromophore is linearly related to its concentration, multispectral photoacoustic imaging can also be used to determine the concentrations of chromophores present in the imaged region [3,4].

These two properties, discrimination of tissues and concentration evaluation, are of major interest, depending on the application at hand. They can be exploited to determine different concentrations of the same chromophore (e.g., estimation of the concentration of a contrast agent in the body [5]), or to distinguish one particular chromophore from all of the other imaged ones without considering its dilution (e.g., determination of the level of vascularization and calculation of the concentration of oxygenated and deoxygenated blood, for oxygenation rate evaluation [6,7], cancer tissue evaluation [8], or imaging connectivity in brain [9]).

In photoacoustic images, each pixel can be interpreted as a mixture of pure chromophores, where the weights represent their average concentrations. The spectra of pure chromophores are called endmembers, and their weights are the abundance coefficients. Unmixing techniques for multispectral photoacoustic data are presented in this article, with the objective being to estimate the concentrations of the chromophores present in the imaged area. In the field of photoacoustic imaging, several unmixing methods have already been used to detect chromophores from the body [10] or contrast agents injected into tissues. For instance, principal component analysis (PCA) and independent component analysis were used in Reference [3]. However, the abundance matrices calculated by these two methods do not allow accurate estimation of concentrations of chromophores. In particular, these matrices can contain negative values that cannot be interpreted as concentrations even after scaling. This issue is one of the challenges addressed in this paper. Indeed, this study highlights experimental results as a proof of concept of the usefulness of photoacoustic imaging to evaluate chromophore concentration.

The concentration of an endmember should be equal to zero if the chromophore is absent in the considered pixel, and equal to one if it is a pure pixel. The endmember will lie in the interval $(0,1)$ if the pixel is partly composed of the chromophore (i.e., diluted chromophore, or a mix of chromophores). The abundance matrix must then be subject to the following constraints: (1) each abundance coefficient must be in the range $[0,1]$; and (2) for a given pixel, the sum of all K abundance coefficients, which correspond to the K endmembers, must be less than 1 (i.e., diluted chromophore) or equal to 1 (i.e., pure chromophore, or a mixture of chromophores). These constraints have never been investigated in the field of photoacoustics, while they are usually encountered in remote sensing, where imaged areas might contain, for example, roads, buildings, rivers, or forests, which correspond to the chromophores to be unmixed. The abundances are, thus, related to the local concentrations of the chromophores in the scene. Spectral unmixing is also a widely explored field for various applications, like food safety [11,12], pharmaceutical process monitoring [13], and industrial and forensic applications [14,15]. However, unmixing of remote sensing data is the most closely related area to our multispectral photoacoustic data. As many efficient algorithms have been devised and analyzed over the last decade, we focus on this particular field in this study. The main challenge of this paper is then to show that photoacoustic imaging can benefit from these innovations, principally considering the constraints presented above, because of the similarities between these two fields. However, photoacoustic imaging contains less spectral bands than remote sensing and is subject to spectrally dependent light attenuation which makes the unmixing even more challenging. Moreover, the methods presented in this study are based on relative comparison between spectra. Even if PCA is also a relative method, as mentioned before, it does not allow proper evaluation of concentration, based on authors preliminary experiments. The study explores, then, other field to adapt unmixing methods that better suit the considered goal. As all spectra are impacted in the same way by the light absorption and ultrasound attenuation for each wavelength, and these methods allow us to avoid taking into account these parameters.

2. Unmixing of Photoacoustic Data

A multispectral photoacoustic image consists of a two-dimensional area that is imaged at L wavelengths. Each pixel ${\mathbf{x}}_i$ is characterized by its coordinates ${\mathrm{s}}_i\in \mathrm{I}{\mathrm{N}}^2$ and is endowed with a spectrum ${\mathit{A}}_i\in \mathrm{I}{\mathrm{R}}^L$. The N samples, ${\mathbf{x}}_i$, of the region of interest are expressed as follows:

$${\mathbf{x}}_i=\left[{\mathrm{s}}_i,{\mathit{A}}_i\right]\in \mathbf{X},$$

(1)

with $i\in \{1,\dots ,N\}$ as the sample index; see Figure 1a. Before any other processing, the data need to be normalized over the range $[0,255]$. This allows the same range of values for data from different acquisition systems. Two types of pixels are illustrated in Figure 1: pixel x${}_1$ contains only noise, called the background, with low and flat photoacoustic amplitudes at all wavelengths (Figure 1b, green curve A${}_1$), and pixel x${}_2$ provides information on the optical absorbers, with a high photoacoustic signal (Figure 1b, purple curve A${}_2$).

2.1. Linear Mixing Model

The acquisition of photoacoustic data over a region of interest at several wavelengths provides a multispectral (i.e., three-dimensional) image, where each pixel is endowed with a spectrum. This spectrum is either pure, and is then considered as a single endmember, or mixed if it is composed of a mixture of endmembers. The linear mixing model (LMM) [16] considers that each mixed pixel is a convex combination of the spectra of the endmembers. More formally, this is defined as follows:

$$\begin{array}{cc}\hfill {\mathit{A}}_i& ={\displaystyle \sum _{k=1}^K}{u}_{ki}{\mathit{E}}_k+{\mathit{g}}_i,\forall i\in \{1,\dots ,N\},\hfill \end{array}$$

(2)

where ${\mathit{A}}_i\in \mathrm{I}{\mathrm{R}}^L$ is the spectrum at the i-th pixel, L is the number of wavelengths, K is the number of endmembers, $u_{ki}$ is the abundance of the k-th endmember in the i-th pixel, ${\mathit{E}}_k$ is the L-dimensional spectrum of the k-th endmember, and ${\mathit{g}}_i$ is a vector of zero-mean white noise that defines the sensor noise and the error of the model. All of the vectors are column vectors. Then, $\mathit{U}$ is used for the matrix of abundances with the generic $(i,j)$-th entry $u_{ij}$. As abundances are relative contributions, they must be positive and their sum has to be equal to one:

$$\left\{\begin{array}{c}{u}_{ki}\ge 0,\forall i\in \{1,\dots ,N\},\forall k\in \{1,\dots ,K\}\hfill \\ {\displaystyle \sum _{k=1}^K}{u}_{ki}=1,\forall i\in \{1,\dots ,N\}.\hfill \end{array}\right.$$

(3)

The LMM is a simple but effective model that is widely used in remote sensing [16]. For medical applications, the LMM can be a powerful tool for estimation of concentration of chromophores. In this context, it is assumed that the observed spectra can either correspond to a fully concentrated chromophore or a diluted chromophore, or to a mixture of several chromophores. The original LMM with the constraints of Equation (3) does not consider dilution. It can, however, be extended to this scenario by relaxing the sum-to-one constraint, which can be performed by considering a zero endmember with the LMM, called a shadow endmember.

2.2. Unmixing Strategy

Data unmixing requires inversion of a mixing model that can be either linear or nonlinear, depending on the hypotheses. The study focus on unsupervised methods to be robust to the various photoacoustic systems that exist. Unsupervised algorithms have been devised in remote sensing to extract endmembers and to estimate the abundance matrix. Unsupervised methods, such as group lasso with unit sum and positivity constraints (GLUP) [17] and vertex component analysis (VCA) [18], can do both of these simultaneously, subject to constraints (3). Other unsupervised methods, such as N-FINDR [19], only extract the endmembers. A supervised unmixing algorithm is then required to calculate the abundances. The fully constrained least-square algorithm (FCLS) [20] usually provides satisfactory performance as long as the endmembers are accurately extracted. FCLS can be used with GLUP, VCA, and N-FINDR [17], or with any other strategy for extraction of the endmembers. In this study, we use FCLS to estimate the abundance matrices.

In the following, we use GLUP, VCA (that was already used in photoacoustic imaging in Reference [21]), and N-FINDR to estimate the endmembers in experimental multispectral photoacoustic data to evaluate the interest of remote sensing methods to process these data. We compare their performances with the spatio-spectral mean-shift (SSM-S) algorithm [22] devised by some of the present authors.

3. Method

3.1. Pre-Processing

The pre-processing introduced in this section aims to discriminate pixels of interest from the background. Among the possible strategies, background pixels can be determined by thresholding filtered data with a Sobel filter. The Sobel threshold ${\lambda }_S$ can be set as follows:

$${\lambda }_S=2\ast \sqrt{{\overline{\nabla }}_S},$$

(4)

where ${\overline{\nabla }}_S$ is the mean of the Sobel gradient magnitude [23]. For our application, only the edges detected by the Sobel filter were used to calculate ${\overline{\nabla }}_S$. Threshold ${\lambda }_S$ was applied to the sum over all of the wavelengths of the multispectral photoacoustic image, to design a mask where the background pixels were set to 0. The estimation methods introduced hereafter were only applied to pixels out of the background.

3.2. Endmember Extraction

3.2.1. GLUP Algorithm

Group lasso with unit sum and positivity constraints assumes that the endmembers are not known but are present in the form of some pure (unmixed) pixels in the image [17]. Based on this assumption, the LMM (2) can be reformulated as follows:

$$\begin{array}{cc}\hfill {\mathit{A}}_i& ={\displaystyle \sum _{j=1}^N}{u_G}_{ji}{\mathit{A}}_j+{\mathit{g}}_i,\forall i\in \{1,\dots ,N\},\hfill \end{array}$$

(5)

where ${u_G}_{ji}$ is the abundance of ${\mathit{A}}_j$ in ${\mathit{A}}_i$. On the one hand, if ${\mathit{A}}_j$ is an endmember, the j-th row ${\mathit{U}}_{G_j}$ of matrix ${\mathit{U}}_G$ calculated with GLUP should have nonzero entries that represent the abundance map of ${\mathit{A}}_i$. On the other hand, if ${\mathit{A}}_j$ is a mixed pixel, ${\mathit{U}}_{G_j}$ should have all of its entries equal to zero. This means that ${\mathit{U}}_G$ should have $N-K$ rows of zero, with the other K rows equal to the rows of $\mathit{U}$.

The premise in GLUP is that ${\mathit{U}}_G$ allows the identification of the endmembers in observation $\mathit{A}$ through its nonzero rows. The resulting unmixing problem requires that ${\mathit{U}}_G$ only has a few rows different from zero, in addition to the non-negativity and sum-to-one constraints. This leads to the following convex optimization problem:

$$\underset{{\mathit{U}}_G}{min}\left(\frac{1}{2}\sum _{j=1}^N\parallel \mathit{A}-\mathit{A}{\mathit{U}}_G{\parallel }_{\mathrm{F}}^2+\mu \sum _{j=1}^N\parallel {\mathit{U}}_{G_j}{\parallel }_2\right)$$

(6)

$$\begin{array}{cc}\mathrm{subject}\mathrm{to}:& \left\{\begin{array}{c}{u}_{G_{ji}}\ge 0,\forall i,j\in \{1,\dots ,N\}\hfill \\ {\displaystyle \sum _{j=1}^N}{u_G}_{ji}=1,\forall i\in \{1,\dots ,N\},\hfill \end{array}\right.\end{array}$$

where $\mu >0$ is a small regularization parameter that is set by the user, and $\mathit{A}=[{\mathit{A}}_1,\dots ,{\mathit{A}}_N]$ is the matrix of all pixels to unmix. The first term in Eqtuation (6) ensures that the observations match the model in Equation (5), and the second term is the group lasso regularizer that induces sparsity by possibly driving several rows of ${\mathit{U}}_G$ to zero [24]. The minimization problem has constraints to ensure that the abundances obey the positivity and the sum-to-one constraints. This optimization problem can be solved with a primal-dual method; see Reference [17] for details.

To conclude, GLUP allows the identification of the endmembers in $\mathit{A}$ by identification of the nonzero rows in ${\mathit{U}}_G$. GLUP also provides the estimated abundances that correspond to the nonzero rows in the estimated matrix ${\mathit{U}}_G$. In this paper, similarly to Reference [17], we did not use the abundances estimated by GLUP. For a fair comparison with other algorithms, we re-estimated the abundances with the FCLS algorithm from the endmember spectra estimated using GLUP.

3.2.2. VCA Algorithm

Vertex component analysis also assumes the presence of pure pixels in the data [18]. The driving principle of the VCA algorithm is to project the data onto a direction orthogonal to the subspace spanned by the endmembers already extracted. The new extracted endmember is the farthest signal in this projection domain. Considering this endmember, a new subspace is calculated, and the same procedure is iteratively performed until the preset number K of endmembers to extract is reached.

The first step of the VCA algorithm consists of determination of the initial subspace. Two methods are recommended, which depend on the signal-to-noise ratio. If the signal-to-noise ratio is greater than ${\lambda }_V$ in Equation (7), this first subspace is calculated using the singular value decomposition algorithm [25]. Otherwise, the subspace considered is constructed using principle component analysis. The threshold ${\lambda }_V$ is defined in Reference [18] as follows:

$${\lambda }_V=15+10{log}_{10}\left(K\right).$$

(7)

Let us use ${\mathit{S}}_k$ to denote the subspace available at iteration k. A vector ${\mathit{v}}_k$ orthonormal to ${\mathit{S}}_k$ is calculated as follows:

$${\mathit{v}}_k=\frac{{\mathit{r}}_k-{\mathit{S}}_k{\mathit{S}}_k^{+}{\mathit{r}}_k}{\parallel {\mathit{r}}_k-{\mathit{S}}_k{\mathit{S}}_k^{+}{\mathit{r}}_k\parallel },$$

(8)

where ${\mathit{r}}_z$ is a zero-mean random vector, and ${\mathit{S}}_z^{+}$ is the pseudo inverse of ${\mathit{S}}_k$. The observed spectra $\mathit{A}$ are then projected onto direction ${\mathit{v}}_k$:

$${\mathit{f}}_k={\mathit{v}}_k^{\top }\mathit{A}.$$

(9)

The largest entry of ${\mathit{f}}_k$ in absolute values allows the designation of the spectrum in $\mathit{A}$ to be considered as a new endmember. This endmember is then added to the set of endmembers that have already been extracted, to define the subspace ${\mathit{S}}_{k+1}$ that is considered at the next iteration. The procedure is stopped when the number K of desired endmembers is reached. They are stored in matrix ${\mathit{E}}_V$, and the abundance matrix ${\mathit{U}}_V$ is calculated by projecting $\mathit{A}$ onto ${\mathit{E}}_V$. For a fair comparison with other algorithms, we estimated the abundances with the FCLS algorithm from the endmember spectra ${\mathit{E}}_V$ estimated with VCA.

3.2.3. N-FINDR Algorithm

N-FINDR also assumes that a pure pixel for each chromophore to unmix is present in the dataset [19,26]. The first step is the random generation of a set of K endmembers to produce the matrix ${\mathit{E}}_0$ of endmembers that are available at iteration 0. At each iteration k, the following volume is calculated:

$${\mathit{V}}_k=\frac{|det({\mathit{E}}_k\left)\right|}{(K-1)!},$$

(10)

$$\begin{array}{cc}\mathrm{with}& {\mathit{E}}_k=\left[\begin{array}{ccc}1& \dots & 1\\ {\mathit{E}}_k^{\left(1\right)}& \dots & {\mathit{E}}_k^{\left(K\right)}\end{array}\right],\end{array}$$

where ${\mathit{E}}_k^{\left(p\right)}$ is the p-th endmember in the matrix of endmembers ${\mathit{E}}_k$ available at the k-th iteration. All of the pixels ${\mathit{A}}_i$ in $\mathit{A}$ take place successively in ${\mathit{E}}_k$, where they are successively substituted to all possible columns p. Volume ${\mathit{V}}_{k+1}$ is updated as follows, for instance:

$${\mathit{V}}_{k+1}=\frac{\left|det\left[\begin{array}{ccccc}1& \dots & 1& \dots & 1\\ {\mathit{E}}_k^{\left(1\right)}& \dots & {\mathit{A}}_i& \dots & {\mathit{E}}_k^{\left(K\right)}\end{array}\right]\right|}{(K-1)!}.$$

(11)

If the new volume is greater than the previous one, ${\mathit{V}}_k$ and ${\mathit{E}}_k$ are updated with the new values. The procedure is stopped when the N pixels of the dataset have been tested. The resulting endmembers are stored in matrix ${\mathit{E}}_N$. No abundance coefficients are calculated with this method, and the FCLS can be used to this end.

3.2.4. SSM-S Algorithm

Spatio-spectral mean-shift is a clustering method that was introduced by Reference [22] for chromophore discrimination in multispectral photoacoustic images. It is based on a spatio-spectral regularization, to cluster the pixels that are spatially and spectrally close.

Consider a pixel ${\mathbf{x}}_i$. Its neighbors in spatial dimensions at a radial distance less than $R_S$ are first considered. These pixels ${\mathbf{x}}_j=\left[{\mathrm{s}}_j,{\mathit{A}}_j\right]$ have to satisfy the following:

$$\frac{1}{R_S}{\parallel {\mathrm{s}}_i-{\mathrm{s}}_j\parallel }_2\le 1.$$

(12)

Assuming that spectra of the same chromophore are close, we also consider that two spectra can be from the same chromophore if their distance is less than $R_{\lambda }$; namely,

$$\frac{1}{R_{\lambda }}{\parallel {\mathit{A}}_i-{\mathit{A}}_j\parallel }_{\infty }\le 1.$$

(13)

Note that we found it interesting in Equation (13) to consider the infinite norm, to discriminate more effectively spectra that might differ only in narrow frequency bands.

The SSM-S only applies to pixels that satisfy both constraints of Equations (12) and (13). It consists of updating their locations and spectra as follows:

$${\mathbf{x}}_i^{[t+1]}=\frac{{\displaystyle \sum _{\begin{array}{c}j=1\end{array}}^N}{g}_S({\mathrm{s}}_i^{\left[t\right]};{\mathrm{s}}_j^{\left[t\right]})·g_{\lambda }({\mathit{A}}_i^{\left[t\right]};{\mathit{A}}_j^{\left[t\right]})·{\mathbf{x}}_j^{\left[t\right]}}{{\displaystyle \sum _{\begin{array}{c}j=1\end{array}}^N}{g}_S({\mathrm{s}}_i^{\left[t\right]};{\mathrm{s}}_j^{\left[t\right]})·g_{\lambda }({\mathit{A}}_i^{\left[t\right]};{\mathit{A}}_j^{\left[t\right]})},$$

(14)

where $g_S({\mathrm{s}}_i^{\left[t\right]};{\mathrm{s}}_j^{\left[t\right]})=1$ and $g_{\lambda }({\mathit{A}}_i^{\left[t\right]};{\mathit{A}}_j^{\left[t\right]})=1$, if Equations (12) and (13) hold, respectively, and 0 otherwise. This iterative procedure is applied to all of the pixels ${\mathbf{x}}_i$ until convergence; see Reference [22] for details. Unlike the GLUP, VCA, and N-FINDR algorithms, note that SSM-S uses spatial information in addition to spectral information. SSM-S is studied here because of this particular capacity which is not often taken into account to process multispectral photoacoustic data. Updating of Equation (14) leads to clusters of pixels, with each defined by a centroid. Each centroid can be considered as an endmember and stored in a matrix, which is denoted by ${\mathit{E}}_S$. As no abundance coefficients are calculated with this method, the FCLS can be used to this end.

Before concluding this section, we want to draw attention to the following: that, unlike the N-FINDR and SSM-S algorithms, GLUP and VCA converge to stationary points that can vary significantly depending on their initialization.

3.3. Abundances Estimation

The FCLS algorithm is widely used in remote sensing to estimate abundances, based on endmember knowledge. It considers the constraints of Equation (3) when solving the following optimization problem:

$$\underset{\mathit{U}}{min}\left(\frac{1}{2}\sum _{j=1}^N{\parallel \mathit{A}-\mathit{U}\mathit{E}\parallel }_{\mathrm{F}}^2\right)\prime$$

(15)

$$\begin{array}{cc}\mathrm{with}& \left\{\begin{array}{c}{u}_{ji}\ge 0,\forall i,j\hfill \\ {\displaystyle \sum _{j=1}^N}{u}_{ji}=1,\forall i\hfill \end{array}\right.\end{array}$$

4. Materials

4.1. Acquisition System

A commercial multispectral photoacoustic system (Vevo LAZR; Visualsonics, Fujifilm) was used for this study; see Figure 2. The optical source of this system is a Nd:YAG pulsed laser with a pulse duration of 5 ns and a 20 Hz repetition rate, which is coupled to an optical parametric oscillator to select the transmission wavelengths [27]. The wavelength can be set from 680 nm to 970 nm. The photoacoustic probe that was used (LZ400) is composed of 256 elements, which can acquire ultrasound signals in the frequency range from 18 MHz to 38 MHz. The imaging depth is approximately $1.5$ cm and the imaged region of interest is around 1 cm. The images are reconstructed with a delay-and-sum algorithm. For every acquisition, the corresponding fluence is saved. For the acquisitions presented in this study, the fluence is varying of only 2% between the minimal and maximal fluences. As the variation is low and the proposed methods are based on unsupervised endmembers extraction with relative comparison between spectra, the images were not normalized by the fluence.

The data were acquired using the full range of available wavelengths $(680$ nm$,970$ nm) with 1 nm steps, which required less than 1 min of acquisition time. The time to switch from one wavelength to the other and to acquire the corresponding image is around 0.2 s. Usually, in multispectral photoacoustic imaging, from 5 to 10 wavelengths are used because the information collected is then sufficient to discriminate or quantify biological chromophores with an acceptable processing time. To get closer to this framework, only eight of the acquired wavelengths were selected to construct our dataset, from 680 nm to 820 nm, with 20 nm steps. These wavelengths were chosen because the used optical absorbers (blue and green inks) can be discriminated in this range (see Figure 3).

The algorithms were computed with MATLAB (The MathWorks Inc., Natick, MA, USA) on a computer with a 16-GB memory and an Intel Core i7-5500U CPU at 2.40 GHz.

4.2. Imaged Phantoms

The homemade colored phantoms used to construct the dataset were composed of three different 4% colored agar parts. The ones on the left side in Figure 4a were prepared using only blue ink (400 g water, 16 g agar, 380 $\mathsf{\mu }$L ink), and the ones on the right with only green ink (400 g water, 16 g agar, 950 $\mathsf{\mu }$L ink). Both of these are considered as relative ink concentrations equal to 1 (i.e., pure chromophores). The blue and green ink quantities were not equal because both of these inks do not show the same maximum photoacoustic signal amplitude. Based on a previous study [28], the ink quantities were chosen to obtain the same maximum photoacoustic signal amplitude, on the range from 400 nm to 1200 nm, for both of these pure chromophores. The used blue and green inks spectra, measured with a Perkin Elmer Lambda900 spectrometer for the study [28], are shown in Figure 3, on the wavelength range of interest for the present study.

In Figure 4b, the sample of shown spectra seems to highlight that the fully concentrated blue and green parts do not exhibit the same, even approximately, spectral intensities for both phantoms. However, this is only due to the spectra shown here. Nevertheless, the Figure shows that the shape of blue or green spectra are quite similar which is the property used by the proposed methods in this study.

4.2.1. Chromophore Dilution

The central part of the left phantom in Figure 4a is composed of a $0.53$ dilution of the blue ink relative to the pure concentration indicated above. Therefore, here 100 g water, 4 g agar, and 50 $\mathsf{\mu }$L blue ink were used. This phantom is referred to as B-Bdil-G(Vevo) in the following.

4.2.2. Mixing of the Chromophores

The central part of the second phantom, shown in Figure 4a-right, is a mix of the blue and green inks (50 g water, 2 g agar, 20 $\mathsf{\mu }$L blue ink, 80 $\mathsf{\mu }$L green ink). This means that the mix corresponds to a blue concentration of 0.42 and a green concentration of 0.67, relative to the pure concentrations indicated above. This phantom will be referred to as B-mix-G(Vevo) in the following.

4.3. Performance Evaluation on Phantoms

The performances of the GLUP/FCLS, VCA/FCLS, N-FINDR/FCLS, and SSM-S/FCLS strategies were evaluated on the phantoms described in Section 4.2. This means that FCLS was used to estimate the concentrations of the endmembers extracted beforehand with the other methods presented in Section 3.2. The resulting abundance maps were compared to the ground truth.

For each abundance map of each endmember, a mean concentration was calculated for each part of the imaged region, as illustrated in Figure 5. The illustrative phantom considered in Figure 5 is composed of a pure blue chromophore, a dilution of this chromophore, and a pure green chromophore (from left to right). Considering the abundances related to the blue endmember, three average values were calculated: one for each part, as circled by the yellow dotted line in Figure 5. As we assume that the pure blue chromophore, which corresponds to the endmember of interest, is present in the image, the mean values were normalized by their maximum, as circled by the orange dashed line in Figure 5. Where the data would be perfectly unmixed, the maximum value should correspond to the pure blue part of the phantom (i.e., the left part), which is then set to 1 after normalization. The other normalized values should correspond to the relative concentrations of the pure chromophore, i.e., the relative concentration of blue in the illustrative example in Figure 5.

This calculation was performed to estimate the abundance maps provided by the different unmixing strategies. Indeed, by considering the light absorption and the diffusion in tissue, each pixel of similar region exhibit similar spectra shape tendency (see Figure 4b). For each pixel, a particular concentration is then evaluated. The normalized mean values were compared with the ground truth to evaluate the performances of the algorithms on each whole part.

4.4. In Vivo Data Acquisitions

Images were acquired with the Vevo LAZR system. For imaging the mice tumors, the BALBc MMTV-Neu transgenic mouse model was used, which expresses the Her-2/neu V664E oncogene [29]. This model develops spontaneous mammary tumors that are visible with the photoacoustic system after 17 weeks (between 5 mm${}^3$ and 10 mm${}^3$ in size). Palpable tumors are detected around 20 weeks old; see Figure 6, where the white arrows indicate some of the tumors, and the green square denote the imaged one. A 20-week-old mouse was used for the imaging, where the region of interest was shaved the day before using a commercial hair remover cream (Veet, Cream hair remover; Reckitt Benckiser, UK) to avoid any imaging noise coming from the mouse hair. The mouse was anesthetized with a mix of 3.5% isoflurane and 96.5% oxygen before the photoacoustic acquisitions, using an anesthesia mask. During the imaging procedure, the mouse was maintained asleep through inhalation of a mix that contained 1.5% isoflurane (Figure 6a), to minimize movements during the acquisitions at the different wavelengths. Heartbeat and breathing were monitored by several devices connected to the mouse paws (Figure 6b, red arrows) and the anesthesia mask, respectively. This allowed monitoring of the health of the mouse and the acquisition of images at the same position. In this way, heartbeat and breathing motions were avoided as much as possible. The mice were maintained in the specific pathogen free animal facility at the Cancer Research Center of Lyon, at Center Léon Bérard. All of the mice experiments were performed in accordance with the animal care guidelines of the European Union and were validated by the local Animal Ethics Evaluation Committee and the French Ministry of Higher Education and Research (C2EA-15 and 02296.02).

The dataset used to test this strategy was composed of 15 images that were acquired from 680 nm to 960 nm, in 20 nm steps. The results were compared to the Vevo LAZR sO${}_2$ map. To this end, acquisitions on the same tumor were also performed using the Vevo LAZR Oxy-Hemo mode.

4.5. sO${}_2$ Calculation with Vevo LAZR Oxy-Hemo Mode

The Oxy-Hemo mode acquires two photoacoustic images, at 750 nm and 850 nm. With both of these images, the concentrations of deoxygenated hemoglobin (Hb) and oxygenated hemoglobin (HbO${}_2$) were measured, which are expressed as [Hb] and [HbO${}_2$], respectively [30]. The oxygen saturation was then calculated as follows:

$${\mathrm{sO}}_2=\frac{{[\mathrm{HbO}}_2]}{\left[\mathrm{Hb}\right]+{[\mathrm{HbO}}_2]}.$$

(16)

This ratio is classically used in biomedical applications to provide oxygenation information [31,32]. Both of the concentration maps were saved, as well as the Vevo LAZR sO${}_2$ displayed map, for further comparisons. The sO${}_2$ acquisition is conducted just after the acquisitions on the whole range of wavelength, allowing then a fair comparison of the estimated concentrations.

5. Results

As indicated in Section 3.2, the N-FINDR and SSM-S algorithms converge to a single set of estimated endmembers, regardless of their initialization. On the contrary, the GLUP and VCA methods can converge to different sets of endmembers, which depend on the initialization. For the acquisitions with the phantoms, the results presented for GLUP and VCA are the best ones that were achieved. For the preliminary in vivo results, the sets of endmembers obtained with the different techniques are presented and discussed.

5.1. Chromophore Dilution

The endmembers estimated with GLUP, VCA, N-FINDR, and SSM-S are presented in Figure 7a–d, respectively. They allow us to discriminate the pure blue and green chromophores even if they do not exactly fit the reference spectra measured with a spectrometer (Figure 3). Indeed, as already mentioned, light absorption and diffusion in tissue impact the pure chromophore spectrum. The shape tendency is, however, kept. The shadow endmember is also plotted in each panel.

The abundance maps and the estimated mean concentrations are presented in Figure 8, for both the endmembers. All of these methods localize both of the pure chromophores well. The green abundance maps clearly highlight the green part, i.e., that on the right. The SSM-S achieves the best performance, as it leads to small green abundances for both parts of the blue phantom. Indeed, the mean concentrations are 0.05 and 0.07 for the blue and the diluted blue parts, respectively. The blue abundance maps are, nevertheless, also characterized by relatively large values in the green part. The SSM-S blue abundance map shows the smallest value of 0.26, and estimates a dilution value of 0.51, as compared to 0.53 for the ground truth. In conclusion, for this dataset, SSM-S/FCLS achieves the best performance.

5.2. Mixing of the Chromophores

The unmixing methods were also tested with mixed chromophores using the B-mix-G(Vevo) dataset. The endmembers extracted by the different methods are presented in Figure 9. The extracted spectra differ from one strategy to another, but all blue spectra have similar evolutions, which matches the ground truth multispectral characteristics, in term of curve evolution, of this blue ink (Figure 3). The shadow endmembers are also plotted.

The abundance maps and the estimated mean concentrations are presented in Figure 10, for both of the endmembers. The blue abundance maps calculated with the endmembers extracted by GLUP and VCA localize the pure blue chromophore part (left) well, while the pure green part is more difficult to delimit. Indeed, the abundance values are similar in the central and right parts. For the N-FINDR and SSM-S blue abundance maps, the three parts are more visible. On the SSM-S abundance maps, the green part is characterized by a low concentration of 0.10, which is close to the ground truth value. For the green abundance maps, all of the methods allow discrimination of the blue part (left) and the mixed part (center). The VCA abundance map does not allow for clear discrimination of the mixed (central) and green (right) parts, while this is possible with the other methods.

For the estimated mean concentrations in Figure 10, the first result to be noted is the blue concentration of the mixed chromophore estimated using VCA, which is equal to the ground truth of 0.42, as circled by the red dashed line. Nevertheless, the other concentrations estimated with VCA are larger than the ground truth concentrations: 0.89 for the green in the mixed chromophore instead of 0.67, 0.27 for the blue in the pure green chromophore instead of 0, and 0.20 for the green in the pure blue chromophore instead of 0. Consider now the performance of the SSM-S/FCLS, as circled with the orange dashed line in Figure 10. The estimated concentrations in the mix of chromophores are 0.35 for the blue and 0.73 for the green, instead of 0.42 and 0.67 for the ground truth, respectively. This corresponds to an error of less than 7%, which is already an interesting result in terms of the targeting of biological applications. Note that the estimated abundances where the ground truth concentration equals 0 were estimated at values less than 0.1 by SSM-S, which is a satisfying performance compared to the other methods. As a conclusion, for this dataset, SSM-S/FCLS achieves the best performance.

5.3. Preliminary In Vivo Results

The analyzed strategies were used on biological tissues to calculate the oxygen saturation rate (sO${}_2$). This measurement is of great interest for various medical applications, such as the follow-up of tumors and the evaluation of tissue aging [33].

The pre-processing steps described in Section 3.1 were followed. Then, the endmember estimation methods were applied to extract two endmembers, which corresponded to the Hb and HbO${}_2$ chromophores (the theoretical Hb and HbO${}_2$ spectra are shown Figure 11). Finally, the abundance maps of both of these endogeneous chromophores were estimated using the FCLS algorithm. The sO${}_2$ at each pixel was determined using the resulting abundance maps and Equation (16). These results were compared to the one provided by the Vevo LAZR. Indeed, the Vevo LAZR can provide abundance maps of the Hb and HbO${}_2$ chromophores. We used these abundance maps with Equation (16) to calculate the oxygen saturation rates; see Figure 12. The authors are aware that Vevo LAZR Oxy-Hemo mode is not the real ground truth. However, this scanner represents a significant reference since it is typically used for sO${}_2$ measurements in photoacoustic imaging.

Figure 13, left column, shows the endmember spectra estimated with GLUP, VCA, N-FINDR, and SSM-S. Recall that the GLUP and VCA algorithms can converge to different sets of endmembers, depending on their initialization. The GLUP and VCA spectra represented in Figure 13 are those that then lead to the best sO${}_2$ maps compared to the Vevo LAZR map, which are considered as the relative ground truth. FCLS was used with all of these endmember spectra, to estimate the abundance maps of Hb and HbO${}_2$. The sO${}_2$ maps were then calculated using Equation (16); see Figure 13, center column. With this strategy, it appears that many pixels have strong sO${}_2$ values, when VCA and N-FINDR algorithms are used. In addition, the endmembers extracted with both these algorithms are really similar, but not exactly the same, even if it is quite hard to see their differences on Figure 13. However, the sO${}_2$ maps calculated with this strategy are different, which highlights the endmembers differences.

To characterize the performances of these methods, the relative deviations between each estimated sO${}_2$ map (see Figure 13, center column) and the Vevo LAZR map calculated with Equation (16) (see Figure 12, right column) were computed. This relative deviation defines the absolute difference between the Vevo LAZR sO${}_2$ values and each of the estimated sO${}_2$ values, for each pixel. This is provided in %, as the sO${}_2$ values are also in %. The relative deviation maps are presented in Figure 13, right column. Note that skin, tumor, and other biological tissues are present in the imaged area. The mean relative deviation was first calculated for all of these tissues. However, as the study was focused on the tumor only, the mean relative deviation was also calculated for the restricted regions delimited by the green rectangles in Figure 13. These relative deviations for all of these tissues and just for the tumor are in the black and green characters on the right in Figure 13, respectively. For both regions, the best estimations were reached with GLUP and SSM-S, which give similar relative deviations: 24.04% and 15.86% for GLUP, and 24.07% and 15.29% for SSM-S.

6. Discussion

This study experimentally evaluates the interest of using multispectral photoacoustic imaging to estimate chromophore concentration. To tackle such objectives, several methods coming from the remote sensing field have been tested to evaluate their potential performances in the photoacoustic domain. Principally, the constraints of Equation (3) can benefit also for chromophore concentration evaluation purpose. The interests and disadvantages of each method, highlighted by the in vivo proof of concept to evaluate sO${}_2$, are here discussed.

First, the variabilities of GLUP and SSM-S have to be analyzed. As already mentioned, GLUP and VCA can converge to different sets of endmembers, which depend on their initialization. Here, we focus only on GLUP as VCA, regardless, gives inaccurate sO${}_2$ values. Figure 14 and Figure 15 illustrate this phenomenon with the results from 200 runs for the GLUP method. If GLUP can lead to small relative deviations, it should be noted that only 17 runs provided mean relative deviations less than $25\%$. In contrast, 59 runs provided mean errors greater than $50\%$. It is, however, interesting to note that the endmembers extracted with VCA and N-FINDR have the closest spectra to the theoretical ones [1]; see Figure 16a. This observation should be considered with caution though, as the theoretical spectra shown do not take into account the characteristics of the acquisition system, which might modify the acquired spectra. As well as for the inks study, light absorption and diffusion in tissue also impact the pure chromophore spectrum. This could explain the bad results given by these two methods even if the extracted spectra seem the similar ones.

Finally, considering its robustness, SSM-S is an interesting solution, as it converges to a unique solution while achieving a small relative deviation. However, it is important to note that the endmembers extracted with SSM-S do not look like the theoretical ones; see Figure 16a. As mentioned earlier, this could be due to the characteristics of the acquisition system and the light absorption and diffusion in tissue. In addition, in tumor tissues, pure Hb and HbO${}_2$ are not the most frequently present chromophores. The Vevo LAZR sO${}_2$ values (Figure 12) confirm this observation as they are mainly comprised between 25% and 75%, and then correspond to a mix of pure chromophores. This means that the endmembers that were actually extracted probably represent mixed chromophores instead of the pure Hb and HbO${}_2$. Regardless, the results obtained with the SSM-S/FCLS strategy are very encouraging in allowing accurate measurements of sO${}_2$ with photoacoustic imaging.

Table 1. Summary of the average sO${}_2$ relative deviation values. The abundances are estimated using the FCLS algorithm.

Endmember Extraction	Whole Image	Limited to Tumor
GLUP	$24.04\%$	$15.86\%$
VCA	$58.79\%$	$38.96\%$
N-FINDR	$76.81\%$	$51.63\%$
SSM-S	$24.07\%$	$15.29\%$
Theoretical spectra	21.52%	13.33%

summarizes all of the results, the abundances are calculated with the FCLS algorithm, and the mentioned methods are the ones used to extract the endmembers. This table shows that the minimum mean relative deviation was obtained with the theoretical spectra, but the mean relative deviations remain large, as these spectra do not take the characteristics of the acquisition system into account as well as the light absorption and diffusion in tissue. Nevertheless, this result is interesting, as it is independent of the Vevo LAZR Oxy-Hemo mode, and it provides a convenient basis for comparing results. In this sense, given the proximity to the performance obtained with SSM-S/FCLS, these complementary experiments confirm the potential of this strategy to quantify photoacoustic in vivo data. GLUP/FCLS also have similar performance but only for some runs, as mentioned before. This strategy could then be more studied and compared to SSM-S/FCLS on various cases.

However, the proposed strategies of hyperspectral imaging suffer from spectral coulouring effect [34,35]. Indeed, in the proposed strategies, the endmembers are evaluated without taking into account any fluence evolution through the sample. Such effect has demonstrated to bias the sO${}_2$ estimation and would required deeper study to be conducted. However, in the context of this work and its comparison to the sO${}_2$ evaluation of the Vevo LAZR, we did not evaluate its impact.

These results are encouraging to be used for various clinical applications requiring chromophore quantification. Indeed, the evaluation of chromophore concentrations is of major interest for various applications, like concentration of contrast agent in the body [5], calculation of the concentration of oxygenated and deoxygenated blood for oxygenation rate evaluation [6,7], or cancer tissue evaluation [8]. The endmembers extraction could be either automatic or manually selected by the user coupled with a fast inverse problem (FCLS).

7. Conclusions and Perspectives

Several methods for characterization of multispectral photoacoustic data have been considered, with the aim to estimate chromophore concentrations. These strategies mainly consist of two steps: (1) an endmember extraction step, to estimate the spectrum of each pure chromophore to be quantified; and (2) an abundance calculation step, to evaluate their concentrations. For the first step, different algorithms were compared (GLUP, VCA, N-FINDR, and SSM-S). The second step was performed using the FCLS algorithm.

For nonbiological tissues and using the Vevo LAZR acquisition system, we demonstrated that the SSM-S/FCLS procedure provides accurate estimates of chromophore concentrations for dilutions and mixtures. The positivity and sum-to-one constraints imposed by the FCLS, and coming from the remote sensing field, are interesting to take into account for chromophore concentration evalution. In some cases, the sum-to-one constraint is, however, questionable. It is a strong assumption, particularly when dilutions are imaged. To relax this constraint, a shadow endmember was introduced in this study, among the other possible strategies.

The methods were tested on a mouse tumor, where the sO${}_2$ calculation was performed and compared to the sO${}_2$ values calculated by the Vevo LAZR Oxy-Hemo mode. GLUP and SSM-S provided equivalent optimal performance to extract endmembers. GLUP algorithm should, however, be used with caution because of its variability. Nonetheless, it again highlights the interest of the constraints of Equation (3) in the photoacoustic field for chromophore concentration evaluation purpose. The SSM-S algorithm, which converges to a unique solution, appears to be more appropriate for endmember extraction, but further study on various cases should be done to quantify the performances of SSM-S/FCLS and GLUP/FCLS on various cases.

The endmember extraction and sum-to-one constraint relaxation still needs to be studied in more detail, but these primary results with biological tissues are encouraging. Moreover, the ground truth of the oxygenation rate in the imaged tumor was not known. It is, thus, difficult to know exactly which method is the most accurate to estimate the sO${}_2$ values. In addition, considering the presence of only Hb and HbO${}_2$ in the imaged area is a strong hypothesis. The water and other tissue absorbances should probably be taken into account. Further, other wavelengths can be tested to characterize their influence on the results. In the future, these experiments should be correlated with blood extraction and sO${}_2$ evaluation using blood measurements to have the real ground truth.