Fast unmixing of multispectral optoacoustic data with vertex component analysis

Optics and Lasers in Engineering 58 (2014) 119–125

Contents lists available at ScienceDirect

Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng

Fast unmixing of multispectral optoacoustic data with vertex component analysis X. Luís Deán-Ben a, Nikolaos C. Deliolanis b, Vasilis Ntziachristos a, Daniel Razansky a,n a Institute for Biological and Medical Imaging, Technische Universität München and Helmholtz Zentrum München, Ingolstädter Landstraße 1, 85764 Neuherberg, Germany b Fraunhofer Institut für Produktionstechnik und Automatisierung IPA, Theodor-Kutzer-Ufer 1-3, 68167 Mannheim, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 1 September 2013 Received in revised form 18 January 2014 Accepted 30 January 2014

Multispectral optoacoustic tomography enhances the performance of single-wavelength imaging in terms of sensitivity and selectivity in the measurement of the biodistribution of specific chromophores, thus enabling functional and molecular imaging applications. Spectral unmixing algorithms are used to decompose multi-spectral optoacoustic data into a set of images representing distribution of each individual chromophoric component while the particular algorithm employed determines the sensitivity and speed of data visualization. Here we suggest using vertex component analysis (VCA), a method with demonstrated good performance in hyperspectral imaging, as a fast blind unmixing algorithm for multispectral optoacoustic tomography. The performance of the method is subsequently compared with a previously reported blind unmixing procedure in optoacoustic tomography based on a combination of principal component analysis (PCA) and independent component analysis (ICA). As in most practical cases the absorption spectrum of the imaged chromophores and contrast agents are known or can be determined using e.g. a spectrophotometer, we further investigate the so-called semi-blind approach, in which the a priori known spectral profiles are included in a modified version of the algorithm termed constrained VCA. The performance of this approach is also analysed in numerical simulations and experimental measurements. It has been determined that, while the standard version of the VCA algorithm can attain similar sensitivity to the PCA–ICA approach and have a robust and faster performance, using the a priori measured spectral information within the constrained VCA does not generally render improvements in detection sensitivity in experimental optoacoustic measurements. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Optoacoustic imaging Photoacoustic imaging Hyperspectral unmixing Molecular imaging

1. Introduction Multispectral optoacoustic tomography (MSOT) extends the capacities of monochromatic (single-wavelength) imaging by making use of the spectral signatures of the various intrinsic tissue chromophores and extrinsically administered contrast agents present in the imaged tissue [1,2]. Thereby, whereas single-wavelength excitation provides high-resolution mapping of the background absorption in biological tissues for depths up to a few centimeters [3,4], a better mapping of the specific absorbers present is attained by recording images at multiple wavelengths. In this way, it is possible to image e. g. the oxygenation level of blood or the distribution of externally administered photo-absorbing agents, which facilitates emergence of a powerful set of functional and molecular optoacoustic imaging applications [5–8].

n

Corresponding author. E-mail address: [email protected] (D. Razansky).

http://dx.doi.org/10.1016/j.optlaseng.2014.01.027 0143-8166 & 2014 Elsevier Ltd. All rights reserved.

Improvements in the quantification capacity and sensitivity of MSOT are the two main challenges presented by this technique for its applicability in molecular imaging. Quantitative MSOT of photoabsorbing probes requires the accurate determination of their spatial distribution within the imaged tissue in relative amounts or actual units of concentration, where imaging in arbitrary units can be sufficient if the system is properly calibrated. Quantification errors are originated in the reconstructed singlewavelength images mainly due to effects of the ultrasonic transducer [9,10], acoustic attenuation [11,12] and optical attenuation [13,14]. If the single-wavelength images are properly corrected to avoid quantification errors, the multispectral processing procedure determines the sensitivity of MSOT, i.e., its capacity to image low concentrations of the probe. The optical absorption within a biological tissue is a linear combination of the absorption contributed by the different tissue chromophores and external contrast agents, whereas each absorbing substance has a characteristic absorption spectrum. In the simplest form, the distribution of a given substance within a biological tissue can in principle

120

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

