1.

Introduction

A robust and accurate quantitative measurement of blood oxygen saturation (${\mathrm{sO}}_2$) with optoacoustic (OA) imaging, also called photoacoustic imaging, is one of the most sought after goals of quantitative OA imaging research due to its wide range of immediate applications. Usually, quantitative OA imaging research aims to achieve an absolute quantification of optical properties, such as the absorption coefficient ${\mu }_{\mathrm{a}}$, from measured OA signals $S(\mathbf{d},t)$ recorded at times $t$ at detector position $\mathbf{d}$.^1,2 In brief, such a quantification of ${\mu }_{\mathrm{a}}$ encompasses a solution of two ill-posed inverse problems. (1) The acoustic inverse problem from $S(\mathbf{d},t)$ to an initial pressure spatial distribution $p_0(\mathbf{x})$. And (2), the optical inverse problem from $H({\mathbf{x}}_0)=p_0({\mathbf{x}}_0)/\mathrm{\Gamma }({\mathbf{x}}_0)=\varphi ({\mathbf{x}}_0,{\mu }_{\mathrm{a}}(\mathbf{x}),{\mu }_{\mathrm{s}}^{\prime }(\mathbf{x}))·{\mu }_{\mathrm{a}}({\mathbf{x}}_0)$ to ${\mu }_{\mathrm{a}}({\mathbf{x}}_0)$, at a location ${\mathbf{x}}_0$, with the Grüneisen parameter $\mathrm{\Gamma }$ and the reduced scattering coefficient ${\mu }_{\mathrm{s}}^{\prime }$. The fluence $\varphi$ is dependent on unknowns such as the absorption and scattering in the tissue surrounding ${\mathbf{x}}_0$. Quantitative OA imaging methods either depend on model-based inversion^2–7 or data-driven approaches.^8–13 These approaches perform well in silico but often struggle with the translation to real measurements in phantoms or in vivo.

In OA imaging, ${\mathrm{sO}}_2$ estimations are derived from multispectral OA measurements by first performing an acoustic reconstruction yielding images of the OA signal

Both these methods can yield promising results on their own but still suffer from a range of constraints, i.e., MI sensing implementations¹⁶ typically assume and use point illuminations, which enables the use of closed-form solutions of the diffusion approximation of light propagation¹⁴ but limits SNR due to the laser safety limit for skin.¹⁷ The resulting long acquisition times make this method difficult to translate to realistic macroscopic applications.¹⁸ Furthermore, MI sensing so far has theoretical limits in highly inhomogeneous scenes due to its reliance on the diffusion approximation. MI sensing implementations usually aim to estimate absolute values of ${\mu }_{\mathrm{a}}$, which goes beyond what is needed for an estimation of ${\mathrm{sO}}_2$. LSD^15,19 and similar spectral approaches³ currently yield accurate in silico estimations and plausible initial results in highly constrained settings, but they have insufficient input to robustly generalize these results over diverse geometries and applications.

We will investigate LSD as a method to analyze MI data. Both MI sensing and LSD are not yet thoroughly validated; partially due to a lack of stable and reliable ${\mathrm{sO}}_2$ phantoms. Even though substantial progress has been made in dynamic blood flow phantoms for OA imaging validation, these blood or red blood cell suspension phantoms require extensive fine tuning and even then yield reference values with limited accuracy.²⁰ At best, a reference measurement of 2% to 4% is achievable with state-of-the-art blood flow phantoms.^21,22 While validating quantitative OA oximetry methods, the validation phantoms are also often restricted to the extreme ${\mathrm{sO}}_2$ values of 0% and 100% because other values cannot be set reliably.²³ This causes an incomplete range and therefore insufficient validation.

Rather than implement such an ${\mathrm{sO}}_2$ flow phantom, we used copper and nickel sulfate solutions in a relative copper sulfate model similar to work by Buchmann et al.²⁴ to mimic absorption spectra of blood with different oxygen saturation. This allowed a reliable sub 1% error in our ground truth and allowed us to rapidly manufacture stable and highly reproducible phantoms with wide variations in optical properties to generate high quality test sets for spectral decoloring methods.

2.

Materials and Methods

We investigated a method combining LSD and MI measurements. To that effect we

2.1.

Multiple Illumination Optoacoustic Imaging

Our MI OA imaging setup is shown in Fig. 1. It uses a fast wavelength-tunable optical parametric oscillator (OPO) laser system (prototype SpitLight, InnoLas Laser GmbH, Krailling, Germany) with 5-ns pulse duration and 100-Hz pulse repetition frequency. The laser pulses were sequentially coupled into four high power fiber bundles (FiberOptic P.+P. AG, Spreitenbach, Switzerland) with NA 0.22 fibers, each bundle with a 2-mm diameter. This was achieved using a galvo mirror system (GVS011/M, Thorlabs Inc., Newton, New Jersey, USA) driven by an arbitrary waveform generator (AWG) (TG5011, Aim-TTi, Cambridgeshire, UK), which was synchronized with the laser system. The fiber bundle output sides were arranged in a line array with 8 mm spacing. The illumination pulses were attenuated to have a maximum energy of 10 mJ per pulse at the fiber output. To comply with ANSI safety limits,^17,25 the beams are widened to 7-mm full-width at half-maximum (FWHM) at the tissue or phantom surface. Illumination and acoustic detection ensue through 18-mm-thick ultrasound gel pad (Parker Laboratories Inc., Fairfield, New Jersey, USA). We measure the 64 center channels of a 128-element linear array transducer (L7-4, Advanced Technology Laboratories Inc., Bothell, Washington, USA) with a center frequency of 5 MHz, a pitch of 0.3 mm, and a fractional bandwidth of 80%. The number of acquisition channels was limited by our 64 channel US data acquisition system (V-1-64, Verasonics, Inc., Kirkland, Washington, USA). For this study, we used the full tuning range of our OPO and acquired OA measurements for 16 equidistant wavelengths from 680 to 980 nm in 20 nm steps, each for four illumination positions. After firing one pulse of one wavelength in each fiber bundle, the wavelength is tuned to the next in sequence. Using this $4\times 16$ sequence, each MI and multispectral stack of OA images takes 640 ms to acquire. We generally recorded the raw data for 30 such stacks for each scan. Live beamforming and visualization with 25 fps was performed using custom MATLAB scripts but this live visualization was solely used for probe positioning and quality control (e.g., avoiding air inclusions under the gel pad).

