Multispectral optoacoustic (photoacoustic) tomography (MSOT) is a hybrid modality that can image through several millimeters to centimeters of diffuse tissues, attaining resolutions typical of ultrasound imaging. The method can further identify tissue biomarkers by decomposing the spectral contributions of different photo-absorbing molecules of interest. In this work we investigate the performance of blind source unmixing methods and spectral fitting approaches in decomposing the contributions of fluorescent dyes from the tissue background, based on MSOT measurements in mice. We find blind unmixing as a promising method for accurate MSOT decomposition, suitable also for spectral unmixing in fluorescence imaging. We further demonstrate its capacity with temporal unmixing on real-time MSOT data obtained in-vivo for enhancing the visualization of absorber agent flow in the mouse vascular system.
© 2011 Optical Society of America
Optoacoustic (or photoacoustic) imaging is a potent imaging method that can offer high resolution visualization of optical absorption deep in living tissues [1, 2, 3]. Being a hybrid technique, it combines the advantages of visualizing highly specific optical absorption contrast with the high spatial resolution of ultrasound. Recording optoacoustic signals emitted by tissues from multiple projections, in response to the absorption of ultrafast light pulses by tissue elements, allows obtaining tomographic reconstructions of light absorbing structures in deep tissues at resolutions of tens of micrometers and below. The potential of this technology has been demonstrated in many applications, in particular imaging of blood vessels , blood flow [5, 6] or using various photo-absorbing agents and nanoparticles . Recently, the use of acoustic detector arrays has eliminated the need for mechanical scanning of a single detector and enabled video rate optoacoustic imaging [8, 9, 10], allowing a wide number of possible applications in both clinical research and drug discovery. While subsurface imaging at the microscopic or macroscopic regimes has been suggested since the early 80’s [11, 12], optoacoustic imaging has significantly improved in recent years. Of particular importance are spectral imaging methods where tissues are excited at two or more wavelengths and generate in response images of the bio-distribution of various tissue molecules or biomarkers, for example in visualizing oxygen saturation , extrinsically administered fluorochromes , fluorescent proteins  or circulating nanoparticles [16, 17].
In many cases, the contributions of various photo-absorbing agents, such as contrast agents, targeted probes or nanoparticles may constitute only small signal variances over the background absorption, depending on their concentration. This may complicate the detection of the photoabsorbing agent on single wavelength images, even if the selected wavelength corresponds to the absorption maximum of the agent of interest. When data at multiple wavelengths are obtained however, it is possible to improve the contrast and detection sensitivity by resolving the spectral signature of the absorber agent used over other non-specific spectral contributions, e. g. from highly absorbing hemoglobin.
Unmixing methods based on differential or fitting algorithms use the known spectral or temporal information to process the image on a pixel-by-pixel basis. This separation can be challenging since the exact spectral profile of the background contribution is not always known in in-vivo tissue imaging. In addition, the spectral signature of the agent of interest may also be not accurately known, for instance the absorption spectrum may change in different biochemical environments. Moreover, light attenuation and ultrasonic dispersion in tissues leads to a corresponding non-linear relationship between the measured optoacoustic signals and the corresponding target concentration, as a function of depth, or target size . This is particularly evident for volumetric tissue imaging since the optoacoustic signals depend on the product of the local light fluence at different depths times the local absorption coefficient of the agent of interest and other background chromophores . Finally, the spectral profile of light propagating in tissue is also altered by depth, since different wavelengths are attenuated by a different rate, which may also contribute to intensity variations in the recorded MSOT signals with depth.
These potential complications in spectral unmixing may be addressed with the use of multivariate data analysis and matrix factorization algorithms, such as principal component analysis (PCA) , non-negative matrix factorization (NNMF) , multivariate curve resolution (MCR) , or independent component analysis (ICA) . Applications of these methods include astronomy , pattern recognition [25, 26], anatomical segmentation , fluorescence imaging [28, 29], and finance . Herein we have investigated PCA and ICA in spectral unmixing of optoacoustic images and compared them with spectral fitting schemes. Contrary to the pixel-by-pixel processing approach in the differential and fitting unmixing methods, the key element in the multivariate approaches is the unaided identification of changes that is common across various pixels, helping to identify contrast agents that have a non-uniform spatial bio-distribution. Of particular interest in the context of unmixing MSOT measurements was the behavior of these algorithms concerning the dependence of the optoacoustic signals on depth and other tissue optical parameters, as well as biologically relevant spectral shifts of the detected fields. We have examined the performance of the unmixing methods in spectral and temporal decomposition of tissue components in in-vivo and ex-vivo experiments.
Generally, the term “unmixing” refers to the unsupervised decomposition of mixed observations into usually a smaller number of endmembers . Its purpose is to improve the spectral or temporal resolution by identifying a small number of underlying sources that are superposed to yield the measured data and to separate them into individual images. For simplicity, in the following description we use the terminology of the spectral domain, but the methodology holds for the temporal domain as well.
2.1. Spectral Fitting
The most commonly employed unmixing method is spectral fitting, i. e. finding the source component that best fits a given absorption spectrum in the least-squares sense. Given the (n × m) multispectral measurement matrix M, where n is the number of image pixels and m is the number of measurements, as well as the (k × m) spectral matrix S with the absorption coefficients of the k components at the m measurement wavelengths, the data can be unmixed with the Moore-Penrose pseudoinverse [32, 33] S+.
2.2. Principal Component Analysis
Principal Component Analysis [34, 35] is a blind source unmixing technique, that is based on the assumption that the source components are statistically uncorrelated. PCA yields a linear orthogonal transformation into a new coordinate system, in which the largest data variance is projected onto the first principal component, the largest remaining variance onto the second one, and so on. Consequently, the correlated measurement data is unmixed by being transformed to uncorrelated principal components. PCA can be calculated as a singular value decomposition  of M or as an eigenvalue decomposition of its covariance matrix . In both cases, not only the unmixed components RPCA but also the corresponding transformation matrix UPCA are returned by the PCA algorithm, i. e.
2.3. Independent Component Analysis
Independent Component Analysis  (ICA) is yet another blind source separation technique, but it is based on a different assumption about the sources than PCA. While the latter assumes uncorrelated sources, ICA finds endmembers that satisfy the more general and therefore stronger condition of statistical independence. The algorithm seeks a transformation of the dependent mixed spectral components into a set of independent source components and also yields the corresponding mixing matrix UICA, i. e.37], a sum of non-gaussian variables is closer to a gaussian distribution than the original variables. Consequently, the mixed multispectral measurements are expected to be more gaussian than the unmixed independent components. Herein, we employ the FastICA algorithm , which finds the independent sources by using a fixed point iteration scheme in order to maximize non-Gaussianity, as measured by kurtosis, a fourth order statistical moment. By maximizing the kurtosis of the components, a non-Gaussian, and thus an independent representation is found. A key difference between ICA and PCA is that no eigenvalues and hence no measure of a component’s significance is obtained.
As a consequence, ICA requires some preprocessing when dealing with large datasets because it attempts to retrieve as many components as the number of available measurements and cannot sort out superfluous signals on its own. Thus, many independent components merely represent noise or minor signal variations, which makes the interpretation of the results more difficult. A logical solution would be using PCA as a preprocessing step. Thereby, the data are initially analyzed by PCA and the insignificant principal components, mainly containing noise, can be discarded. This step not only serves to eliminate the noise but also reduces the dimension of the data, presenting the ICA algorithm with a better conditioned problem. When R′PCA is the reduced PCA subset, the unmixing is denoted byEquation (5) also shows that the estimated spectral characteristics of the unmixed components can be retrieved as the product of the two mixing matrices and .
3. Experimental Methods
3.1. Experimental system and data processing
Experimental MSOT measurements were performed with a real-time MSOT scanner, previously described in Ref. . Briefly, a tunable (680 – 950 nm) pulsed (< 10 ns) optical parametric oscillator laser (Opotek Inc., Carlsbad, CA, USA) with a 10 Hz repetition rate is used for signal generation. The laser beam is guided into a silica fused-end fiber bundle (CeramOptec GmbH, Bonn, Germany) creating a ring shaped illumination pattern of ∼7 mm width upon the surface of the animal. Signal collection is based on a custom-made 64-element focused transducer array (Imasonic SaS, Voray, France) covering a solid angle of 172° around the imaged object, while the detection plane coincided with the center of the illumination ring. The individual detection elements are manufactured using piezocomposite technology with a central frequency of 5 MHz, a bandwidth (−6 dB) of more than 50% sensitivity of ∼ 18μV/Pa and are shaped to create a cylindrical focus at 40 mm. The detected signals are digitized at a sampling frequency of 60 MHz by 8 multi-channel analog to digital converters (Model PXI5105, National Instruments, Austin, TX, USA) with noise floor of . A linear stage (NRT150, Thorlabs GmbH, Karlsfeld, Germany) allows linear translation of the animal holder in the axial z direction. The data acquisition is synchronized so that the signals are acquired only when the stage comes to a complete rest. For all experiments the optoacoustic pressure distribution was reconstructed with a filtered backprojection algorithm .
3.2. Animal preparation and imaging
For spectral unmixing, a nude mouse cadaver phantom with two inclusions of fluorophores simulating tumors was prepared. The inclusions consisted of 1% agar solution and 1% intralipid (Sigma-Aldrich, St. Louis, MO, USA), which was initially heated to 90 °C. While cooling down to 50 °C, solutions of 12.9 μM of ICG (peak absorption ≈ 800 nm, Pulsion Medical, Munich, Germany) and 3.7 μM of Cy7 Cyanine Dye (peak absorption ≈ 750 nm, GE Healthcare, Little Chalfont, UK) were prepared to both have 2 cm−1 optical density at their absorption peaks. At 35 °C, 50 μL of each solution was injected subcutaneously in the upper back / neck area of a euthanized mouse and were left to equilibrate with room temperature and solidify to a shape that is similar to a subcutaneous tumor.
The mouse was positioned in supine position inside a water-impenetrable yet transparent membrane, providing a wide tomographic view of ∼ 180°. A total of 31 spectral measurements in the 700 – 850 nm range were acquired in steps of 5 nm and reconstructed with the filtered backprojection algorithm . Interwavelength variations have been normalized by powermeter readings (Model FieldMaxII-TOP, Coherent GmbH, Dieburg, Germany) for quantitative multispectral reconstructions.
For imaging and unmixing in the time domain, a white CD1 mouse was anesthetized with ketamine and xylazine mixture, catheterized in the right tail vein and injected with 50 μL of black india ink saline solution as contrast medium. An axial slice in the pelvic region, approximately 2 cm away from the catheter towards the heart, was imaged with MSOT for 80 s, using herein data at 1 Hz frame rate, after performing 10 signal averages at a wavelength of 800 nm. After imaging, the mouse was euthanized in situ.
4. Spectral domain unmixing
We showcase the described methods on the reconstructions from the absorber implantation experiment, where the presence of two contrast agents with an overlapping spectrum presents a challenging unmixing problem. Reconstructions of an axial slice in the neck area at 4 representative wavelengths are shown in Fig. 1
The two implantations can be easily identified in the left and the right side of the mouse neck under the skin. It can be observed that the reconstructed optoacoustic signal at the Cy7 implantation area has a peak at about 750 nm while at the ICG area the peak is at around 800 nm, which correlates well with the expected absorption peaks of the corresponding dyes. The ICG, Cy7 and background components, unmixed by fitting (pseudoinverse), PCA, and ICA are displayed in Fig. 2. The spectral absorption data of the two fluorochromes and deoxygenated hemoglobin, that were used for the fitting unmixing method, as well as the corresponding spectra calculated by ICA, are shown in Fig. 3. Deoxygenated hemoglobin was chosen as the background absorber since it is assumed to be the dominant absorber in the dead tissue.
The images show that all three methods successfully managed to isolate the two fluorochrome components from the background. However, fitting with the pseudoinverse algorithm yielded crosstalk of the background signal into the ICG and Cy7 components. Unmixing with PCA also exhibited crosstalk, which is attributed to the orthogonality criterion that is imposed. On the other hand, ICA accurately unmixed the components and the spectra calculated are very close to the expected spectra of the employed fluorochromes. The correspondence between the ICA spectra and the known absorption curves was determined by visual inspection and the ICA components were allocated accordingly. Using Pearson’s product moment correlation coefficient ρ, the correlation between the known spectra and the ones estimated by ICA was calculated to be ρ = 0.91 for the background, ρ = 0.83 for ICG and ρ = 0.96 for Cy7. In order to evaluate the performance of the unmixing methods the signal-to-background ratios (SBR) were calculated using the standard deviation of the pixel values for the two fluorochrome inclusions. The results, as shown in Table 1 coincide with the visual examination that the ICA unmixing method outperforms fitting and PCA.
The proposed preprocessing of ICA by PCA was also applied to this dataset. While it reduced the number of components and, thus, made the identification and allocation of the components easier, the unmixed images did not improve upon ICA alone, since the two results were virtually indistinguishable.
The results demonstrate the usefulness of blind source unmixing methods in the analysis of MSOT measurements. In particular, ICA appears more suited for MSOT imaging as it yields better separation of components than PCA or spectral fitting methods, exhibiting a significantly higher SBR for both ICG and Cy7.
5. Time domain unmixing
The imaging of circulation dynamics by MSOT in combination with time domain unmixing is demonstrated on the experimental data shown in Fig. 4. In the field of view, various anatomical structures can be seen, such as the bladder, spine, and the extensions of the right and the left tail veins in the torso (Fig. 4a). The right tail vein appears to be less absorbing because of the lower blood volume due to the restriction of flow from catheterization.
Images acquired at different time points are shown in Figs. 4b–d. ICG was injected at t = 10 s leading to a notably increasing and then gradually decreasing signal in the right vein. A few seconds later the signal in the left vein showed an oscillatory pattern that eventually faded out, indicating that the absorber was fully diluted in the blood stream after 50 s. To investigate the unmixing methods in the time domain, we examined the ability to separate the tissue components based on the dynamic changes of the images. The first three principle components from PCA on the full 80-frame time-series and the corresponding temporal profiles (eigenvectors) are shown in Figs. 5a–d.
The first principal component contains mainly the parts of the image that remain constant, the second one shows the most prominent changes in the images, whereas the third one shows the changes with respect to the two previous ones and so forth. Due to the orthogonality criterion, the PCA components do not disentangle the two types of veins and the temporal profiles also show negative values that make it difficult to identify a specific temporal pattern. Principal components from 8 to 80 mostly contain noise. On the other hand, performing ICA directly on the full set of 80 images yielded a large number of independent components (not shown), which were not found useful as they represented minor changes in the images, such as variations introduced by laser power fluctuations, mouse movement and noise. A third combinational method proved to be most viable, where PCA was employed as a preprocessing step for ICA. As stated above, the first eight PCA components practically contain all the useful information on the temporal variations. This reduced subset can subsequently be used as an input into ICA, feeding the algorithm with a better-conditioned dataset. After noise reduction achieved by PCA pre-processing, ICA yielded the components shown in Figs. 5e–g. The temporal profiles can be computed as a product of the PCA and ICA mixing matrices and are shown in Fig. 5f.
6. Discussion and Conclusion
With the emergence of real-time MSOT [10, 17], spectral decomposition of various tissue biomarkers and dynamic monitoring of internal or external photo-absorbing molecules has become possible in small animals. The post-processing of the generated spectral and temporal data requires however a robust and reliable method for decomposing the constitutive contributions from the measurements. In this paper, we have suggested and experimentally examined the blind source unmixing approach as a very promising alternative to separating mixed component data by fitting. Fitting requires the knowledge of the spectral or temporal profiles of all components to compute a generalized inverse of the mixing matrix. However, these are often not available for all sources with the necessary accuracy. The primary advantage of blind unmixing is that no a priori information is needed, making the technique suitable for a wide range of applications. It can separate contrast agents of a common spectral or temporal bio-distribution from background absorbers such as hemoglobin. The unmixing quality of PCA and ICA demonstrated herein on mouse measurements was found to be superior to spectral fitting, further generating components that can lead to spectral identification of specific tissue biomarkers such as hemoglobin or the contributions of intrinsically expressed or administered molecular probes. Additionally, the use of blind unmixing was successfully demonstrated in the time domain by separation of blood veins based on dynamic changes of an optoacoustic image sequence. It is evident, that the sensitivity of the imaging systems directly benefits from an improved unmixing performance.
Among the discussed blind source methods ICA appears to be the superior method, yielding a cleaner unmixing compared to PCA in the spectral domain. Contrary to spectral unmixing, it is unlikely to have an a priori knowledge of temporal profiles of biological processes or the contrast agent’s temporal biodistribution in living organisms. Therefore fitting procedures are not common in dynamic measurements. In response, we found that temporal unmixing based on blind decomposition methods required a combination of PCA and ICA to lead to clearly perceived time components.
Future research will aim to further evaluate the performance of blind source techniques in the unmixing of MSOT data and to define suitable pre- and post-processing steps that can help to automate the procedure and also guide the unmixing with a priori information that may be available. Applications of these unmixing methods are envisioned also in the field of optical fluorescence imaging, where they can help to increase detection sensitivity, and in autofluorescence removal.
The authors would like to thank Dr. George Themelis for the helpful discussions, and Dr. Niels Harlaar for the help in the surgical procedures. N. C. D. is supported by a Marie Curie Intra-European Fellowship within the 7th European Community Framework. D. R. acknowledges support from the German Research Foundation (DFG) Research Grant (RA 1848/1). V. N. acknowledges support from an ERC Advanced Investigator Award and BMBF’s Medizintechnik award.
References and links
2. M. Xu and L. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77, 041101 (2006). [CrossRef]
4. C. G. A. Hoelen, F. F. M. de Mul, R. Pongers, and A. Dekker, “Three-dimensional photoacoustic imaging of blood vessels in tissue,” Opt. Lett. 23, 648–650 (1998). [CrossRef]
6. P.-C. Li, S.-W. Huang, C.-W. Wei, Y.-C. Chiou, C.-D. Chen, and C.-R. C. Wang, “Photoacoustic flow measurements by use of laser-induced shape transitions of gold nanorods,” Opt. Lett. 30, 3341–3343 (2005). [CrossRef]
8. H.-P. Brecht, R. Su, M. Fronheiser, S. A. Ermilov, A. Conjusteau, and A. A. Oraevsky, “Whole-body three-dimensional optoacoustic tomography system for small animals,” J. Biomed. Opt. 14, 064007 (2009). [CrossRef]
9. J. Gamelin, A. Maurudis, A. Aguirre, F. Huang, P. Guo, L. V. Wang, and Q. Zhu, “A real-time photoacoustic tomography system for small animals,” Opt. Express 17, 10489–10498 (2009). [CrossRef]
11. G. Busse and A. Rosencwaig, “Subsurface imaging with photoacoustics,” Appl. Phys. Lett. 36, 815–816 (1980). [CrossRef]
12. A. Rosencwaig, “Potential clinical applications of photoacoustics,” Clin. Chem. 28, 1878–1881 (1982). [PubMed]
13. X. Wang, X. Xie, G. Ku, L. V. Wang, and G. Stoica, “Noninvasive imaging of hemoglobin concentration and oxygenation in the rat brain using high-resolution photoacoustic tomography,” J. Biomed. Opt. 11, 024015 (2006). [CrossRef]
15. D. Razansky, M. Distel, C. Vinegoni, R. Ma, M. Perrimon, R. W. Koster, and V. Ntziachristos, “Multispectral opto-acoustic tomography of deep-seated fluorescent proteins in vivo,” Nat. Photonics 3, 412–417 (2009). [CrossRef]
16. P.-C. Li, C.-R. C. Wang, D.-B. Shieh, C.-W. Wei, C.-K. Liao, C. Poe, S. Jhan, A.-A. Ding, and Y.-N. Wu, “In vivo photoacoustic molecular imaging withsimultaneous multiple selective targeting using antibody-conjugated gold nanorods,” Opt. Express 16, 18605–18615 (2008). [CrossRef]
17. A. Taruttis, E. Herzog, D. Razansky, and V. Ntziachristos, “Real-time imaging of cardiovascular dynamics and circulating gold nanorods with multispectral optoacoustic tomography,” Opt. Express 18, 19592–19602 (2010). [CrossRef]
19. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative optoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imaging 28, 1997–2006 (2009). [CrossRef]
20. I. T. Jolliffe, Principal Component Analysis, 2nd ed. (Springer, 2002).
21. A. Cichocki, R. Zdunek, A. H. Phan, and S. I. Amari, Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi-way Data Analysis and Blind Source Separation, 1st ed. (Wiley, 2009). [PubMed]
22. R. Tauler, B. Kowalski, and S. Fleming, “Multivariate curve resolution applied to spectral data from multiple runs of an industrial process,” Anal. Chem. 65, 2040–2047 (1993). [CrossRef]
23. A. Hyvärinen, J. Karhunen, and E. Oja, Independent Component Analysis, Adaptive and Learning Systems for Signal Processing, Communications, and Control, 1st ed. (Wiley InterScience, 2002).
25. B. A. Draper, K. Baek, M. S. Bartlett, and J. R. Beveridge, “Recognizing faces with pca and ica,” Comput. Vis. Image Underst. 91, 115 – 137 (2003). [CrossRef]
26. J. Yang, D. Zhang, A. F. Frangi, and J. Y. Yang, “Two-dimensional pca: a new approach to appearance-based face representation and recognition,” IEEE Trans. Pattern Anal. Mach. Intell. 26, 131–137 (2004). [CrossRef]
27. E. M. C. Hillman and A. Moore, “All optical anatomical co registration for molecular imaging of small animals using dynamic contrast,” Nat. Photonics 1(9) 526–530 (2007). [CrossRef]
28. H. Xu and B. W. Rice, “In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique,” J. Biomed. Opt. 14, 064011 (2009). [CrossRef]
29. A.-S. Montcuquet, L. Hervé, F. Navarro, J.-M. Dinten, and J. I. Mars, “Nonnegative matrix factorization: a blind spectra separation method for in vivo fluorescent optical imaging,” J. Biomed. Opt. 15, 056009 (2010). [CrossRef]
30. S. Clémençon and S. Slim, “On portfolio selection under extreme risk measure: the heavy-tailed ICA Model,” Int. J. Theor. Appl. Finance 10, 449–474 (2007). [CrossRef]
31. N. Keshava, “A survey of spectral unmixing algorithms,” Lincoln Lab. J. 14, 55–78 (2003).
32. E. Moore, “On the reciprocal of the general algebraic matrix,” Bull. Am. Math. Soc. 26, 394–395 (1920).
33. R. Penrose, “A generalized inverse for matrices,” in Proceedings of the Cambridge Philosophical Society , (1955), Vol. 51, pp. 406–412. [CrossRef]
34. K. Pearson, “On lines and planes of closest fit to a system of points in space,” London Edinburgh Dublin Philos. Mag. J. Sci.6, 559–572 (1901).
35. S. M. Kay, Fundamentals of Statistical Signal Processing, 1st ed. (Prentice Hall PTR, 1993), Vol. 1.
36. J. Nash, “The singular-value decomposition and its use to solve least-squares problems,” in Compact Numerical Methods for Computers: Linear Algebra and Function Minimization, 2nd ed. (Inst. of Physics Pub., 1990), pp. 30–48.
37. L. Le Cam, “The central limit theorem around 1935,” Stat. Sci. 1, 78–91 (1986). [CrossRef]
38. A. Hyvrinen, “Fast and robust fixed-point algorithms for independent component analysis,” IEEE Trans. Neural Netw. 10, 626–634 (1999). [CrossRef]
39. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71, 016706 (2005). [CrossRef]