be imaged by exciting with an optical wavelength close to its peak absorption [15]. However, this procedure is limited to those cases in which the absorption of the substance of interest is considerably higher than the background absorption. If the optical absorption of the different chromophores is comparable, the spectral profile of each must be taken into account in order to efficiently distinguish the spatial distribution of the substance of interest using images taken at different optical wavelengths. This procedure is further referred to as spectral unmixing [16]. The simplest spectral unmixing method consists in calculating the difference of the images taken at two wavelengths corresponding to high and low absorption of the substance of interest [17]. This method is especially convenient for unmixing fluorescent dyes having a sharp drop in their optical absorption in a relatively narrow spectral window. However, the detection sensitivity of the various substances can be substantially increased by considering images at more wavelengths. If the spectral signature of the different components is known, which is usually the case as the background absorption is mainly due to blood (oxygenated and deoxygenated hemoglobin), this information can subsequently be used in direct spectral inversion procedures such as least-square fitting unmixing [6,18]. A different approach, the so-called blind unmixing, does not include any a priori spectral information but the spectra are instead recovered as a part of the inversion procedure. Blind unmixing in MSOT can be achieved by e.g. combining principal component and independent component analysis methods (PCA–ICA) [19], which was shown to perform better in terms of sensitivity compared to the spectral fitting approach. The reasons for this enhanced performance are mainly related to errors in the reconstructed spectral profiles of the images and can be multifarious. For instance, errors in the reconstruction due to wrong assumptions, wrong normalization of the laser energy incident upon the tissue or the spectral dependence of the light attenuation might influence the spectral signature of the optoacoustic tomographic reconstructions. The computational time of the unmixing procedure has also become an important factor with the emergence of fast (real-time) optoacoustic imaging systems [15,20–23], so that fast unmixing of multispectral data is essential for the visualization of spectrallydistinctive chromophores. In this paper we suggest a different blind unmixing approach termed vertex component analysis (VCA) [24] for unmixing of multispectral optoacoustic data. Its performance in terms of sensitivity is compared with the previously reported PCA–ICA method. Since in this approach the spectral components are determined individually (separately), we further investigate the effect of including some of the a priori information on the spectral signatures of absorbers in the spectral inversion process, which is herein termed the constrained VCA or semi-blind unmixing approach.

2. Theory 2.1. Spectral fitting For each given laser wavelength λi, the reconstructed optoacoustic image represents the spatial distribution of the optical absorbed energy Hðr; λi Þ. If m different substances contribute to the absorption, Hðr; λi Þ in arbitrary units is given by Hðr; λi Þ ¼ Uðr; λi Þϵ1 ðλi ÞC 1 ðrÞ þ Uðr; λi Þϵ2 ðλi ÞC 2 ðrÞ þ ⋯ þUðr; λi Þϵm ðλi ÞC m ðrÞ;

ð1Þ

with i ¼ 1; 2; …; n and n being the total number of wavelengths employed. Uðr; λi Þ represents the space and wavelength dependent light fluence, ϵj ðλi Þ and C j ðrÞ are the wavelength-dependent molar

extinction coefficient for the j-th absorber and its spatial distribution (concentration) respectively. Here, a homogeneous distribution of the Grueneisen parameter was assumed. For each pixel k of the image, Eq. (1) can be expressed in a matrix form as Hk ¼ ϵk Ck ;

ð2Þ

where H is a column vector with elements H ðiÞ ¼ Hðr k ; λi Þ, C is a column vector with elements Ck ðjÞ ¼ C j ðr k Þ and ϵk is a matrix with elements ϵk ði; jÞ ¼ Uðr k ; λi Þϵj ðλi Þ. Generally, the spectral variation of the light fluence is also spatially dependent. Even if the illumination profile on the surface of the object is kept constant for different wavelengths, the changes in the optical attenuation create a wavelength dependence of the light fluence for deeper locations. Thus, estimation of absorber concentration must be done on a per-pixel basis by means of e.g. spectral fitting, which consists of a least square minimization of the form k

k

Cksol ¼ argmin‖Hkexp  ϵk Ck ‖22 ;

k

ð3Þ

Ck