Fig. 1

MI OA imaging setup. Illumination via fast tunable OPO laser sequentially illuminating fiber bundles using a galvo mirror system driven by an AWG. OA signals were measured with a linear array ultrasound (US) probe and recorded by a 64-channel US data acquisition (DAQ) system. An US gel pad is used to allow for in plane illumination.

2.3.

Phantoms

The phantoms used consisted of arrays of polythene tubing (Smiths Medical International Ltd., Kent, UK) with 0.58-mm inner diameter and 0.96-mm outer diameter. These tubes were filled with a relative copper sulfate model solution (as detailed in Sec. 2.3.1) and arranged as shown in Sec. 2.3.3. The relative copper (rCu) in this model is mimicking blood oxygenation (${\mathrm{sO}}_2$).

Selecting the same small size tubes allowed us to rapidly assemble and modify phantoms with many target structure locations. The small size of the tubes was also chosen because the rCu solution in the tube did not include a scattering agent.

For all the phantom experiments, the background scattering medium was a fat emulsion (SMOFlipid 20%, Fresenius Kabi, Switzerland) diluted to 1.5% fat content. To avoid errors introduced by interbatch variations in the scattering properties of stock fat emulsions, such as intralipid or SMOFlipid, the optical properties of the used stock emulsion were assessed with a time-correlated single photon counting (TCSPC) technique as detailed in Sec. 2.3.2.

2.3.1.

Relative copper sulfate model

The relative copper sulfate model solution was based on a 2.2-molar nickel sulfate (${\mathrm{NiSO}}_4$) water solution, produced using nickel(II) sulfate hexahydrate ($>98\%$, Sigma-Aldrich), and on a 0.25 molar copper sulfate (${\mathrm{CuSO}}_4$) water solution, produced using copper(II) sulfate pentahydrate ($>98\%$, Sigma-Aldrich).²⁷ As shown in Fig. 2, these chromophores are mimicking the NIR absorption spectra of oxy- and deoxyhemoglobin in average whole blood with a hemoglobin concentration $c_{\mathrm{wb}}(\mathrm{HbT})=150{\mathrm{gl}}^{-1}$.²⁸ Copper and nickel sulfate were also chosen for their temporal stability and resistance to bleaching.

Fig. 2

Absorption coefficient ${\mu }_{\mathrm{a}}$ spectra. (a) Oxy- and deoxyhemoglobin at whole blood concentrations $c_{\mathrm{wb}}(\mathrm{HbT})=150{\mathrm{gl}}^{-1}$. (b) Copper and nickel sulfate in aqueous solution in whole blood equivalent solutions using a relative copper sulfate (rCu) model. The reference measurements of the five rCu mixtures used in our phantoms are plotted as “+.” Sulfate spectra were measured with a spectrophotometer.

The spectra of the sulfate solutions absorption coefficients ${\mu }_{\mathrm{a}}$ in whole blood mimicking concentrations are defined as $c_{\mathrm{wb}}$ (${\mathrm{NiSO}}_4$) :=2.2 M and $c_{\mathrm{wb}}$ (${\mathrm{CuSO}}_4$) :=0.25 M. These solutions were measured using a 2-mm quartz cuvette (QS Hellma, Müllheim, Germany) in a UV–VIS–NIR spectrophotometer (Perkin Elmer Lambda 750, Waltham, Massachusetts, USA), in the range of 680 to 980 nm. The scattering in this wavelength range is negligible.²⁷ The initial reference measurements were done in 2 nm steps, with a 10-s integration time and using a photomultiplier tube sensor. The absorption spectroscopy measurements were repeated on the solutions after 70 days to verify their stability over time. Whenever new batches of the sulfate solutions were produced, their absorption spectra were checked against the spectra of the first batch. The solutions were corrected when they deviated from the reference spectra by more than 1%.

The relative copper (rCu) in this model aims to mimic blood oxygenation (${\mathrm{sO}}_2$) and is therefore similarly defined as

Eq. (2)

$$\mathrm{rCu}=\frac{c_{\mathrm{r}}({\mathrm{CuSO}}_4)}{c_{\mathrm{r}}({\mathrm{CuSO}}_4)+c_{\mathrm{r}}({\mathrm{NiSO}}_4)},$$

with the respective concentrations of the sulfate solutions relative to their blood mimicking base solutions

Eq. (3)

$$c_{\mathrm{r}}({\mathrm{CuSO}}_4)=\frac{c({\mathrm{CuSO}}_4)}{c_{\mathrm{wb}}({\mathrm{CuSO}}_4)}\text{and}{c}_{\mathrm{r}}({\mathrm{NiSO}}_4)=\frac{c({\mathrm{NiSO}}_4)}{c_{\mathrm{wb}}({\mathrm{NiSO}}_4)}.$$

For comparison, the definition of blood oxygen saturation is

Eq. (4)

$${\mathrm{sO}}_2=\frac{c({\mathrm{HbO}}_2)}{c({\mathrm{HbO}}_2)+c(\mathrm{Hb})}.$$

