Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues

Light propagating in tissue attains a spectrum that varies with location due to wavelength-dependent fluence attenuation, an effect that causes spectral corruption. Spectral corruption has limited the quantification accuracy of optical and optoacoustic spectroscopic methods, and impeded the goal of imaging blood oxygen saturation (sO2) deep in tissues; a critical goal for the assessment of oxygenation in physiological processes and disease. Here we describe light fluence in the spectral domain and introduce eigenspectra multispectral optoacoustic tomography (eMSOT) to account for wavelength-dependent light attenuation, and estimate blood sO2 within deep tissue. We validate eMSOT in simulations, phantoms and animal measurements and spatially resolve sO2 in muscle and tumours, validating our measurements with histology data. eMSOT shows substantial sO2 accuracy enhancement over previous optoacoustic methods, potentially serving as a valuable tool for imaging tissue pathophysiology.

T he assessment of tissue oxygenation is crucial for understanding tissue physiology and characterizing a multitude of conditions, including cardiovascular disease, diabetes, cancer hypoxia 1 or metabolism. Today, measurements of the partial pressure of oxygen in tissue (pO 2 ) and hypoxia measurements remain challenging and often rely on invasive methods that may change the tissue physiology, such as singlepoint needle polarography or immunohistochemistry 2 . Noninvasive imaging methods have been also considered, underscoring the importance of assessing pO 2 , but come with limitations. Positron emission tomography or single-photon emission computed tomography assess cell hypoxia through the administration of radioactive tracers 2 , but are often not well suited for quantifying pO 2 , suffer from low-spatial resolution and are unable to provide longitudinal or dynamic-imaging capabilities. Electron paramagnetic resonance imaging 3 can measure tissue pO 2, but is not widely used and offers limited spatial and temporal resolution. Imaging methods using tracers may be further limited by restricted tracer bio-distribution, in particular to hypoxic areas. Tracer-free modalities have also been researched, in particular blood-oxygen-level dependent MRI 4 , which however primarily assesses only deoxygenated haemoglobin and, therefore, presents challenges in quantifying oxygenation and blood volume 5 .
Measurement of blood oxygen saturation levels (sO 2 ) is a vital tissue physiology measurement and can provide an alternative way to infer pO 2 and hypoxia. Arterial sO 2 is widely assessed by the pulse oximeter, but this technology cannot be applied to measurements other than arterial blood. Optical microscopy methods like phosphorescence quenching microscopy 6 or optoacoustic (photoacoustic) microscopy 7 can visualize oxygenation in blood vessels and capillaries but are restricted to superficial (o1 mm depth) measurements. Diffuse optical methods received significant attention in the last two decades for sensing and imaging oxy-and deoxygenated haemoglobin deeper in tissue 8 . Despite recent progress 9 , diffuse optical methods are inherently limited in spatial resolution and accuracy, due to photon scattering. Owing to the high heterogeneity of blood sO 2 in tissue, the values reported by diffuse optical methods are often hard to interpret.
Multispectral optoacoustic tomography (MSOT) detects the spectra of oxygenated and deoxygenated haemoglobin in high resolution deep within tissues, since signal detection and image reconstruction are not significantly affected by photon scattering 10,11 . Despite the principal MSOT suitability for noninvasive imaging of blood oxygenation, accuracy remains limited by the dependence of light fluence on depth and light colour. Unless explicitly accounted for, the wavelength-dependent light fluence attenuation with depth alters the spectral features detected and results in inaccurate estimates of blood sO 2 12,13 . Despite at least two decades of research in optical imaging, the problem of light fluence correction has not been conclusively solved. To date, this problem has been primarily studied from an optical property quantification point of view 13,14 . However, it is not possible today to accurately image tissue optical properties in vivo, in high resolution, and compute light fluence 13 . Therefore, quantitative sO 2 measurement deep in tissue in vivo remains an unmet challenge. Conventional spectral optoacoustic methods 15,16 typically ignore the effects of light fluence and employ linear spectral fitting with the spectra of oxy-and deoxy-haemoglobin for estimating sO 2 (linear unmixing), a common simplification that can introduce substantial errors in deep tissue.
In this work, we found that the spectral patterns of light fluence expected within the tissue can be modelled as an affine function of a few reference base spectra, independently of the specific distribution of tissue optical properties or the depth of the observation. We show how this principle can be employed to solve the spectral corruption problem without knowledge of the tissue optical properties, and significantly increase the accuracy of spectral optoacoustic methods. The proposed method, termed eigenspectra Multispectral Optoacoustic Tomography (eMSOT), can provide quantitative estimation of blood sO 2 in deep tissue. We demonstrate the superior performance of the method with 42,000 simulations, phantom measurements and in vivo controlled experiments. Then, using eMSOT, we image oxygen gradients in the skeletal muscle in vivo, previously only accessible through invasive methods. Furthermore, we show the application of eMSOT in quantifying blood oxygenation gradients in tumours during tumour growth or O 2 challenge, and relate label-free noninvasive eMSOT readings to tumour hypoxia; demonstrating the ability to measure quantitatively the perfusion hypoxia level in tumours, as confirmed with invasive histological gold standards.