where Hkexp corresponds to the experimental measurement of the column vector Hk . As a first order approximation, one can assume that the spectral dependence of the light fluence is less significant than the spectral variation of the molar extinction coefficients of the different substances, so that Eq. (2) can be approximated as Hk  U ðr k ÞϵCk ;

ð4Þ

where U ðr k Þ is the mean value of Uðr k ; λi Þ for all the wavelengths and ϵ is a matrix with elements ϵði; jÞ ¼ ϵj ðλi Þ. In this case, Eq. (3) is transformed in arbitrary units to U ðr k ÞCksol ¼ argmin‖Hkexp  ϵCk ‖22 :

ð5Þ

Ck

The minimization in Eq. (5) can be done with the pseudoinverse ϵ þ of matrix ϵ for all the pixels of the image, i.e., CU;sol ¼ ϵ þ Hexp ;

ð6Þ

being CU;sol and Hexp two matrices whose columns are U ðr k ÞCksol and Hkexp respectively. As indicated in Eq. (5), the resulting images obtained by spectral fitting are affected by light fluence variations, so that they must be corrected in order to obtain a quantitative distribution of the absorber concentration. Even though several methods can be used to correct for the light fluence variations [13,25,26], this work is only concerned with the procedure to retrieve CU;sol . 2.2. Vertex component analysis In the previous section, it was shown that the spectral dependence of the light fluence may induce errors in the unmixed images if the spectral variations in the images are directly fitted to the molar extinction coefficients of the absorbing substances, especially for absorbers located deeper in tissue where spectrally-dependent light attenuation may become significant. In practice, the spectral variations of the reconstructed images may also be affected by other factors. For example, even if the signals for each wavelength are normalized to the respective laser pulse energy, the illumination profile may also change with the laser wavelength due to unstable or otherwise low laser beam quality. Furthermore, errors in the tomographic reconstructions may also influence the spectral variations in the images. The blind unmixing methods, including VCA, are expected to show improved sensitivity in those cases. The VCA algorithm consists of a two step procedure [24]. First, the spectral signatures of the different components are calculated whereas distribution of the different components is subsequently obtained by spectral fitting. The spectral signatures correspond to

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

the columns of the matrix ϵVCA given by 0 1 ϵVCA;1 ðλ1 Þ ϵVCA;2 ðλ1 Þ ⋯ ϵVCA;m ðλ1 Þ B C B ϵVCA;1 ðλ2 Þ ϵVCA;2 ðλ2 Þ ⋯ ϵVCA;m ðλ2 Þ C C: ϵVCA ¼ B B C ⋮ ⋮ ⋱ ⋮ @ A

ð7Þ

ϵVCA;1 ðλn Þ ϵVCA;2 ðλn Þ ⋯ ϵVCA;m ðλn Þ

Once the matrix ϵVCA is determined, distribution of the different components is obtained as described in the previous section, i.e., þ CU;sol ¼ ϵVCA Hexp :

ð8Þ

Two hypotheses are made in deriving the algorithm. The first one assumes that there is at least one pixel of the image containing a pure component, i.e., for every component j, there is a position k in the reconstructed image that fulfills

ϵVCA;j ¼ Hkexp ;

ð9Þ

where ϵVCA;j is the j-th column of the matrix ϵVCA . A maximum of n components can be obtained with the VCA algorithm, where the second hypothesis assumes that a linear combination of the absorption at different wavelengths is constant for every pixel in the image, i.e., for every pixel k of the image, the following condition is fulfilled: ∑ai Hðr k ; λi Þ ¼ 1;

ð10Þ

i

where ai are wavelength-dependent coefficients. Fig. 1 illustrates the principle of the VCA algorithm. Here it is considered that two absorbing components (with absorption spectra corresponding to ! ! vectors a and b ) are present and the sample is imaged at two different wavelengths. The black dots correspond to a representation of the absorbed energy at the imaged wavelengths for each pixel of the image. We further assumed that the first approximation above is verified, i.e., there are two pixels A and B whose ! ! absorption spectra correspond to the vectors a and b respectively. If the second assumption mentioned above is further accurate, the pixels lie in a hyperplane defined by Eq. (10), represented as a dashed line in Fig. 1a. In general, this hypothesis may turn inaccurate in cases where the light fluence variations are strong. The procedure to obtain the spectral signatures works as ! follows. First, a random direction is considered (m 1 in Fig. 1a). The vectors representing the absorption of the different pixels (black dots) are projected onto this direction while the pixel corresponding to the maximum value of the projection is considered (point A in Fig. 1a). The direction defined by this point ! (vector a in Fig. 1a) is assigned to the spectral signature of the first component. Subsequently, a random direction normal to the