While of course not following hemoglobin spectra exactly, this sulfate model is a good qualitative fit to hemoglobin and is much easier to accurately control and reproduce than the saturation of oxygen in hemoglobin. It is highly stable over time; i.e., over 70 days only changes $<1\%$ in absorption were observed. Mimicking the blood volume fraction (BVF) in tissue, we define a sulfate volume fraction (SVF) in our model as $\mathrm{SVF}=c_{\mathrm{r}}({\mathrm{CuSO}}_4)+c_{\mathrm{r}}({\mathrm{NiSO}}_4)$. The SVF within the blood vessel mimicking tubing was always 100% mimicking whole blood, whereas the SVF in the background was varied as detailed in Sec. 2.3.3.

2.3.2.

Optical property reference measurements of phantoms

In the background medium, the scattering comparable to tissue (i.e., ${\mu }_{\mathrm{s}}^{\prime }=15{\mathrm{cm}}^{-1}$ at 750 nm) was obtained using a 1.5% fat emulsion (diluted from SMOFlipid 20%, Fresenius Kabi, Switzerland).

To ensure a reproducible and tissue mimicking scattering, the background medium was analyzed with TCSPC spectroscopy. The TCSPC instrument used for the spectral analysis of the emulsions optical properties consisted of a white light supercontinuum laser (SuperK Extreme, NKT Photonics, Birkerød, Denmark) with $\approx 100\mathrm{ps}$ pulse duration (varying with wavelength), running at 39 MHz with $<4\mathrm{mW}$ laser output. This white light was filtered by a tunable filter (SuperK Varia, NKT Photonics, Birkerød, Denmark), which was tuned in a range from 600 to 840 nm in 20 nm steps, with a bandwidth of 10 nm; 840 nm being the maximum of the tunable filter’s range. A single-photon avalanche diode (MDP PDM Series, Micro Photon Devices, Bolzano, Italy) was used to detect single photons. The diode has a prolonged dead time of $\approx 80\mathrm{ns}$ after a photon detection. Because of that, the photon detection rate was kept sufficiently low to make photon detection events during the dead time unlikely. We ensured a detection rate lower than ${10}^5{\mathrm{s}}^{-1}$ ($\ll 1/80\mathrm{ns}$), making a correction for missed photons during the dead time unnecessary. The distributions of times of flight were recorded with single photon counting electronics (SPC-160, Becker & Hickl GmbH, Berlin, Germany). Source and detector fibers were fixed in blunted hypodermic needles for stability. The laser pulse shape, temporal dispersion in the optical fibers, and response of the detector were characterized in the overall instrument response function (IRF), yielding an FWHM of $\approx 140\mathrm{ps}$ overall, varying with wavelength. The source and detection fibers were placed perpendicular to the surface of the sample medium and immersed in the medium by 0.5 mm. To reduce the detection of early arriving photons, a carbon fiber mesh blocker was placed into the direct path, at a distance of 6 mm from the source fiber (dimension: 1 mm depth, 4 mm width, 0.4 mm thickness). We measured the SMOFlipid 1.5% medium in an 8 cm radius, 10 cm deep beaker, with the fibers at the center. This is a sufficiently large volume to be approximated as a semi-infinite medium for the analytic diffusion model. The resulting media were both measured with a source detector separation $\rho =20\mathrm{mm}$, for each wavelength until at least ${10}^7\text{photons}$ were detected. For some wavelengths, the laser needed to be attenuated to keep the photon detection rate below ${10}^5{\mathrm{s}}^{-1}$. This acquisition protocol ensured a high signal-to-noise ratio (SNR) and allowed us to fit our diffusion model only to late arriving photons where the diffusion approximation is more accurate. For the phantom experiments, two bottles of a new batch of SMOFlipid were used—both batches and bottles were measured independently prior to experiments to avoid hidden variations in the background medium.

An analytic diffusion model²⁹ with an extrapolated boundary condition for a semi-infinite medium^30,31 was convolved with the corresponding IRF for each wavelength $\lambda$. The results were then fitted to the measured histograms of the single photon arrival times, yielding a series of tuples (${\mu }_{\mathrm{s}}^{\prime \mathrm{SPC}}(\lambda )$, ${\mu }_{\mathrm{a}}^{\mathrm{SPC}}(\lambda )$). Our tunable filter was limited in range to a maximum wavelength 840 nm but we needed credible ${\mu }_{\mathrm{s}}^{\prime }$ values up to 980 nm for the optical forward simulations. Therefore, a generic tissue model [Eq. (5)] from the mcxyz framework³² was used to expand and define the scattering properties within the optical forward simulation:

Eq. (5)

$${\mu }_{\mathrm{s}}^{\prime }(\lambda )={\mu }_{\mathrm{s}500}^{\prime }·(f_{\mathrm{ray}}·{(\lambda /500\mathrm{nm})}^{-4}+(1-f_{\text{ray}})·{(\lambda /500\mathrm{nm})}^{-b_{\mathrm{mie}}}),$$

With ${\mu }_{\mathrm{s}500}^{\prime }=42.4{\mathrm{cm}}^{-1}$ the initial guess for ${\mu }_{\mathrm{s}}^{\prime }$ at 500 nm, $f_{\text{ray}}=0.62$ the initial guess for fraction of Rayleigh scattering at 500 nm, and $b_{\mathrm{mie}}=1.0$ the initial guess for the scatter power for Mie scattering. This was fitted to the TCSPC data with a least squares fit—the entire data processing pipeline with all parameters is part of the open source code supplement (see the Code, Data, and Materials Availability section). The resulting fits are shown in Fig. 3.

Fig. 3

Optical properties of the uncolored phantom background medium. Single data points are diffusion model results from measurements with the TCSPC spectroscopy instrument of a 1.5% fat emulsion (diluted from SMOFlipid 20%). Each data point corresponds to a fit on a TCSPC histogram of at least ${10}^7\text{photons}$ collected over at least 100 s. Validation phantoms cf. Fig. 4(a) were constructed with the “second batch, first bottle.” A generic tissue scattering model [Eq. (5)] fitted to this first bottle measurement was used to set the background scattering properties (a) for the Monte Carlo simulations. The “second bottle” was used for the background media in the test phantoms. The absorption results (b) are shown together with a literature spectrum of water absorption³³ (dashed line).