Results
eMSOT concept and application. A new concept of treating light fluence in diffusive media/tissues is introduced, based on the observation that the light fluence spectrum at different locations in tissue is the result of a cumulative light absorption operation by tissue chromophores, such as haemoglobin. We, therefore, hypothesized that there exists a small number of base spectra that can be combined to predict any fluence spectrum present in tissue; therefore, avoiding the unattainable task of knowing the distribution of tissue optical properties at high resolution. To prove this hypothesis, we first applied principal component analysis (PCA) on 1,470 light fluence spectral patterns, which were computed by simulating light propagation in tissues at 21 different (uniform) oxygenation states of haemoglobin and 70 different discrete depths (Methods). PCA analysis yielded four significant base spectra, that is, a mean light fluence spectrum (Fig. 1a) and three fluence Eigenspectra (Fig. 1b-d). PCA was used due to its optimality in modelling the spectral variability of light fluence in a linear manner (see Methods). We found that the selection of three Eigenspectra offers a simple model with relatively high-modelling accuracy (Fig. 1e).
We then postulated that light fluence spectra in arbitrary and non-uniform tissues can be modelled as a superposition of the mean fluence spectrum (F M ) and the three Eigenspectra (F i (l), i ¼ 1..3) multiplied by appropriate scalars m 1 , m 2 and m 3 , termed Eigenfluence parameters. To validate this hypothesis we computed the light fluence in 4500 simulated tissue structures of different and non-uniform optical properties and haemoglobin oxygenation values (Supplementary Note 1). For each pixel, we fitted the simulated light fluence spectrum to the Eigenspectra model and derived the Eigenfluence parameters (m 1 , m 2 , m 3 ) and a fitting residual value. The residual value represents the error of the Eigenspectra model in matching the simulated data and typically assumed values below 1% (see Supplementary Note 1, Supplementary Fig. 1), indicating that three Eigenspectra can accurately model all simulated fluence spectra generated in tissues of arbitrary structure. We further observed that the values of m 2 vary primarily with tissue depth while the values of m 1 , m 3 also depend on the average levels of background blood sO 2 (see Fig. 1fh). Intuitively, this indicates that the second Eigenspectrum F 2 (l) is mainly associated with the modifications of the light fluence spectrum due to depth and the average optical properties of the surrounding tissue, while the first Eigenspectrum F 1 (l) is also associated with the 'spectral shape' of light fluence that relates to the average oxygenation of the surrounding tissue. Localized measurements of light fluence spectra obtained in vivo and post mortem corroborated these observations ( Supplementary Fig. 2).
Following these observations, we propose eigenspectra MSOT (eMSOT), based on three eigenspectra F 1 (l), F 2 (l), F 3 (l), as a method that formulates the blood sO 2 estimation problem as a nonlinear spectral unmixing problem (see Methods), that is: where P(r,l) is the multispectral optoacoustic image intensity obtained at a position r and wavelength l, e HbO 2 ðlÞ and e Hb (l) are the wavelength-dependent molar extinction coefficients of oxygenated and deoxygenated haemoglobin, c 0 HbO 2 r ð Þ and c 0 Hb (r) are the relative concentrations of oxygenated and deoxygenated haemoglobin (proportional to the actual ones with regard to a common scaling factor, see Methods), and F 0 (r,l) ¼ F M (l) þ m 1 (r)F 1 (l) þ m 2 (r)F 2 (l) þ m 3 (r)F 3 (l). Equation (1) defines a nonlinear inversion problem, requiring measurements at 5 wavelengths or more for recovering the five unknowns, that is, c 0 HbO2 r ð Þ, c 0 Hb (r), m 1 (r), m 2 (r), m 3 (r) and is solved as a constrained optimization problem (see Methods, Supplementary Note 2, Supplementary Fig. 3). Since the light fluence varies smoothly in tissue, we only compute the Eigenfluence parameters on a coarse grid subsampling the region of interest (Fig. 1i), for computational efficiency. Then, cubic interpolation is employed to compute the Eigenfluence parameters in each pixel within the convex hull of the grid (Fig. 1j) and calculate a fluence spectrum F 0 (r,l) for each pixel using these parameters. Equation (1) is then solved for c 0 HbO 2 r ð Þ and c 0 Hb (r), and sO 2 is computed (see Methods).
The performance of eMSOT was validated using simulated data obtained from a light propagation model (finite element solution of the diffusion equation) applied on 42,000 randomly created maps of different optical properties, simulating different tissue physiological states (Supplementary Note 3). Figure 1k depicts a representative example of a simulated blood sO 2 map and visually showcases the differences between the eMSOT sO 2 image (middle), the sO 2 image obtained using linear unmixing (left), and the original sO 2 simulated image (right). eMSOT offered significantly lower sO 2 estimation error with depth, compared with the linear fitting method (Fig. 1l). A substantially improved sO 2 estimation accuracy was observed using eMSOT over linear unmixing when we analysed the complete simulation data set (Fig. 1m). In particular, for imaging tissue depths of 45 mm eMSOT offered a 3-8-fold sO 2 estimation improvement over linear unmixing ( Supplementary Fig. 4n). A thorough validation of eMSOT performance across different data sets, optical properties and grid densities is presented in For experimentally assessing the accuracy of eMSOT, we performed a series of blood phantom experiments that suggest an up to 10-fold more reliable sO 2 estimation derived by eMSOT, as compared with conventional linear unmixing (Supplementary Note 4, Supplementary Fig. 5). In addition, controlled mouse measurements (n ¼ 4) were performed in vivo, under gas anaesthesia, by rectally inserting capillary tubes containing blood at 100 and 0% sO 2 levels (Methods). The mice were imaged in the lower abdominal area under 100% O 2 and 20% O 2 breathing conditions (Fig. 2a). Figure 2a presents the eMSOT grid applied on the images processed (left column), the sO 2 maps obtained with linear unmixing (middle column) and with eMSOT (right column). The spectral fitting of linear unmixing (left) and eMSOT (right) corresponding to a pixel in the area of the capillary tube (yellow arrows in a) are presented in Fig. 2b along with the estimated sO 2 values. In the controlled in vivo experiments, the mean linear unmixing error ranged from 16 to 35% while eMSOT offered a mean sO 2 error ranging from 1 to 4%, indicating an order of magnitude improved accuracy (Fig. 2c).
Imaging blood oxygenation gradients in muscle and tumour. Blood oxygenation and oxygen exchange in the microcirculation have been traditionally studied through invasive, single-point polarography or microscopy measurements in vessels and capillaries of the skeletal muscle 17 . Research for macroscopic methods that could non-invasively resolve muscle oxygenation was broadly pursued in the past two decades by considering nearinfrared spectroscopy and diffuse optical tomography, which, however can only report bulk tissue sO 2 values 18,19 . In a next set of experiments we, therefore, studied whether eMSOT could non-invasively quantify the oxygenation gradient in the skeletal muscle, and we compared this performance to conventional spectral optoacoustic methods. eMSOT was applied in the area of the hindlimb muscle of mice undergoing an O 2 challenge as described in Supplementary Note 5 (n ¼ 6 animal experiments); three of the mice were then killed with an overdose of CO 2 , the latter binding to haemoglobin and deoxygenating blood.
eMSOT applied in the hindlimb muscle area (grid shown in Fig. 3a) resolved oxygenation gradients as a function of breathing conditions in vivo (Fig. 3 b,c) and post mortem after CO 2 breathing (Fig. 3d)  herein as a control experiment and was also analysed with linear unmixing for comparison (Fig. 3e). In the post-mortem case, linear unmixing overestimated the sO 2 as a function of tissue depth (Fig. 3e) and yielded large errors in matching the tissue spectra ( Fig. 3f-upper row). Conversely, eMSOT offered sO 2 measurements in agreement with the expected physiological states ( Fig. 3b-d) and consistently low-fitting residuals ( Fig. 3f-lower row, Supplementary Fig. 6). Figure 3d-f demonstrates the prominent effects of spectral corruption with depth that impair the accuracy of conventional spectral optoacoustic methods, but are tackled by eMSOT. The estimated blood sO 2 values corresponding to a deep tissue area (yellow rectangle in Fig. 3b) are tabulated in Fig. 3g for eMSOT and linear unmixing, and depict that the latter demonstrated unrealistically small sO 2 changes between the normoxic in vivo and anoxic post-mortem (after CO 2 breathing) states. In addition to physiological tissue features, MSOT also reveals tissue morphology. MSOT images at a single wavelength (900 nm) captured prominent vascular structures likely corresponding to femoral vessels or their branches (Fig. 3h) with implicitly co-registered eMSOT blood oxygenation images (Fig. 3i). This hybrid mode enables the study of physiology at specific tissue areas. We selected to study blood oxygenation measurements at a region of interest around large vessels (ROI-1; Fig. 3h) and a region of interest within the muscle presenting no prominent vascular structures (ROI-2; Fig. 3h) for the 100% O 2 , 20% O 2 and CO 2 breathing conditions. Average tissue sO 2 was typically measured at 60-70% saturation under medical air breathing and at 70-80% saturation under 100% O 2 breathing near large vessels (Fig. 3j). Average tissue blood oxygenation away from large vessels (ROI-2) was estimated consistently lower, at 35-50% saturation under normal breathing conditions, and 45-60% saturation under 100% O 2 breathing (Fig. 3k).
The low blood saturation values in tissue (35-50%) cannot be explained by considering arterial and venous blood saturation. However, previous studies based on direct microscopy measurements in vessels and capillaries through polarography, haemoglobin spectrophotometry and phosphorescence quenching microscopy have revealed similar oxygenation gradient in the skeletal muscle 17 with haemoglobin saturation in the femoral artery found to range between 87 and 99% sO 2 17,20 , while rapidly dropping down to 50-60% sO 2 in smaller arterioles 20,21 . The average oxygen saturation in venules and veins has been found to range between 45% and 60% sO 2 under normal breathing conditions, reaching up to 70% at 100% O 2 breathing 21,22 . Average capillary blood oxygenation has been estimated at 40% sO 2 with a large standard deviation 22 , often reported lower, at an average, than venular oxygenation 17 . Therefore, the eMSOT values measured at ROI-1 possibly relate to a weighted average of arterial/arteriolar and venous/venular sO 2 in the skeletal muscle, while the values measured at ROI-2, which anatomically presents no prominent vasculature, relate more to capillary sO 2 measurements.
The improved accuracy observed in eMSOT over previous approaches and general agreement with invasive tissue measurements prompted the further study of perfusion hypoxia emerging from the incomplete delivery of oxygenated haemoglobin in tissue areas. We hypothesized that measurements of blood saturation could be employed as a measure of tissue hypoxia, assuming natural haemoglobin presence in hypoxic areas. To examine this hypothesis, we applied eMSOT to measure blood oxygenation in 4T1 solid tumours orthotopically implanted in the mammary pad of eight mice (Methods, Supplementary Note 6). MSOT revealed the tumour anatomy and eMSOT exposed tumour heterogeneity, which was found consistent to anatomical features identified through cryoslice colour photography and haematoxil and eosin    Supplementary Fig. 7c-g). Furthermore, imaging tumours at different time-points revealed the progression of hypoxia during tumour growth (Fig. 4a,b). The spread of hypoxia, that is, the percentage of the hypoxic area (area presenting sO 2 values below a threshold which varied from 50 to 25% sO 2 ) over the total tumour area also increased during tumour progression (Fig. 4c). Following the in vivo measurements we harvested the tumour tissue and related the non-invasive eMSOT findings to the histological assessment of tumour hypoxia (see Supplementary Note 6 and Supplementary Fig. 7). Tumour tissue was stained by Hoechst 33342 (ref.

23) (indicating perfusion) and
Pimonidazole 24 (indicating cell hypoxia). The results indicated close correspondence between the hypoxic areas detected by eMSOT using haemoglobin as a hypoxia sensor (Fig. 4b) and the histology slices (Fig. 4d). We found that eMSOT could not only quantitatively distinguish between high-and low-hypoxia levels in the tumours, but the spatial sO 2 maps further presented congruence with the spatial appearance of hypoxia spread and reduced perfusion seen in the histology slices ( Fig. 4e-g). A quantitative congruence analysis is shown in Supplementary  Figure 7. Finally, clear differences were also observed between the hypoxic tumour and healthy tissue response to an O 2 breathing challenge ( Fig. 4h; Supplementary Fig. 8), with areas in the core of the tumour presenting a limited response to such external stimuli, likely due to the presence of non-functional vasculature.