121

spectra of the already determined components is taken, so that the new component is obtained by projecting the absorption of each pixel onto the new direction as described above. In the example shown in Fig. 1a, the spectral signature of the second ! ! component (vector b ) is retrieved by taking a direction m 2 ! normal to a . In Fig. 1a, the correct spectra of the components are retrieved. However, wrong spectra may be rendered in some cases. For example, the maximum value of the projection ! ! along direction m 1 in Fig. 1b corresponds to vector c , which does not correspond to the spectrum of an actual component. This kind of errors may appear in case that two or more components with a high concentration are located in the same region or in case that areas with a high light fluence (hot spots) are present. 2.3. Constrained vertex component analysis In many practical applications, the purpose of the unmixing in MSOT is to determine the distribution of a specific substance with known absorption spectrum. In particular, in molecular imaging applications one often seeks after a specific bio-marker targeting a certain molecular pathway. In such a case, the spectral signature of the absorption of the substance of interest can be a priori measured (for example with an absorption spectrometer). This information can be then included in the VCA algorithm by taking advantage of the fact that the spectral signatures are determined separately. Thus, one can fix the first l components and use the VCA algorithm to only determine the spectral signatures of other l  m components, so that the matrix ϵCVCA is given by 0 1 ϵ1 ðλ1 Þ ⋯ ϵl ðλ1 Þ ϵVCA;l þ 1 ðλ1 Þ ⋯ ϵVCA;m ðλ1 Þ B C B ϵ1 ðλ2 Þ ⋯ ϵl ðλ2 Þ ϵVCA;l þ 1 ðλ2 Þ ⋯ ϵVCA;m ðλ2 Þ C C: ð11Þ ϵCVCA ¼ B B ⋮ C ⋮ ⋮ ⋮ ⋱ ⋮ @ A ϵ1 ðλn Þ ⋯ ϵl ðλn Þ ϵVCA;l þ 1 ðλn Þ ⋯ ϵVCA;m ðλn Þ This variation of the suggested algorithm, which is herein termed constrained VCA, may help to improve the outcome in cases where errors are produced due to the two simplified assumptions mentioned above.

3. Materials and methods Numerical simulations and experimental measurements were done in order to test the performance of the VCA and the constrained VCA algorithms in MSOT imaging. A numerical phantom was utilized in order to simulate theoretical improvement of constrained VCA with respect to the standard version of the algorithm. The phantom was constructed using

Fig. 1. Principle of the vertex component analysis algorithm.

122

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

1 SO 2 0

1

0

( ) (A .U .)

( ) (A .U .)

1

700

800 (nm)

900

1

700

800 (nm)

900

800 (nm)

900

1 ( ) (A .U .)

( ) (A .U .)

0

0

700

800 (nm)

900

0

700

Fig. 2. (a) The simulated numerical phantom. (b) Random distribution of blood oxygenation levels as applied to the simulated data. The crosses and diagonal crosses in (b) show the two scenarios for placement of the AlexaFluor 750 absorbers. (c) and (d) Unmixing results attained with the standard version of the vertex component analysis algorithm. (e) and (f) Unmixing results for the constrained version of the algorithm.

cryosection image of a mouse in the kidney region [27], which mimics realistic tissue morphology to be imaged by optoacoustic tomography. For this, it was considered that the blood concentration is given by the image in Fig. 2a while a random map was taken to represent blood oxygenation levels (Fig. 2b). Also, two insertions containing AlexaFluor 750 (AF750) fluorochrome were included. In this way, it was assumed that absorption of light arises solely from the linear combination of optical absorption by oxygenated hemoglobin, deoxygenated hemoglobin and AF750. A uniform light fluence was considered throughout the specimen, so that optoacoustic image is directly represented by the distribution of the optical absorption coefficient. Thereby, no wavelengthdependent light fluence was taken into account, which would probably generate further distortion to the spectral variations of the images. Random noise was added to the images. The peak absorption of AF750 was taken equivalent to the absorption of blood at 797 nm, so that the optical absorptions of the three chromophores are in the same order. Two different scenarios were considered. In the first scenario, the AF750 insertions were located in two regions with a relatively high blood concentration but