2.3.3.

Phantom data sets

Three sets of phantoms (A,B,C) were produced, with different layout as shown in Fig. 4. All phantoms use polythene tubing filled with the relative copper sulfate model solution as target structures. The phantom backgrounds consist of a 1.5% fat emulsion with added sulfates.

Fig. 4

Cross sections of the phantom data sets, with denoted parameters: relative copper sulfate model in the tubes ${\mathrm{rCu}}_{\text{tube}}$, in the background medium ${\mathrm{rCu}}_{\mathrm{bg}}$ and sulfate model volume fraction in the background medium ${\mathrm{SVF}}_{\mathrm{bg}}$. (a) The validation phantoms with five single tubes. (b) The main test phantoms. (c) The test phantoms in longitudinal scan direction and thereby somewhat OOD of the training data. The shown two-dimensional cross sections correspond to the imaging plane. In sets A and B, the tubes run perpendicular to the imaging plane. Phantom test C has the same geometry as set B, with the imaging plane rotated by 90 deg to yield longitudinal scans instead of transversal scans of the tubes.

Phantom layout A was measured as a validation data set for hyperparameter tuning of the machine learning models and validation of image reconstruction as well as parameter tuning in the Monte Carlo simulations. Layouts B and C were exclusively measured as test data sets. Phantom test set B is expected to be within the distribution of the simulation parameters (cf. Fig. 5). Phantom test set C, however, consists only of longitudinal scans w.r.t. the tube orientation. Because the orientation of the illumination positions changes with the imaging plane, set C was illuminated along the tubing. The measurements in set C are therefore expected to be out-of-distribution (OOD) with respect to the Monte Carlo simulated training sets. As detailed in the next section, the simulations were exclusively performed for transversal orientation of the tubing.

Fig. 5

The in silico training data set consists of 4000 volumes, simulated with Monte Carlo simulations, modeling the geometry of the real MI setup. An additional 1000 volume test set is kept separate. Each volume has two sets of tubes each with a random number of tubes, uniformly distributed as specified. Tube and background relative copper (rCu) as well as background SVF are also drawn from uniform random distributions $U$. The ${\mathrm{sO}}_2$ training and test sets are simulated identically, substituting rCu absorption spectra for hemoglobin spectra, cf. Fig. 2.

The phantom data sets contain 164 multispectral MI OA scans from 115 scan configurations as follows:

• A. 30 scan configurations as laid out in Fig. 4(a): six phantom configurations, one with only a 1.5% SMOFlipid background solution and five with an added 1% SVF background with relative copper ${\mathrm{rCu}}_{\mathrm{bg}}$ set to {0, 25, 50, 75, 100}%. On each of these six configurations, five MI multispectral scans were performed centering the transducer on each of the tubes with ${\mathrm{rCu}}_{\text{tube}}=\{\mathrm{0,25,50,75,100}\}\%$.

• B. 55 scan configurations as laid out in Fig. 4(b): eleven phantom configurations, one with only a 1.5% SMOFlipid background solution, five with an added 1% SVF background with ${\mathrm{rCu}}_{\mathrm{bg}}=\{\mathrm{0,25,50,75,100}\}\%$, and five with a 0.5% SVF. On these 11 phantom configurations, MI multispectral scans were performed centering on each of the five four-tube-arrays with ${\mathrm{rCu}}_{\text{tube}}=\{\mathrm{0,25,50,75,100}\}\%$. For each four-tube-array, two regions of interest (ROI) (one containing the two lower and one the two upper tubes) were analyzed separately. The imaging plane was positioned for transversal scans of the tubes.

• C. 30 scan configurations as laid out in Fig. 4(c): three phantom configurations, one with only a 1.5% SMOFlipid background solution, two with an added 1% SVF background with ${\mathrm{rCu}}_{\mathrm{bg}}=\{\mathrm{0,100}\}\%$. On these three phantom configurations, MI multispectral scans were performed with each of the five shallowest tubes, and each of the five deepest tubes of the four-tube-arrays in the imaging plane, with ${\mathrm{rCu}}_{\text{tube}}=\{\mathrm{0,25,50,75,100}\}\%$. The imaging plane was positioned for longitudinal scans of the tubes.

All scan configurations were scanned for 19.2 s yielding 30 MI and multispectral sequences. Due to the limited field of view of our US system (parallel read-out of 64 channels on a 19.2-mm linear array), we repositioned the probe between acquisitions—i.e., measuring five scan positions for phantom geometries A and B. The center of the linear transducer was always placed above the center of the targeted tubes. Scans with technical difficulties such as frame drops or wrong positioning were discarded in postprocessing. This affected one of the 115 scan configurations: the ${\mathrm{rCu}}_{\text{tube}}=100\%$, ${\mathrm{rCu}}_{\mathrm{bg}}=50\%$, $\mathrm{SVF}=0.5\%$ was discarded for erroneous positioning. All scans of phantom geometry C were performed twice. The $\mathrm{SVF}=0$ scans on phantom geometry B were performed five times on different days as a baseline measurement. The total phantom data set consists of 164 scans.

It is important to note that both copper and nickel sulfate act as a demulsifier when mixed with the diluted SMOFlipid background or any other fat in water emulsion. Phases will form and the bulk optical properties will change significantly within tens of seconds. To avoid the forming of phases, the background medium with added sulfates was continuously stirred with a magnetic stirrer during all the measurements.

2.5.

Machine Learning Algorithms