Discussion
Spectral corruption has so far limited the potential of optical and optoacoustic methods to offer accurate, quantitative assessment of sO 2 deep inside tissues. Conventional computational methods in optical imaging propose to invert a light transport operator to recover tissue optical properties (absorption and scattering) 13 ; then use these properties for calculating tissue physiological parameters. However, the very-high numerical complexity and ill-posed nature of the inversion problem has not allowed so far accurate, high-resolution sO 2 imaging. We hereby followed an alternative approach that describes the spectral features of light fluence as a combination of spectral base functions. Using this principle, we formulated the sO 2 quantification problem as a nonlinear spectral unmixing problem that does not require knowledge of tissue optical properties or the inversion of a light transport operator. Effectively, eMSOT converts sO 2 imaging from a problem that is spatially dependent on light propagation and optical properties, as common in traditional optical methods, to a problem solved in the spectral domain. Therefore, sO 2 can be directly quantified without estimating tissue optical properties. eMSOT requires theoretically at least five excitation wavelengths for resolving spectral domain parameters and the relative oxygenated and deoxygenated haemoglobin concentrations. We hereby utilized 21 wavelengths for ensuring high accuracy. The recent evolution of video-rate MSOT imaging systems, based on fast tuning optical parametric oscillator lasers 25 , allows the practical implementation of the method. Modern MSOT systems offer five wavelength scans at 20 Hz or better. Therefore eMSOT is a technology that optimally interfaces to a new generation of fast and handheld spectral optoacoustic systems 26 .
The method developed demonstrated quantitative, noninvasive blood oxygenation images in phantoms and tissues in vivo (muscle and tumour) in high resolution, showing good correlation with the expected physiological state or the histologically observed spatial distribution of perfusion and hypoxia. eMSOT measures blood oxygenation. We hypothesized that a  correlation exists to tissue oxygenation and hypoxia measurements by assuming a wide presence of haemoglobin in tissues. We demonstrated congruence (Supplementary Note 6) between traditional invasive histological assays resolving tissue hypoxia and eMSOT analysis. Importantly, not only average values are resolved, but there is a close spatial correspondence between hypoxia patterns resolved by eMSOT non-invasively and histological analysis (Fig. 4, Supplementary Fig. 7). High-resolution non-invasive imaging of blood oxygenation across entire tissues and tumours offers novel abilities in studying physiological and pathological conditions. This goal has been pursued for decades with near-infrared methods, but the strong effects of photon scattering and photon diffusion on the signals detected limited imaging resolution and often impeded accurate quantification 27 . Optoacoustic imaging improves the resolution achieved, over diffuse optical imaging methods but its sO 2 estimation accuracy has been limited so far by depth-dependent fluence attenuation and spectral corruption effects. We showed that conventional spectral optoacoustic methods employing linear unmixing can significantly misestimate blood oxygen saturation values in several cases, including simulations and controlled animal measurements. eMSOT was tested on a vast data set consisting of 42,000 tissue simulations and was consistently found to provide from a comparable to substantially better sO 2 estimation accuracy over linear unmixing. (Supplementary Note 3). The large number of simulations was necessary to validate eMSOT, which presents a non-convex optimization problem. eMSOT was further tested on tissue mimicking blood phantoms (Supplementary Note 4) and controlled in vivo experiments (Fig. 2, Supplementary Note 5). In all cases tested, eMSOT offered from comparable to significantly more accurate performance over conventional spectral optoacoustic methods.
A particular challenge in this study was the confirmation of the eMSOT values obtained in vivo. Polarography measurements are invasive, disrupt the local microenvironment and do not allow to recover spatial information. Nuclear methods using tracers are not well suited for longitudinal studies and utilize tracers that need to distribute in hypoxia areas, that is, areas with problematic supply. Therefore, the results may not directly compare to eMSOT, even though such study is planned as a next step. Blood-oxygen-level dependent MRI only resolves the effects of deoxygenated haemoglobin but cannot observe oxygenated haemoglobin. For this reason, we selected to utilize traditional histology methods, using cryoslicing, which allows to maintain spatial orientation so that eMSOT and histological results could be compared not only in terms of quantity but also in regard to the spatial appearance.
eMSOT proposes a solution to a fundamental challenge in optical and optoacoustic imaging. In the absence of established and reliable methods that can image blood oxygenation, it may be that eMSOT becomes the gold standard method in blood sO 2 studies. Its congruence with tissue hypoxia may also allow a broad application in tissue and cancer hypoxia studies. Nevertheless eMSOT performs optimally when applied on well-reconstructed parts of optoacoustic images (Supplementary Note 5). For this reason, it was selectively applied herein to the part of the image that is within the optimal sensitivity field of the detector employed. An eMSOT advantage is that it is insensitive to scaling factors such as the Grüneisen coefficient or the spatial sensitivity field of the imaging system (Methods). However, due to its scale invariance eMSOT only allows for quantifying blood sO 2 and not absolute blood volume, a goal that will be interrogated in future studies. Next steps further include the eMSOT validation with a larger pool of tissue physiology interrogations spanning from cancer, cardiovascular and diabetes research, relation of physiological phenotypes to metabolic and '-omic' outputs and in clinical application.