different oxygenation levels (diagonal crosses in Fig. 2b). In the second scenario, the insertions were in regions with low background absorption (crosses Fig. 2b). Since the first assumption of the VCA algorithm is only valid when the background blood concentration at the location of the insertion is low, the performance of the algorithm is expected to be better for the first scenario. Experimental measurements on a post-mortem mouse were also performed with two polyethylene tubings inserted subcutaneously. The tubings contained different chromophoric (fluorescent) substances, namely AF750 and indocyanine green (ICG) dyes. The tubing material is transparent and does not introduce multispectral noise. Furthermore, the inner and outer diameters are 0.6 mm and 1 mm respectively, so that no significant distortion in the images is produced due to acoustic mismatch. The concentration of the two absorbers was varied so that the peak absorption coefficient ranged between μa ¼ 2 cm  1 and μa ¼ 0:2 cm  1 as was also verified with an absorption spectrophotometer. The overall purpose was to experimentally compare performance of the VCA versus PCA–ICA algorithms in terms of sensitivity. Also, the absorption spectra of the two absorbers, as measured with the spectrophotometer, were used in the constrained VCA algorithm in order to analyze the potential improvement in experimental measurements. A whole-body optoacoustic scanner was used to perform the experimental measurements [15]. The scanner consists of a tunable short-pulsed (nanosecond) laser, in which a fiber bundle was used in order to create a ring-type illumination around the imaged crosssection of the mouse. The laser wavelength was varied between 700 nm and 900 nm in a 10 nm step to acquire multispectral data. The generated ultrasonic waves were detected with a phased-array transducer consisting of 64 elements covering an angle of 1721 around the imaged cross-section. Each element of the array is cylindrically focused with a focal length of 40 mm and a central frequency of 5 MHz. The signals acquired for each wavelength were averaged 20 times and a band-pass filter with cut-off frequencies 0.1 MHz and 7 MHz was applied to the recorded signals. The subsequent image reconstruction was performed with a modelbased algorithm described in [28] and the images normalized with the measured per-pulse energy of the laser.

4. Results The simulation results for the two different scenarios presented in Fig. 2b are shown in Fig. 2c,e and d,f. The images and spectra rendered with the VCA algorithm by considering 3 components are displayed in Fig. 2c and d. It is shown that a distorted spectrum is retrieved if the background blood absorption is high (first column). The unmixed image in that case is also affected by cross-talk noise. On the other hand, a more faithful spectrum and overall a better quality image are rendered when the background blood concentration is low. The corresponding results obtained with the constrained VCA algorithm by including the spectral signature of AF750 are displayed in Fig. 2e and f. The cross-talk noise present in Fig. 2c is reduced, which indicates a better performance of the constrained version of the VCA. It is important to take into account that the results obtained with both algorithms depend on directions along which the projections are done, which are generally random. As a result, the representative component may or may not be obtained, yet the constrained VCA algorithm generally provides better results in numerical simulations. The single-wavelength tomographic reconstructions for the mouse experiment described in Section 3 are displayed in Fig. 3 overlaid with the full spectra of the two probes employed. The locations of the tubings are indicated with white arrows in the

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