The estimation of an ${\mathrm{sO}}_2$ or rCu value from a measured spectrum is a regression problem. The usual approach to this problem in OA imaging is linear spectral unmixing (LU).^26,35 For one pixel, the OA signal spectrum $\mathbf{S}(\mathbf{\lambda })$ is measured at a set of wavelengths $\mathbf{\lambda }$. This sampled OA signal spectrum $\mathbf{S}$ is then fitted to a linear combination of known absorption spectra. Here, LU is performed numerically using an iterative least squares solver implemented in Python’s scipy.optimize submodule. These LU estimations (${\mathrm{rCu}}_{\mathrm{est}}^{\mathrm{LU}}$) are given throughout the results section as a reference.

We also compare our results to LSD, a type of machine learning algorithm. LSD also aims to estimate ${\mathrm{sO}}_2$ or rCu from the same single illumination OA signal spectra $\mathbf{S}$ measured at wavelengths $\mathbf{\lambda }$. Similar to prior implementations, our modified LSD models are machine learning algorithms that are trained on large amounts of simulated absorbed energy spectra labeled with ground truth rCu. Before training, each absorbed energy spectrum is normalized with its $L1$ norm to $\widehat{\mathbf{H}}(\mathbf{\lambda })$. This normalization makes them equivalent to a normalized OA signal spectrum $\widehat{\mathbf{S}}(\mathbf{\lambda })$. This is because we can assume that for a signal spectrum $\mathbf{S}$ at a position ${\mathbf{x}}_0$

Eq. (6)

$$\mathbf{S}({\mathbf{x}}_0,\mathbf{\lambda })=\mathrm{\Gamma }({\mathbf{x}}_{\mathbf{0}})·A({\mathbf{x}}_0)·\mathbf{H}({\mathbf{x}}_0,\mathbf{\lambda })$$

Eq. (7)

$$\Rightarrow \widehat{\mathbf{S}}({\mathbf{x}}_0,\mathbf{\lambda })\approx \widehat{\mathbf{H}}({\mathbf{x}}_0,\mathbf{\lambda }).$$

Assuming a linear acoustic reconstruction such as DAS, $A({\mathbf{x}}_0)$ is a spatially varying but wavelength-independent factor introduced by the imperfect acoustic reconstruction, the instrument response, and the calibration. $\mathrm{\Gamma }({\mathbf{x}}_{\mathbf{0}})$, as a material property is also independent of the illumination wavelength.²⁷ The LSD model, which was trained on the in silico training set tuples $(\widehat{\mathbf{H}},{\mathrm{rCu}}_{\text{tube}})$, is then presented (1) unseen in silico test set spectra $\widehat{\mathbf{H}}$ or (2) spectra $\widehat{\mathbf{S}}$ from an unseen phantom data test set to estimate the corresponding ${\mathrm{rCu}}_{\mathrm{est}}^{\mathrm{LSD}}$.

Note that $A$ actually does depend on the fluence distribution $\mathbf{\varphi }(\mathbf{x},{\mu }_{\mathrm{a}}(\mathbf{x},\mathit{\lambda }),{\mu }_{\mathrm{s}}^{\prime }(\mathbf{x},\mathit{\lambda }))$. A varying optical wavelength may lead to different acoustic spectra of the OA signal corresponding to the same structure, due to different spatial distributions in the absorbed energy. Our assumption is that this effect is small compared with the spectral coloring introduced directly by the fluence term in

Eq. (8)

$$\mathbf{H}({\mathbf{x}}_0,\mathbf{\lambda })=\mathbf{\varphi }({\mathbf{x}}_0,{\mu }_{\mathrm{a}}(\mathbf{x},\mathbf{\lambda }),{\mu }_{\mathrm{s}}^{\prime }(\mathbf{x},\mathbf{\lambda }))·{\mu }_{\mathrm{a}}({\mathbf{x}}_0,\mathbf{\lambda }).$$

For MI-LSD, we have multiple such normalized spectra $\widehat{\mathbf{S}}$ as input variables. For illustration, Fig. 6 shows spectra of the same pixel in an absorber with ${\mathrm{rCu}}_{\text{tube}}=100\%$ with two example illuminations $I_0$, $I_1$ and for two backgrounds with ${\mathrm{rCu}}_{\mathrm{bg}}=\{0\%,100\%\}$ and $\mathrm{SVF}=1\%$. The difference in background absorption causes a different spectral coloring but so does a variation of the illumination position. We hypothesize that training our machine learning algorithms on, i.e., four such spectra will allow us a more accurate and/or more robust estimation compared with LU and LSD.

Fig. 6

Examples for spectral coloring in L1 normalized spectra of an absorber with ${\mathrm{rCu}}_{\text{tube}}=100\%$. Comparison between (a) in silico absorbed energy spectra $\widehat{\mathbf{H}}$ and (b) phantom OA signal spectra $\widehat{\mathbf{S}}$. Spectra for two background media are shown: ${\mathrm{rCu}}_{\mathrm{bg}}=100\%$ (yellow), ${\mathrm{rCu}}_{\mathrm{bg}}=0\%$ (dark); with $\mathrm{SVF}=1\%$. The spectra for two illumination positions $I_0$ (line) and $I_1$ (dotted) are shown as an example for two of the four illuminations. Systematic changes in the spectral coloring can be seen for different background media and illumination positions. These changes are qualitatively similar for $\widehat{\mathbf{H}}$ and $\widehat{\mathbf{S}}$. On the validation phantom examples LU estimations for single spectra ${\mathrm{rCu}}_{\mathrm{est}}^{\mathrm{LU}}$ are listed—spectral coloring can cause large estimation errors relative to the ${\mathrm{rCu}}_{\text{tube}}=100\%$ ground truth.

Two types of machine learning algorithms were employed for spectral decoloring: feed forward neural networks (NN) and random forests (RF). Training of the MI-LSD models includes mirrored illumination positions for each volume as a minor data augmentation. Sorting the training data illumination position spectra stacks by their L1 norm before training was also investigated but did not prove beneficial on the validation data.