Methods
Animal preparation and handling. All procedures involving animal experiments were approved by the Government of Upper Bavaria. For the preparation of orthotopic 4T1 tumour models, 8-week-old, adult female, athymic, Nude-Foxn1 mice (Harlan, Germany) were orthotopically inoculated in the mammary pad with cell suspensions (0.5 million 4T1 cells (CRL-2539)). Animals (n ¼ 8) were imaged in vivo using MSOT after the tumours reached a suitable size. All imaging procedures were performed under anaesthesia using 1.8% isoflurane. In the O 2 challenge experiment, the mouse was initially breathing 100% O 2 and in the following medical air (20% O 2 ). During the O 2 Challenge, the mice were stabilized for a period of 10 min under each breathing condition before MSOT acquisition. For controlled mouse measurements (n ¼ 4), MSOT acquisition was performed on mice under gas anaesthesia and breathing 100% O 2 or 20% O 2 by rectally inserting a capillary tube containing pig blood at 100 or 0% sO 2 oxygenation levels. Mice were killed during MSOT imaging with an overdose of CO 2 or after MSOT acquisition by a ketamine/xylazine overdose. In the following the mice were stored at À 80°C for further analysis.
4T1 cell line was acquired from American type culture collection (ATCC-CRL-2539, #5068892). The cells were authenticated by the Americal type culture collection (ATCC) by several analysis tests: post-freeze viability, morphology, mycoplasma contamination, post-freeze cell growth, interspecies determination; bacteria and fungal contamination. Additional mycoplasma contamination tests were also performed. For the animal studies no randomization, blinding or statistical methods were performed.
MSOT imaging. Optoacoustic imaging was performed using a real-time whole-body mouse imaging scanner, MSOT In Vision 256-TF (iThera-Medical GmbH, Munich, Germany). The system utilizes a cylindrically focused 256-element transducer array at 5 MHz central frequency covering an angle of 270 degrees. The system acquires crosssectional (transverse) images through the animal. The animals are placed onto a thinclear polyethylene membrane. The membrane separates the animals from a water bath, which is maintained at 34°C and is used for acoustic coupling and maintaining animal temperature while imaging. Image acquisition speed is at 10 Hz 28 . Imaging was performed at 21 wavelengths from 700 to 900 nm with a step size of 10 nm, and at 20 consecutive slices with a step size of 0.5 mm. Image reconstruction was performed using a model-based inversion algorithm 29,30 with a non-negativity constraint imposed during inversion and with Tikhonov regularization.
eMSOT method and sO 2 maps. All optoacoustic images P(r,l) obtained over excitation wavelength l were calibrated to correct for the intensity of laser power per pulse, and for the absorption of water surrounding the tissue. With HbO 2 and Hb being the main tissue absorbers in the near-infrared, multispectral optoacoustic images can be related to the concentrations of oxy-and deoxy-haemoglobin through equation (2).
In equation (2), F(r,l) is the space and wavelength-dependent optical fluence, C(r) is a spatially varying scaling factor corresponding to the effects of the system's spatial sensitivity field and the Grüneisen coefficient, e HbO2 ðlÞ and e Hb (l) are the wavelength-dependent molar extinction coefficients of oxygenated and deoxygenated haemoglobin, while c HbO2 r ð Þ and c Hb (r) the associated concentrations at a position r. U(r) is a vector corresponding to the light fluence spectrum at position r, and ||U(r)|| 2 is its norm across all excitation wavelengths at a position r. We define F 0 (r,l) ¼ F(r,l)/||U(r)|| 2, which corresponds to the normalized wavelength dependence of light fluence at a specific position (that is, normalized spectrum of light fluence).
The space-only dependent factors C(r) and ||U(r)|| 2 do not affect the estimation of blood sO 2, which is calculated as a ratio once the relative concentrations of HbO 2 and Hb are known. We define c 0 HbO2 r ð Þ ¼ C 0 (r)c HbO2 r ð Þ and c 0 Hb (r) ¼ C 0 (r)c Hb (r), respectively, where C 0 (r) ¼ C(r)||U(r)|| 2 is a common, space-only dependent scaling factor. Using this notation, equation (2) can be transformed into equation (1). Given c 0 HbO2 r ð Þ and c 0 Hb (r), sO 2 can be computed as in: For the accurate quantitative extraction of the relative concentrations c 0 HbO2 r ð Þ and c 0 Hb (r), accounting for, or estimating the wavelength dependence of the light fluence F 0 (r,l) is further required.
The Eigenspectra model. eMSOT is based on the observation that the spectral patterns of light fluence present in tissue can be modelled as an affine function of only a few base spectra, independently of tissue depth and the specific distribution of optical properties of the tissue imaged. This hypothesis stems from the notion that the spectrum of light fluence is the result of the cumulative light absorption by haemoglobin; thus the spectrum of light fluence will always be related to the spectra of haemoglobin in a complex nonlinear manner. This complex relation can be linearized using a data-driven approach, that is, through the application of PCA on a selected set of light fluence spectra.
The wavelength dependence of the light fluence was herein modelled as a superposition of a mean fluence spectrum F M (l) and a linear combination of a number of light fluence Eigenspectra F i (l). This model was derived by applying PCA on a training data set comprised of a set of light-fluence spectral patterns. Briefly, a training data set was formed through the creation of multispectral light fluence simulations using the 1D analytical solution of the diffusion equation for infinite media. A set of light-fluence spectral patterns F z,ox (l) were computed for high physiological tissue optical properties (m a ¼ 0.3 cm À 1 , m s 0 ¼ 10 cm À 1 ), tissue depths ranging from z ¼ 0 to z ¼ 1 cm with a step size of 0.143 mm (70 discrete depths in total) and for 21 different uniform background blood sO 2 levels (ox E {0%, 5%, 10, y, 100%}). The so computed set of light fluence spectra F z,ox (l) was normalized (F' z,ox (l) ¼ F z,ox (l)/||U z,ox || 2 ) and used in the following as training data in the context of PCA in order to create an affine, 3-dimensional model consisting of a mean fluence spectrum F M (l) and three Eigenspectra F i (l). PCA was used for offering a minimum square error property in capturing the spectral variability of the simulated light fluence spectra, in a linear manner. Three components were selected for providing a relatively simple model while also offering a small model error with respect to the training data set (Fig. 1e). The wavelength dependence of the light fluence F 0 (r,l) at any arbitrary tissue position r can thus be modelled as a superposition of the mean fluence spectrum and three fluence Eigenspectra multiplied by appropriate scalar parameters m 1 (r), m 2 (r) and m 3 (r), (hereby referred to as Eigenfluence parameters) as per equation (4): The so created 3-dimensional affine forward model of the wavelength dependence of light fluence was tested with regard to light-fluence spectral patterns produced in completely heterogeneous media with varying and randomly distributed optical properties and oxygenation values and demonstrated high accuracy ( Supplementary Fig. 1). The forward model was further tested through in vivo and ex vivo light fluence measurements, obtained from controlled experiments ( Supplementary Fig. 2).
Through simulations, it was observed that the values of the m 2 Eigenfluence parameter relate primarily to tissue depth and the average tissue optical properties. This trend was observed both in the case of tissue simulations with uniform optical properties (Fig. 1g), as well as in complex and randomly created tissue simulations described in Supplementary Note 1, 3. Conversely, the values of the Eigenfluence parameters m 1 and m 3 relate both to tissue depth, as well as to tissue background oxygenation. Specifically both m 1 and m 3 present a trend of increasing absolute values with depth and a sign that relates to background blood sO 2 . These observations were confirmed with in vivo and ex vivo light fluence measurement experiments (Supplementary Note 1).
Fluence correction and sO 2 quantification. The minimization of cost function f grid (equation (6)) under the constraints of equation (7) yields an estimate of m i (r p,l ) for each Eigenfluence parameter i and each grid point r p,l . The Eigenfluence parameters in the convex hull of the grid are in the following estimated by means of cubic interpolation. We note that due to the nature of diffuse light propagation, the Eigenfluence parameters are expected to vary smoothly in tissue and thus their interpolation is not expected to introduce large errors in the result (see Supplementary Note 3 and Supplementary Table 2). The wavelength dependence of light fluence is computed for each pixel within the convex hull of the grid as in F 0 (r,l) ¼ F M (l) þ m 1 (r)F 1 (l) þ m 2 (r)F 2 (l) þ m 3 (r)F 3 (l), where F i (l) is the ith fluence Eigenspectrum. Finally, a spectrally corrected eMSOT image is obtained after dividing the original image P(r,l) with the normalized wavelength-dependent light fluence F 0 (r,l) at each position r and wavelength l, i.e., P eMSOT (r,l) ¼ P(r,l)/ F 0 (r,l). The relative concentrations of HbO 2 and Hb (c 0eMSOT HbO2 r ð Þ, c 0 Hb eMSOT (r)) are computed for each pixel of P eMSOT (r,l) image independently through non-negative constrained least squares fitting with the spectra of oxygenated and deoxygenated haemoglobin. Thus, the eMSOT blood sO 2 maps retain the original resolution of the MSOT imaging system.
We note that both the Eigenspectra model and the inversion scheme were hereby optimized for the application of small-animal imaging. The Eigenspectra model was trained for a maximum depth of 1 cm and the inversion scheme was designed with respect to the same tissue depth and optical properties within the physiological range (Supplementary Note 2, 3).
Linear unmixing. Under the simplifying assumption that the light fluence attains a flat spectrum irrespective of the tissue position F(r,l) ¼ F(r) and by assuming haemoglobin as the major absorber in tissue, optoacoustic spectra can be modelled as a linear combination of the spectra of oxy-and deoxy-haemoglobin. The term linear unmixing refers hereby to the computation of the relative concentrations of HbO 2 and Hb (c lu HbO2 r ð Þ, c Hb lu (r)) and subsequently blood sO 2 , through non-negative constrained least squares fitting of the original image P(r,l) with the spectra of Hb and HbO 2 .
Blood phantom preparation. For validating the accuracy of eMSOT in quantifying blood oxygenation in deep tissue, we prepared tissue mimicking phantoms, containing blood at known oxygenations levels. Specifically, for simulating tissue background, 2-cm-diameter cylindrical solid phantoms were created by using 1.5% Agarose Type I, Sigma-Aldrich (solidifying in o37°), 2% intralipid and 3-5% freshly extracted pig blood diluted in NaCl. Different blood oxygenation levels were achieved by diluting oxygen in whole blood (oxygenation process) or by mixing the blood with different amounts of sodium dithionite (Na 2 O 4 S 2 ) (deoxygenation