figure. The optoacoustic images were acquired at wavelengths corresponding to the peak absorption of the probes, for which the absorption coefficient was μa ¼ 2 cm  1 for both probes, as measured with a spectrometer. Fig. 4 shows comparison of unmixing results obtained with the combined PCA–ICA (first and second

1

0

4mm

( ) (A .U .)

( ) (A .U .)

1

700

800

(nm)

0

900

4mm

700

800

900

(nm)

Fig. 3. Optoacoustic tomographic images at 750 nm (a) and 800 nm (b). The absorption spectra of AF750 and ICG measured with a spectrometer are also overlaid. The locations of the tubings containing AF750 and ICG are indicated with white arrows in (a) and (b) respectively.

123

column) versus the VCA algorithm (third and fourth column). In this case, four different concentrations of the probe were used so that each row corresponds to a different value of optical absorption in the tubings. The peak absorptions measured with a spectrometer are indicated within the figure. Both algorithms depend on random initializations and on the number of components assumed, so that representative components cannot always be efficiently retrieved. The results shown in Fig. 4 correspond to the most representative solutions obtained in both cases. For each unmixed image, the signal-to-noise ratio (SNR) is estimated as the maximum value in the region marked with a white square in the first row of Fig. 4 divided by the standard deviation outside this region. For the first column the resulting values of the SNR are 11.5, 6.0, 3.9 and 3.5, corresponding to absorption coefficients values of 2 cm  1, 1 cm  1, 0.4 cm  1 and 0.2 cm  1 respectively. The equivalent values are 8.8, 9.2, 11.1 and 9.3 for the second column, 9.2, 7.9, 5.2 and 3.8 for the third column and 9.1, 10.6, 9.0 and 9.7 for the fourth column. Generally, both algorithms present a similar performance in terms of sensitivity, whereas the computational effort for the VCA algorithm is significantly lower. Specifically, MATLAB implementations of the VCA and the PCA–ICA algorithms for 21 images with 200  200 pixels take approximately 0.2 and 1 s respectively on a Intel Core2 Duo CPU E8400 computer operating

Fig. 4. Unmixing results for different concentrations of AF750 and ICG with PCA–ICA (columns one and two) and with VCA (columns three and four). Each row corresponds to a peak absorption of the probe indicated in the figure.

124

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

1

800

(nm)

800

(nm)

900

0 700

800

900

(nm)

1

( ) (A.U.)

( ) (A.U.)

( ) (A.U.)

800

0 700

900

(nm)

1

1

0 700

0 700

900

( ) (A.U.)

( ) (A.U.) 0 700

1

( ) (A.U.)

1

800

(nm)

900

0 700

800

(nm)

900

Fig. 5. Unmixing results obtained with the constrained VCA by considering 6 components in which the first 2 spectral signatures (AF750 in (a) and ICG in (b)) are a priori fixed.

at 3.00 GHz with 8 GB of RAM. It is also important to emphasize that the VCA algorithm can be readily parallelized so that a more efficient implementation of the code can further substantially accelerate the unmixing procedure. Finally, performance of the constrained VCA algorithm was tested experimentally. As a representative example, Fig. 5 shows the results obtained with constrained VCA by assuming 6 components with the two first spectra (AF750 and ICG) being a priori known and used in the unmixing procedure. The absorption peak of the probes measured with a spectrometer is μa ¼ 2 cm  1 in this case. It is expected that the constrained VCA algorithm retrieves the distribution of AF750 and ICG for the first two components, where the spectra are a priori fixed. However, the distribution of these two absorbers is retrieved for the sixth and fourth components respectively, i.e., the spectral variation of the reconstructed absorbed energy fits better to the blindly-determined spectra than to the theoretical ones and the constrained approach fails to retrieve the correct results. This limits the capability of the constrained VCA algorithm to unmix probes with low concentrations in practical cases.

5. Discussion and conclusions Derivation of efficient unmixing procedures for sensitive detection of tissue bio-markers in the presence of strongly absorbing background is a challenging yet an important task for the development of efficient molecular imaging applications by means of multispectral optoacoustic tomography. In this work, applicability of the vertex component analysis algorithm in unmixing multispectral optoacoustic reconstructions has been evaluated. Specifically, the performance of vertex component analysis in terms of sensitivity was compared to another blind unmixing procedure based on a combination of principal component analysis and independent component analysis. The possibility of using a priori known spectra during the unmixing procedure was also suggested in a variation of the algorithm termed constrained vertex component analysis. The performance of this suggested semi-blind approach was also evaluated in numerical simulations and experimental measurements. The experimental results indicate that the overall sensitivity of the VCA method in optoacoustic unmixing is similar to that of the combination of PCA–ICA, whereas its computational time is substantially lower, an important advantage when dealing with