2.5.1.

Feedforward neural networks

Feedforward NNs were previously used for LSD implementations.^15,19 We used this state-of-the-art NN architecture as a starting point and further tuned the hyperparameters on the training and validation sets. Doing so we mainly found the dropout layers of previous implementations to be counterproductive—dropout leading to a much lower precision on the validation set. The two final NNs used for both LSD and MI-LSD consisted of four hidden layers with twice the size of the input layer (16 for LSD and 64 for MI-LSD), all with leaky ReLu activation layers (and for comparison and dropout layers). For comparison to the previous implementation, additional results for a dropout in the dropout layers with probability $p=0.2$ are presented in Figs. S35–S50 and S66–S72 in the Supplementary Material. In the main results, no dropout was used ($p=0$). We segmented all vessels in the 4000 volume training set and trained on the segmented 1,052,152 simulated MI signal spectra for 100 epochs. As in the previous implementation, we used a batch size of ${10}^5$ and a learning rate of ${10}^{-2}·{0.9}^{\text{epoch}/2}$. All implementations are documented in the open source appendix. The trained models are also available in the open data appendix. We trained the algorithms on an RTX 2060 Super GPU (Nvidia, Santa Clara) and used the CPU for inference.

2.5.2.

Random forest regression

We also investigated RF regression,³⁶ usually a highly accurate learning algorithm for regression problems with few dimensions.⁸ RFs are also usually less impacted by noise models. In particular, they should not overfit to noise.³⁶ This should prove useful as we did not try to model a realistic wavelength-dependent noise. We used the Python scikit-learn v0.23 implementation of RF regressors using 100 trees with a maximum depth of 30 to limit memory consumption. Further parameters were set to default.

3.

Results

We first show some qualitative comparisons between in silico rCu and phantom data and then present the performance of our trained RF and NN models on our in silico rCu test set and the two phantom test sets.

The hyperparameters of the machine learning models were tuned on the phantom validation set. The rCu machine learning models that performed best on our validation data were used to estimate rCu from the test sets. These models are presented in the results. For further information, all estimations for all models (on the validation set and for every single test measurement) can be found in the figures in the Supplementary Material; representative examples are shown here.

We trained the same RF and NN models, using the same hyperparameters, on an additional ${\mathrm{sO}}_2$ training set.

We compare MI-LSD with LSD and LU. Comparing a method based on a single measurement with a method based on multiple such measurements, the multiple measurement method should generally be more accurate simply due to an increase in SNR. To more fairly compare MI-LSD with the single spectrum methods such as LU and LSD, we estimated LSD and LU results on the reconstructed signals, averaged over the four illuminations. Using this averaged illumination spectrum as input for LU and LSD, we can compare methods for the same delivered energy during the same time—giving no method an inherent SNR advantage. LSD was also trained on in silico data using the same averaged illumination spectra from four simulated illuminations.

In addition to using the validation data set for hyperparameter tuning, we also qualitatively compared a set of our measurements with Monte Carlo simulations of one of the validation phantoms, creating an exact in silico representation of the light propagation in the validation phantom. Figure 7 serves as a qualitative (phantom to in silico) comparison for some of the averaged illumination spectra.

Fig. 7

Qualitative comparison between spectra in a validation phantom (for $\mathrm{SVF}=0$) and its digital twin from Monte Carlo (MC) simulations, showing the effects of various spectral coloring on the mean illumination spectra. Relative copper ${\mathrm{rCu}}_{\text{tube}}$ is varied in the target tube (up-down) and the background medium ${\mathrm{rCu}}_{\mathrm{bg}}$ (left-right). For reference, linear unmixing (LU) rCu estimates are given for each spectrum.

We report the estimation error distributions on the three distinct test sets. Reported are rCu estimation errors $\mathrm{\Delta }{\mathrm{rCu}}_{\mathrm{est}}={\mathrm{rCu}}_{\mathrm{est}}-{\mathrm{rCu}}_{\text{tube}}$ and their absolutes $|\mathrm{\Delta }{\mathrm{rCu}}_{\mathrm{est}}|$, with ${\mathrm{rCu}}_{\text{tube}}$ being the ground truth rCu in the tube. In the in silico test set, all selected models are in close agreement. As shown in Fig. 8, both LSD and MI-LSD estimations of rCu with both RFs and NNs yield median $Q_2$ absolute estimation errors of $<3$ percentage points (pp). The same can be seen in Fig. 9 for the in silico ${\mathrm{sO}}_2$ test set. As expected, estimation with all used models is very fast compared to LU. Inference on CPU for all the 266,105 samples in the in silico test sets took 1.6 s for RF MI-LSD, 1.3 s for RF LSD, 0.2 s for NN MI-LSD, 0.04 s for NN LSD, compared to 642 s for LU.

Fig. 8

Error distributions of the in silico rCu test set cf. Fig. 5. (a) rCu estimation errors $\mathrm{\Delta }{\mathrm{rCu}}_{\mathrm{est}}$ and (b) their absolutes. Blue shows the rCu estimators using MI-LSD, orange the estimators using LSD, and gray is the LU reference. Medians $Q_2$ of the error distributions are shown, together with interquartile ranges (IQR) and 90th percentiles $P_{90}$. The two feedforward NN models and the two RF models all have median absolute errors below 3 pp.

Fig. 9

Error distributions of the in silico ${\mathrm{sO}}_2$ test set cf. Fig. 5. (a) ${\mathrm{sO}}_2$ estimation errors $\mathrm{\Delta }{\mathrm{sO}}_2^{\mathrm{est}}$ and (b) their absolutes. Blue shows the ${\mathrm{sO}}_2$ estimators using MI-LSD, orange the estimators using LSD and gray is the LU reference. Medians $Q_2$ of the error distributions are shown, together with IQR and 90th percentiles $P_{90}$. The two feedforward NN models and the two RF models all have median absolute errors below 3 pp.

From phantom test set B, tube signal was segmented by thresholding. In each reconstructed MI-multispectral OA image stack, two ROIs were chosen: one containing the two upper tubes and one containing the lower two tubes. One such lower tubes ROI is shown in Fig. 10. Each ROI has a fixed size of $3.75\mathrm{mm}\times 3.3\mathrm{mm}$. The 15% highest mean (over all wavelengths and illuminations) OA signal pixels in each ROI were segmented as tube and rCu was estimated from the MI-multispectral OA signals in all pixels of these tube signal areas. The ROIs were thresholded separately to get an equal number of lower tube samples into the test set. A thresholding on the entire image or a larger ROI, using a lower cut-off percentage, would lead to more clutter and noise in the test set and the lower tubes being underrepresented in the test set.

Fig. 10

Example ROI in the phantom test set B with the estimation results for MI-LSD, LSD, and LU. Shown are the lower two tubes of a four tube phantom with $\mathrm{SVF}=0.5\%$ and ${\mathrm{rCu}}_{\mathrm{bg}}=25\%$. To indicate the content of the ROI, the mean OA signal in the ROI is shown left, with the ground truth ${\mathrm{rCu}}_{\text{tube}}$ annotated. The brightness of the OA signal is independently and linearly autoleveled for each ROI. The mean rCu estimate $\overline{\mathrm{rCu}}$ over the ROI is noted for the three estimators.

From phantom test set C, the tube signal was segmented in a similar fashion: from each reconstructed MI-multispectral OA signal image stack an ROI of fixed size ($7.5\mathrm{mm}\times 1.5\mathrm{mm}$) was selected, containing either the upper tube or the lower tube. Two such lower tubes example ROIs are shown in Figs. 11(a) and 11(b) for varying reference ${\mathrm{rCu}}_{\text{tube}}$. Within these ROIs, the 50% highest mean (over all wavelengths and illuminations) OA signal pixels were segmented as tube. rCu was then estimated from the MI-multispectral OA signals in all pixels within these tube signal locations.

Fig. 11

Example ROI 8 mm deep in the phantom test set C with estimation results for MI-LSD, LSD, and LU. Shown are two ranges of five imaged tubes with their ${\mathrm{rCu}}_{\text{tube}}$ annotated above their mean OA signal. The brightness of the OA signal is independently and linearly autoleveled for each ROI. The mean rCu estimate $\overline{\mathrm{rCu}}$ over the ROI is noted for the three estimators. The ROIs for two sets of phantoms are shown. (a) A representative result (${\mathrm{rCu}}_{\mathrm{bg}}=100\%$, 1% SVF background), (b) a result with outlier estimation errors (0% SVF background). LSD has highest estimation errors in deep vessels and in phantoms with no added sulfates in the background medium, i.e., in (b) for ${\mathrm{rCu}}_{\text{tube}}=100\%$.

The estimated rCu image examples from the test sets are shown for the RF models, because with the exception of the in silico test set, NN models performed similarly or worse than RF models. For all estimated rCu images from all models, see the figures in the Supplementary Material. The error distributions for phantom test set B are shown in Fig. 12 and for phantom test set C in Fig. 12. Descriptive statistics of the relative error distributions in the estimated rCu data are reported in Table 1 for the two phantom test sets B and C.

Fig. 12

Error distributions of the phantom test set B [cf. Fig. 4(b)]. (a) rCu estimation errors $\mathrm{\Delta }{\mathrm{rCu}}_{\mathrm{est}}$ and (b) their absolutes. Blue shows the rCu estimators using MI-LSD, orange the estimators using LSD, and gray is the LU reference. Medians $Q_2$ of the error distributions are shown, together with IQR and 90th percentiles $P_{90}$. The feedforward NN models performed similar to the RF models for the LSD method but underperformed for MI-LSD.

Table 1

Relative rCu estimation errors (ΔrCuest) and absolute rCu estimation errors (|ΔrCuest|) for the RFs, NNs, and linear unmixing. Mean, median Q2, first and third quartiles Q1 and Q3, and the 90th percentile P90 are listed for the phantom test sets B (transversal tubes) and C (longitudinal tubes).

4.

Discussion

The qualitative comparison of the absorbed energy spectra from the Monte Carlo simulations and the phantom OA signal spectra reveals a general agreement between the simulations and the phantom results. The existing variations between the normalized spectra of the two domains are likely due to discrepancies in the simulation, e.g., the beam profiles and the optical properties of the gel pad. The gel pad for example is currently simulated as water but also has some low-level scattering properties, which was omitted in the simulation. In addition, the realistic laser noise was not simulated and the phantom positioning was only accurate to a millimeter. An acoustic forward simulation (e.g., using k-wave) was also not included in the simulation pipeline due to computational time constraints. While there are some acoustic artifacts (e.g., reflection artifacts) in the real OA image reconstructions, it is sensible to assume that they do not vary for different wavelength illumination, therefore, their effect on spectral coloring should be negligible. Variations in the Grüneisen parameter were also not part of the training set, even though it does vary significantly with rCu, because $\mathrm{\Gamma }(c_{\mathrm{wb}}({\mathrm{NiSO}}_4))\approx 0.21$ and $\mathrm{\Gamma }(c_{\mathrm{wb}}({\mathrm{CuSO}}_4))\approx 0.14$ at room temperature.²⁷ This results in a systematically higher SNR for low rCu—an effect not present in ${\mathrm{sO}}_2$,³⁷ which may explain why high rCu estimations are systematically worse in all of our phantom test sets. Laser noise levels are also wavelength dependent, which is reinforced by the pulse energy correction, e.g., resulting in a factor two SNR when measuring at 800 nm compared with 680 nm.