large amounts of data (e.g. in three-dimensional imaging) or trying to perform unmixing in real time [21,29]. A more accurate comparison between the sensitivity of the methods would involve considering different regions in mice, different probes and different relative concentrations of the probe with respect to blood background, as well as a statistical analysis on the interpretation of the results by independent users. We have also shown that, whereas constrained VCA performs better in numerical simulations, it turned less sensitive in experimental measurements where it was found that the spectral profile of the reconstructed absorbed energy fits better to the spectra of the components which were determined blindly than to the absorption spectrum measured with the spectrophotometer. This result might be attributed to the fact that the wavelength-dependent beam profile of the laser and errors in the reconstructions have also affected spectral signatures of probes located close to the surface of the object. Another possible reason is the fact that some chromophores can present an intrinsic optoacoustic spectrum differing substantially from the absorption spectrum measured with a spectrophotometer [30]. Moreover, the spectral dependence of the light attenuation further affects the spectral signatures of absorbers located deeper in tissue while fluctuations in the laser energy or otherwise motion artefacts may also contribute to the spectral variations in the reconstructions. Accurate modeling of the wavelength-dependent optical attenuation is then an important issue to address in future work. As was previously shown in [19], the PCA–ICA approach can generally achieve higher sensitivity compared to spectral fitting, and the results presented here confirmed the better performance of blind unmixing. Reconstructions obtained with the constrained VCA are nevertheless always easier to interpret and reproduce than those yielded with the standard version of the algorithm, for which the correct component may not be retrieved correctly depending on random initial settings. The constrained approach may also be beneficial when the probe is mixed with other tissue chromophores e.g. haemoglobin, which is common in blood-pool agents. However, since the spectrum measured with the spectrophotometer may not be reliable, an important next step is the accurate estimation of the intrinsic optoacoustic spectra of photoabsorbing probes with a dedicated optoacoustic spectrometer. The determination of the intrinsic optoacoustic spectra may also improve the sensitivity of the constrained VCA algorithm. In conclusion, the vertex component analysis algorithm can be employed in multispectral optoacoustics to successfully unmix the distribution of relatively low concentrations of spectrally-distinct bio-markers and optical probes. The fast computational time of the algorithm anticipates its applicability for real-time unmixing and handling of large multi-spectral datasets.

Acknowledgments D.R. acknowledges support from the European Research Council under Starting Independent Researcher Grant ERC-SG-LS7 (Dynamit). V.N. acknowledges support from European Union project FAMOS (FP7 ICT, Contract 317744). The authors thank Sarah Glasl, Uwe Klemm and Florian Jurgeleit for excellent technical assistance. References [1] Wang LV, Hu S. Photoacoustic tomography: in vivo imaging from organelles to organs. Science 2012;335(6075):1458–62. [2] Beard P. Biomedical photoacoustic imaging. Interface Focus 2011;1(4):602–31. [3] Brecht HP, Su R, Fronheiser M, Ermilov SA, Conjusteau A, Oraevsky AA. Wholebody three-dimensional optoacoustic tomography system for small animals. J Biomed Opt 2009;14(6):064007.

X. Luís Deán-Ben et al. / Optics and Lasers in Engineering 58 (2014) 119–125