Our MI-LSD method with RF estimators was highly accurate with median absolute estimation errors of only 2.9 and 4.5 pp in the two phantom test sets, respectively. Our NN models, however, failed to give accurate estimates for MI-LSD. LSD estimates using NN were only improved over the LU reference and only in the phantom test set B. When testing on the OOD test set C, our NN models showed no clear improvement over LU. This leads us to the initial conclusion that the overly complex NN models are prone to overfitting to the in silico data, even when optimizing their hyperparameters with simple phantom data. The attempt to remedy this with dropout layers lead to overall inaccurate estimations.

It is not surprising that the overall quantification performance was worse in deeper tubes. SNR in 8 mm deep tubes was very low, i.e., the longer distance illuminations with 980 nm light often yielded no detectable OA signal. This is due to background water absorption in combination with the high scattering, even when adding no sulfates to the background medium. We therefore also investigated omitting these higher wavelengths—training and testing with fewer wavelengths from 680 to 920 nm spaced 20 nm. This yielded obviously worse model performance overall, which either suggests that (these) 13 wavelengths are insufficient for accurate estimation or suggests that spectral coloring due to water absorption can be useful for a pixelwise correction of spectral coloring, as it can give implicit information on the optical path length. For further investigation, it may be useful to add explicit information on the pixel position to the input features. It may also be interesting to perform similar experiments with a wider range of and/or more lower wavelength measurements and then optimize the wavelength selection on these oversampled multispectral sequences. This was not done in this initial proof-of-concept work because it risks overoptimizing on unrealistic aspects of the rCu model (e.g., the difference in Grüneisen parameter of the two sulfate solutions) or setup specific aberrations (e.g., wavelength-dependent SNR). Simulating more wavelengths also prolongs the already computationally expensive, one GPU year, simulation time for the necessary training data (Fig. 13).

Fig. 13

Error distributions of the phantom test set C [cf. Fig. 4(c)]. (a) rCu estimation errors $\mathrm{\Delta }{\mathrm{rCu}}_{\mathrm{est}}$ and (b) their absolutes. Blue shows the rCu estimators using MI-LSD, orange the estimators using LSD, and gray is the LU reference. Medians $Q_2$ of the error distributions are shown, together with IQR and 90th percentiles $P_{90}$. The feedforward NN models performed similar to the RF models for the LSD method but underperformed for MI-LSD.

A final somewhat surprising observation was that estimation of both LSD and to a lesser extent MI-LSD is poorest in phantoms with no added sulfates in the background medium. Figure 14(b) shows the worst estimation results in the lower tubes of phantom test set B—combining three detrimental circumstances: (1) great depth, (2) high rCu, and (3) only spectral coloring of water. Though even in this worst case, MI-LSD is more accurate than the LSD or LU estimations.

Fig. 14

Worst estimation example in phantom test set B: deep tubes (b) compared with more shallow tubes (a) in the same phantom with an $\mathrm{SVF}=0\%$. Estimation results for MI-LSD, LSD, and LU. To indicate the content of the ROI, the mean OA signal in the ROI is shown left, with the ground truth ${\mathrm{rCu}}_{\text{tube}}$ annotated. The brightness of the OA signal is independently and linearly autoleveled for each ROI. The mean rCu estimate $\overline{\mathrm{rCu}}$ over the ROI is noted for the three estimators. Shallow tubes can be estimated very accurately while lower SNR in deep, high rCu tubes correlates to poor estimation accuracy.

One of the main shortcomings of the presented phantom validation is that it did not model melanin absorption in skin. Spectral coloring by melanin still causes large errors in standard of care pulse oximetry devices^38,39 and needs to be addressed for quantitative OA imaging. We were, however, not able to reproducibly include a skin mimicking layer with varying melanin absorption in our liquid phantoms—future work will address this, using solid, layered gel wax phantoms.⁴⁰

We showed a proof-of-concept setup with comparably poor image quality due to the US DAQ and transducer. An in vivo applicable system should make use of state-of-the-art US components and further engineering improvements to sensitivity and SNR, as this currently further limits the achievable estimation accuracy for deep ROI using our setup. Wavelength selection and illumination geometry are suitable but their optimal choice was not the aim of this work. Lastly, while the rCu model is a very useful tool for the investigation and thorough validation of a quantitative OA oximetry method and while the MI-LSD approach shows similar results for in silico ${\mathrm{sO}}_2$, an explicit translation to actual ${\mathrm{sO}}_2$ estimation in vivo must be the next step. One of the main challenges for this translation will be the adequate modeling of additional chromophore distributions, such as melanin. Melanin will strongly affect both overall SNR and spectral coloring.

5.

Conclusions

We presented MI-LSD, a quantitative OA oximetry method using MIs and machine learning; and presented a real-time MI OA imaging setup with a linear ultrasound transducer. We used this setup to image 115 phantom configurations by employing a highly reliable, reproducible, and easily scalable phantom model.

MI-LSD with RFs was able to accurately and quickly estimate blood oxygen saturation modeled by copper and nickel sulfate. Compared with LU, MI-LSD approximately halved the magnitude of the relative estimation error, achieving median absolute estimation errors of only 2.9 and 4.5 pp in our two phantom test sets, respectively. To investigate such ML regression methods, thorough phantom validation is critical, as in silico tests do not give sufficient data to validate a method, and in vivo measurements lack a reliable ground truth. This is further illustrated by the fact that previously reported LSD NN models, which were only validated on in silico data, slightly outperformed RF models on in silico data (as was previously reported) but underperformed RF models in phantom tests while simply breaking on OOD phantom data.

The results of this study give us a high degree of confidence that the domain gap from in silico spectral decoloring to real data can be bridged using MI-LSD, paving the way to a clinical application of quantitative OA oximetry imaging.