[4] Ma R, Distel M, Deán-Ben XL, Ntziachristos V, Razansky D. Non-invasive whole-body imaging of adult zebrafish with optoacoustic tomography. Phys Med Biol 2012;57(22):7227–37. [5] Zhang HF, Maslov K, Stoica G, Wang LHV. Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging. Nat Biotechnol 2006;24 (7):848–51. [6] Razansky D, Distel M, Vinegoni C, Ma R, Perrimon N, Koster RW, et al. Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo. Nat Photon 2009;3(7):412–7. [7] Wang L, Maslov K, Wang LV. Single-cell label-free photoacoustic flowoxigraphy in vivo. Proc Natl Acad Sci 2013;110(15):5759–64. [8] Deán-Ben XL, Razansky D. Adding fifth dimension to optoacoustic imaging: volumetric time-resolved spectrally enriched tomography. Light: Science & Applications 2014;3:e137 http://dx.doi.org/10.1038/lsa.2014.18. [9] Wang K, Su R, Oraevsky AA, Anastasio MA. Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography. Phys Med Biol 2012;57(17):5399–423. [10] Queiros D, Deán-Ben XL, Buehler A, Razansky D, Rosenthal A, Ntziachristos V. Modeling the shape of cylindrically focused transducers in three-dimensional optoacoustic tomography. J Biomed Opt 2013;18(7):076014. [11] Deán-Ben XL, Razansky D, Ntziachristos V. The effects of acoustic attenuation in optoacoustic signals. Phys Med Biol 2011;56(18):6129–48. [12] Treeby BE. Acoustic attenuation compensation in photoacoustic tomography using time-variant filtering. J Biomed Opt 2013;18(3):036008. [13] Jetzfellner T, Rosenthal A, Buehler A, Dima A, Englmeier KH, Ntziachristos V, et al. Optoacoustic tomography with varying illumination and non-uniform detection patterns. J Opt Soc Am A 2010;27(11):2488–95. [14] Cox B, Laufer JG, Arridge SR, Beard PC. Quantitative spectroscopic photoacoustic imaging: a review. J Biomed Opt 2012;17(6):061202. [15] Buehler A, Herzog E, Razansky D, Ntziachristos V. Video rate optoacoustic tomography of mouse kidney perfusion. Opt Lett 2010;35(14):2475–7. [16] Farkas DL, Du CW, Fisher GW, Lau C, Niu W, Wachman ES, et al. Non-invasive image acquisition and advanced processing in optical bioimaging. Comput Med Imaging Graphics 1998;22(2):89–102. [17] Razansky D, Vinegoni C, Ntziachristos V. Multispectral photoacoustic imaging of fluorochromes in small animals. Opt Lett 2007;32(19):2891–3.

125

[18] Keshava N. A survey of spectral unmixing algorithms. Linc Lab J 2003;14(1): 55–78. [19] Glatz J, Deliolanis NC, Buehler A, Razansky D, Ntziachristos V. Blind source unmixing in multi-spectral optoacoustic tomography. Opt Express 2011;19(4): 3175–84. [20] Gamelin J, Maurudis A, Aguirre A, Huang F, Guo P, Wang LV, et al. A real-time photoacoustic tomography system for small animals. Opt Express 2009;17(13): 10489–98. [21] Buehler A, Deán-Ben XL, Claussen J, Ntziachristos V, Razansky D. Threedimensional optoacoustic tomography at video rate. Opt Express 2012;20(20): 22712–9. [22] Xiang L, Wang B, Ji L, Jiang H. 4-d photoacoustic tomography. Sci Rep 2013;3(1):1113. [23] Deán-Ben XL, Razansky D. Portable spherical array probe for volumetric real-time optoacoustic imaging at centimeter-scale depths. Opt Express 2013;21(23): 28062–71. [24] Nascimiento JMP, Dias JMB. Vertex component analysis: a fast algorithm to unmix hyperspectral data. IEEE Trans Geosci Remote Sens 2005;43(4): 898–910. [25] Shao P, Cox B, Zemp RJ. Estimating optical absorption, scattering and Grueneisen distributions with multiple-illumination photoacoustic tomography. Appl Opt 2011;50(19):3145–54. [26] Ripoll J, Ntziachristos V. Quantitative point source photoacoustic inversion formulas for scattering and absorbing media. Phys Rev E 2005;71(3):031912. [27] Sarantopoulos A, Themelis G, Ntziachristos V. Imaging the bio-distribution of fluorescent probes using multispectral epi-illumination cryoslicing imaging. Mol Imaging Biol 2011;13(5):874–85. [28] Deán-Ben XL, Ntziachristos V, Razansky D. Acceleration of optoacoustic model-based reconstruction using angular image discretization. IEEE Trans Med Imaging 2012;31(5):1154–62. [29] Deán-Ben XL, Ozbek A, Razansky D. Volumetric real-time tracking of peripheral human vasculature with GPU-accelerated three-dimensional optoacoustic tomography. IEEE Trans Med Imaging 2013;32(11):2050–5. [30] Laufer J, Jathoul A, Pule M, Beard P. In vitro characterization of genetically expressed absorbing proteins using photoacoustic spectroscopy. Biomed Opt Express 2013;4(11):2477–90.