Open Access
Add to BibSonomy
Abstract
Optoacoustic (OA) imaging is a promising modality for quantifying blood oxygen saturation (sO2) in various biomedical applications – in diagnosis, monitoring of organ function, or even tumor treatment planning. We present an accurate and practically feasible real-time capable method for quantitative imaging of sO2 based on combining multispectral (MS) and multiple illumination (MI) OA imaging with learned spectral decoloring (LSD). For this purpose we developed a hybrid real-time MI MS OA imaging setup with ultrasound (US) imaging capability; we trained gradient boosting machines on MI spectrally colored absorbed energy spectra generated by generic Monte Carlo simulations and used the trained models to estimate sO2 on real OA measurements. We validated MI-LSD in silico and on in vivo image sequences of radial arteries and accompanying veins of five healthy human volunteers. We compared the performance of the method to prior LSD work and conventional linear unmixing. MI-LSD provided highly accurate results in silico and consistently plausible results in vivo. This preliminary study shows a potentially high applicability of quantitative OA oximetry imaging, using our method.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Optoacoustic (OA) imaging (also named photoacoustic imaging) is a promising biomedical imaging modality for visualization of tissue structure and function based on the optical absorption of endogenous or exogenous chromophores. [1,2]. By combining optical absorption contrast with ultrasonic detection [3] OA imaging is able to image diagnostically relevant chromophores. Endogenous chromophores like lipids, collagen and melanin can be imaged centimeters-deep in tissue for various clinical purposes. Out of these, preclinically and clinically, OA imaging is most frequently investigated for its ability to image blood. It has been used to image inflammation [4,5] and other blood volume correlated biomarkers like angiogenesis [6,7], but also blood oxygen saturation (sO$_2$) [8] to assess organ function [9], wound healing [10], or to follow cancer therapy [11,12].
Despite a wide range of approaches, a persistent challenge in OA imaging is the quantification of actual physiological sO$_2$ deep in tissue [13]. The main obstacle is the ill-posed inverse problem to the optical transport forward problem. Tissue is a turbid medium with unknown distributions of optical absorption and scattering. Light transport through this medium leads to unknown (illumination geometry dependent and wavelength dependent) optical fluence distributions in the imaged regions. This wavelength dependent optical fluence causes so-called spectral coloring [14,15] in the measured OA spectra. Currently machine learning and data driven approaches are frequently investigated in order to move towards quantitative OA oximetry and imaging in general [16–18]. Learned spectral decoloring [19] is one such approach that aims to estimate sO$_2$ pixel-wise from multispectral (MS) OA signals, with each pixel being evaluated individually as a 1D vector – information from other pixels is not taken into account for the estimation. This is intended to break down spectral coloring, which is effected by a large region, to much fewer dimensions. In order to add information and constraints to this ill-posed problem, physical additions to the OA sensing apparatus have been investigated. These include diffuse tomography [20,21] or multiple illumination geometries. [22–25]
Here, we propose that the application of machine learning to local spectra alone is insufficient to perform quantitative OA oximetry – the ill-posed inverse problem remains (too) ill-posed to arrive at accurate or robust estimates of sO$_2$ without further constraints. We therefore propose to incorporate additional information which may allow us to more accurately estimate sO$_2$. We do this by combining multiple illumination (MI), multispectral sensing with a modified LSD technique.
In earlier work [26] we validated such an approach on a highly controlled copper and nickel sulfate [27] phantom model of sO$_2$. In this work we:
- 1. Developed a hybrid real-time MI MS OA imaging setup with ultrasound (US) imaging capability for human applications.
- 2. Designed and performed Monte Carlo simulations of generic human tissue; generating MI spectrally colored absorbed energy spectra.
- 3. Designed and trained gradient boosting machines on these simulations and used the trained models to estimate sO$_2$ on in vivo human OA measurements.
2. Imaging system
We developed a real-time MI OA imaging system with US imaging capability (shown in Fig. 1). The system is based on an US data acquisition system (Vantage128, Verasonics, Inc., Kirkland, Washington, USA) with a linear probe (L7-4, Advanced Technology Laboratories Inc., Bothell, Washington, USA) with 128 transducer elements, 5 MHz center frequency, 0.3 mm pitch and 80% fractional bandwidth. The OA excitation is performed with a fast wavelength-tunable optical parametric oscillator (OPO) laser system (SpitLight, InnoLas Laser GmbH, Krailling, Germany) using 5 ns laser pulses with 100 Hz pulse repetition frequency. Laser pulses are sequentially coupled into four custom fiber bundles (FiberOptic P.+P. AG, Spreitenbach, Switzerland) with NA 0.22 fibers, each bundle 3 mm in diameter. Pulse energy at the fiber exit does not exceed 12 mJ. The sequential coupling is done with a galvo mirror (GVS011/M, Thorlabs Inc., Newton, New Jersey, USA) and a 25 Hz ramp generated by an arbitrary waveform generator (AWG) (TG5011, Aim-TTi, Cambridgeshire, UK) – synchronized with the OPO.
Fig. 1. Hybrid real-time multiple illumination (MI) optoacoustic (OA) imaging setup with ultrasound (US) imaging. Tissue is illuminated through four fiber bundles by a tunable optical parametric oscillator (OPO) laser system firing laser pulses with 100 Hz. A galvo mirror system is driven at 25 Hz by an arbitrary waveform generator (AWG) to illuminate the fiber bundles in sequence. A linear array US probe detects OA signals which are recorded by a 128-channel US data acquisition (DAQ) system. Conventional US imaging and interleaved OA and US imaging can be performed with the same system. A gold layer is applied to the transducer surface to reduce OA transducer artifacts. A durable copolymer-in-oil gel pad is used to enable in plane illumination.
The fiber bundle outputs are arranged in a line array with 8 mm spacing. To comply with American National Standards Institute (ANSI) safety limits for skin [28], the beams are widened to 8 mm full-width at half-maximum (FWHM) at the tissue surface, staying below 25 mJ cm$^{-2}$. The fiber bundle outputs are attached to the handheld linear probe together with a custom made 10 mm thick gel pad through which illumination and acoustic detection are performed. The gel pad is designed to deliver the illumination pulses to the imaging plane. It is based on a recipe for a copolymer in mineral oil optoacoustic phantom by Hacker et al. [29] – omitting Titanium Dioxide (TiO$_2$) and Low-density Polyethylene (LDPE) from the phantom recipe to avoid optical scattering in the material. Compared to water based gel pads it can be molded into shape and fixed to the linear probe much easier while being temporally stable for weeks [29]; it also has a higher optical transmittance for some relevant wavelengths as seen in Fig. 2. The speed of sound in the gel pad was measured to be $\sim$1460 ms$^{-1}$. To reduce the artifact caused by strong optical absorption in the rubber layer of the linear probe, a gold layer (leaf gold $\sim$100 nm thick, 23.75 Karat Blattgold Rosenobel Doppel) was applied directly to the surface of the linear probe.
Fig. 2. Absorption coefficient spectra $\mu _\textrm {a}(\lambda )$ and wavelength selection. Oxy- (HbO) and deoxyhemoglobin (Hb) spectra are shown at whole blood concentrations $c_\textrm {wb}$(Hb or HbO) = 150 gl$^{-1}$ [30]. Absorption spectra of a gel pad sample and the base material mineral oil were acquired by transmission spectrometry. Water absorption is shown for reference. Grey vertical lines indicate the wavelengths selected for the multispectral imaging sequence. The dotted line is 920 nm which was part of the recorded wavelength sequence but only intended for an analysis of systematic errors – see Supplemental Document 1.
The MI-MS imaging sequence consisted of four illumination positions, each illuminated by a sequence of ten laser pulses over a range of wavelengths – (760, 780, 800, 820, 840, 860, 880, 900, 920, 980) nm. While aiming to measure the widest possible range available by the OPO (680 nm to 980 nm), we omitted wavelengths under 760 nm for the following reasons: (1) Variation between pulse energies is highest in the lower tuning range 680–740 nm of the OPO while (2) pulse energy is low, too; (3) Melanin absorption is highest for low wavelengths. Wavelengths between 900 nm and 980 nm were omitted due to the absorption in the gel pad material and the absorption in fat (a similar peak at 930 nm [31]). This allowed us to reduce the complexity of our simulations by neither simulating absorption of the gel pad material nor fat.
In addition to the chosen MS sequence – (760, 780, 800, 820, 840, 860, 880, 900, 980) nm – we measured at an additional wavelength: 920 nm. This was not used for sO$_2$ estimation in the main body of this work, in the supplement it is used to investigate MI-LSD’s robustness to systematic errors in the optical forward model.
One 4 $\times$ 10 pulse sequence is acquired in 0.4 s. For all experiments, we recorded the raw data for 300 such sequences for each subject. Live OA or US beamforming and visualization at 25 fps (showing B-mode US and mean OA signal over all illumination positions) was performed using custom MATLAB scripts. The live visualization was solely used for probe positioning and image quality control (e.g. avoiding air inclusions in the gel). The following data processing was performed on recorded raw data.
3. Optoacoustic data processing
For MI-LSD we train a gradient boosting machine with in silico data and estimate sO$_2$ on in vivo OA data. An overview of these processes is illustrated in Fig. 3. The training pipeline and our sO$_2$ estimation method in general is detailed in subsection 3.1 and the in vivo data processing pipeline in subsection 3.3. Both pipelines are fully open-source (see Data and Code Availability).
Fig. 3. Overview of the multiple illumination learned spectral decoloring (MI-LSD) machine learning method using a Gradient Boosting Machine (GBM) algorithm [32]. Monte Carlo simulations of randomized generic tissue volumes generate absorbed energy spectra $\boldsymbol {H_I}(\boldsymbol {\lambda })$ for four illuminations $\boldsymbol {I}$. The GBM is then trained with tuples of L1 normalized $\boldsymbol {\hat {H}_I}(\boldsymbol {\lambda })$ and labels sO$_2^\textrm {gt}$ extracted from these simulations. Raw optoacoustic data is preprocessed and beamformed with a naive delay and sum (DAS) algorithm. L1 normalized signal spectra from these beamformed images $\boldsymbol {\hat {S}_I}(\boldsymbol {\lambda })$ are then presented to the trained GBM which estimates sO$_2^\textrm {est}$.
3.1 Blood oxygen saturation estimation method
Our method is an extension of learned spectral decoloring (LSD) [19] and is based on our prior work validating multiple illumination LSD (MI-LSD) on copper and nickel sulfate phantoms [26]. In LSD we estimate sO$_2$ in a pixel, for that we perform a regression on a 1D vector – the measured PA signal spectrum $\boldsymbol {S}(\boldsymbol {\lambda })$, with $\boldsymbol {\lambda }$ being the set of measured wavelengths in that pixel. The naive approach to this regression problem is linear spectral unmixing (LU) [33,34]. In LU $\boldsymbol {S}(\boldsymbol {\lambda })$ is fitted to a linear combination of known spectra of known absorbers (i.e. HbO and Hb). In this work, the performed reference LU estimations are based on singular value decomposition using a pseudo inverse absorption matrix. This fast approach was adapted from the simpa toolkit [35]. While LSD aims to estimate sO$_2$ from a single illumination OA signal spectrum $\boldsymbol {S}(\boldsymbol {\lambda })$, MI-LSD simply uses multiple such spectra with a set of illuminations $\boldsymbol {I}$ as inputs. As in prior work both our LSD and MI-LSD models are machine learning algorithms that are trained on simulated absorbed energy spectra $\boldsymbol {H}(\boldsymbol {\lambda })$ labeled with ground truth sO$_2$. These absorbed energy spectra are each normalized with their respective L1 norm resulting into $\boldsymbol {\hat {H}}_I(\boldsymbol {\lambda })$; this is done for each illumination $I$ separately. As shown in previous work [26], this normalization makes them equivalent to normalized measured OA signal spectra $\boldsymbol {\hat {S}}_I(\boldsymbol {\lambda })$.
The MI-LSD model is trained on an in silico training set consisting of tuples of the absorbed energy spectra $\boldsymbol {\hat {H}}_I(\boldsymbol {\lambda })$ for all four illuminations $\boldsymbol {I}$ and the ground truth sO$_2$ label ($\boldsymbol {\hat {H}}_{\boldsymbol {I}}(\boldsymbol {\lambda }$), sO$_2^\textrm {gt}$). The trained model is then presented (1) unseen in silico spectra $\boldsymbol {\hat {H}}_{\boldsymbol {I}}(\boldsymbol {\lambda })$ from a separate test set or (2) OA signal spectra $\boldsymbol {\hat {S}}_{\boldsymbol {I}}(\boldsymbol {\lambda })$ from an unseen in vivo measurement. The trained model then estimates the corresponding sO$_2^\textrm {est}$. For the LSD model, illumination is averaged by $\boldsymbol {\hat {H}}_\textrm {LSD} = \sum _{\boldsymbol {I}} \boldsymbol {\hat {H}}_I$ and $\boldsymbol {\hat {S}}_\textrm {LSD} = \sum _{\boldsymbol {I}} \boldsymbol {\hat {S}}_I$. As machine learning algorithms we used gradient boosting machine regressors – unmodified – as implemented by the LightGBM framework [32]. In initial experiments on data from previous work [26] they proved faster and less memory intensive than random forests, while providing the same precision. Specifically we used a LGBMRegressor with the regression_l1 objective function and 300 estimators (boosting trees) with a maximum of 200 leaves. These were the only tuned hyperparameters. The hyperparameter tuning was performed exclusively on in silico data sets from previous work, and all other hyperparameters were set to default values (also see the open-source code).
3.2 In silico data processing
To generate in silico training and test data sets for the machine learning algorithms we used Monte Carlo simulations of light transport. These simulations were performed with the GPU accelerated open-source mcx toolkit [36] and followed previous simulation studies [26]. The open sourced simpa [35] framework, was used for illumination modeling and data management. The aim of these simulations is not to exactly model the specific investigated tissue geometry like in previous work [19] but to generate a simplified generic data set that contains the broad range of spectral coloring that is expected for transcutaneous OA imaging.
For that we model a broad range of random blood absorbers, background absorption of blood and water, and a randomized epidermis with variable melanin content. The randomization of our tissue layout is shown in Fig. 4. The in silico training data set consists of 4000 volumes and the separate in silico test set of 1000 volumes. Each volume was simulated with ten wavelengths and four illumination positions – mirroring the MI OA imaging sequence with our system.
Fig. 4. Generic tissue geometry used for the Monte Carlo simulations. In order to generate the training set 4000 random tissues are simulated, for the test set 1000 random tissues – each using the shown parameters. Each volume has two sets of tubes, each set with a random number of tubes, uniformly distributed as specified by $\boldsymbol {x}^\textrm {tube}$, and an epidermis layer with thickness $d^\textrm{skin}$. Melanin volume fraction (mvf$^\textrm {skin}$) in the epidermis layer, blood oxygen saturations in the tube set and background tissue (sO$_2^\textrm {bg}$ and sO$_2^\textrm {tube}$) as well as background blood volume fraction (bvf$^\textrm {bg}$) are again drawn from uniform random distributions $U$. The illumination geometry is modeled using the geometry of the setup shown in Fig. 1.
Each volume has two sets of tubes with the tube count drawn from a discrete uniform distribution $U${3, 9}. Each tube has a radius of 0.4 mm, is perpendicular to the imaging plane, and is uniform-randomly placed in the volume as specified in Fig. 4. The sO$_2$ in each tube set and the background sO$_2$ were picked from a continuous uniform random distribution $U$(0%, 100%). All tubes within a tube set have the same random sO$_2$. The background blood volume fraction bvf$_\textrm {bg}$ was picked from $U$(0%, 3%). The wavelength-dependent background scattering parameters were set using an analytic soft tissue scattering approximation published by Jacques [30]. The relevant scattering parameters for background tissue were fixed to average values for soft tissue derived by that meta analysis ($\mu _\textrm {s}'^{500}$ = 19.1, $f^\textrm {ray}$ = 0.153 $b^\textrm {mie}$ = 1.091). For epidermis the scattering was defined using data from Salomatina et al. [37] ($\mu _\textrm {s}'^{500}$ = 66.7, $f^\textrm {ray}$ = 0.29, $b^\textrm {mie}$ = 0.689). The thickness of the epidermis layer was based on Oltulu et al. [38] $U$(1.2 mm, 3.2 mm). The melanin volume fraction (mvf) in the epidermis layer was sampled from $U$(1.2%, 3.2%) roughly modeling the distribution determined by Alaluf et al. [39] over all ethnicities for photoexposed and photoprotected skin. For each single wavelength and single illumination position, the Monte Carlo simulation was performed with $2\times 10^{7}$ photon packets. We used 1080 GTX (NVIDIA, Santa Clara) GPUs on a high performance computing cluster. A single-wavelength and single-illumination-position simulation took an average of $\sim$30 s, so that all simulations for the test and training sets used a combined 70 days of GPU time. The absorbed energy density generated by the simulation was averaged over the thickness ($\sim$2 mm) of the experimental imaging plane to increase in silico SNR per photon packet. This greatly reduced the number of photons needed for a precise simulation, when compared to previous work. The in silico training and test data sets contain only the tube pixels (not the background pixels) and are open data (see the Data and Code Availability Statement).
3.3 In vivo data processing pipeline
The recorded raw OA data was preprocessed, including a correction for the mean OPO pulse energy spectrum. This OPO pulse energy spectrum was measured at the fiber bundles outputs before the experiments. For a single wavelength, the variation of pulse energy was $\sim$3%. For noise and artifact reduction, the data was then bandpass filtered using a Hann window from 0 to 9.6 MHz (corresponding to the frequency range recorded by the DAQ). The beamforming of OA images was performed with a python implementation of delay and sum (DAS) beamforming (using the numba library to increase processing speed). This beamforming implementation was validated against a previous DAS implementation [40]. The delay and sum (DAS) algorithm assumed a fixed speed of sound of 1480 ms$^{-1}$. This value was chosen intentionally lower than the typically assumed tissue speed of sound (1540 ms$^{-1}$), to compensate for the lower speed of sound in the 1 cm thick gel pad material (1460 ms$^{-1}$). Furthermore a Hann apodization was used over an apodization angle of $\pm 30$ degrees. The B-mode images were formed using a Hilbert transform filter for envelope detection and then downsampling the result to an isometric resolution of 0.15 mm. See the Data and Code Availability Statement for the open-source image processing scripts used. Limited by the pulse repetition rate of the laser, one imaging sequence took 0.4 s to measure. On this time scale, tissue motion causes spatial deregistration of the images acquired at different wavelengths, which leads to artifacts in the sO$_2$ estimates. To reduce these unavoidable motion artifacts and also to further increase SNR we averaged our OA B-Mode images with a simple moving average – calculating the mean signal over five adjacent frames for each frame.
4. Experiments and Results
We tested the performance of MI-LSD for oximetry on human subjects and compared it to LSD and LU sO$_2$ estimation performance. For this we performed (1) an in silico validation experiment, (2) an in vivo experiment, imaging radial arteries and accompanying veins of five healthy human volunteers and (3) a hemodynamic challenge on one volunteer.
4.1 In silico
During the in silico validation of the MI-LSD and LSD methods we tested our trained GBM regressors on the in silico test set which was kept separately. We compared our results to LU using the same simulated imaging sequence as for LSD. Results are presented in Fig. 5 and Table 1, for the imaging sequence with the wavelengths (760, 780, 800, 820, 840, 860, 880, 900, 980) nm.
Fig. 5. Estimates of blood oxygenation saturation (sO$_2$) on the in silico test set. Estimated value (sO$_2^\textrm {est}$) with a linear unmixing (LU) and b with multiple illumination learned spectral decoloring (MI-LSD) is plotted against ground truth (sO$_2^\textrm {gt}$) in 2D confusion plots with a logarithmic color scale; in c zoomed in. d shows a histogram over the estimation errors for LU in grey, LSD in light-blue and MI-LSD in blue.
Table 1. The median absolute estimation error $|\Delta \textrm {sO}_2^\textrm {est}|$in silico for MI-LSD, LSD and LU in percentage points (pp) with its interquartile range (IQR) from quartile 1 (Q$_{1}$) to quartile 3 (Q$_{3}$) and its 90th percentile P$_{90}$.
Regressors were also trained and tested on sequences either including 920 nm illumination or excluding 980 nm (water absorption). The validation results for these modified in silico sequences were similar and can be found in the supplemental figures S1 and S2.
The estimation for all 265979 samples in the test set on a high-end consumer CPU (Intel i9-9900KF) was computationally inexpensive with total estimation times of 20 ms for LU, 700 ms for LSD and 800 ms for MI-LSD. Training the GBM regressors on the 1050673 sample training set took 13 s for LSD and 20 s for MI-LSD on the same CPU.
4.2 In vivo baseline
For this experiment we aimed to image easily accessible vasculature with known references for sO$_2$. We therefore imaged the radial artery and surrounding vessels in five subjects. All subjects gave their consent after having been thoroughly informed about the study, and the voluntary nature of participation. Before these measurements on the left forearms of the subjects we took pulse oximeter (OxiMax NPB-40, Nellcor Puritan Bennett Inc., Pleasanton, CA USA) readings (SpO$_2$) from the index finger of the left hand. For each recording of each subject, we imaged a transverse plane of the left forearm eight to twelve centimeters from the wrist. We located the radial artery with the integrated US imaging (see US in Fig. 6) and positioned the probe such that the radial artery was clearly visible in the center third of the US image – saving that live US image for reference. We then switched to a live OA view, confirmed that the artery was visible in the OA image and proceeded with a two minute OA raw data recording of the region of interest. This raw data was processed offline using the pipeline described in section 3.3. We then segmented bounding boxes for the radial artery, an accompanying vein and a superficial vessel on the first of the 290 OA images (see boxes for OA and US in Fig. 6 ). The vessels were then automatically tracked with fixed size ROIs ($5 \times 3$ pixels) by centering on the maximum intensity OA pixel in each segmented bounding box. All estimates are given as averages over these tracked ROIs. Note that the same regressors trained exclusively on in silico data were used for in vivo estimation. The LU, MI-LSD and LSD estimates in the pixels of these ROIs were averaged and are shown as time series in Fig. 6. The time series averages and standard deviations are listed in Table 2 together with reference SpO$_2$ readings. Videos showing estimates over the entire image over time for all estimation methods and subjects are part of the open-data supplement. Interleaved US was also recorded for subject 1 and 2. With each of the methods the sO$_2$ estimation can be performed in real time.
Fig. 6. Blood oxygen saturation (sO$_2$) estimation on optoacoustic (OA) measurements of the radial artery, radial vein and a superficial vessel in five healthy human volunteers. Ultrasound (US) imaging was used to position the probe and then verify the position of the vessels. The arterial, venous and superficial vessel bounding boxes are marked in the US and OS images in red, blue and grey. Estimation was performed with linear unmixing (LU), multiple illumination learned spectral decoloring (MI-LSD) and mean-illumination learned spectral decoloring (LSD), for comparison.
Table 2. Baseline blood oxygen saturation (sO$_2$) estimation on optoacoustic (OA) measurements of the radial artery, radial vein and an additional superficial vessel in five healthy human volunteers. sO$_2$ values are averaged over 290 sequences recorded in two minutes and over a $5\times 3$ pixel region of interest (ROI), automatically tracking the vessels. Full sO$_2$ time series are plotted in Fig. 6. Pulse oximetry (SpO$_2$) measurements on the index finger were taken as a reference before the OA measurement.
4.3 in vivo hemodynamic challenge
As an additional experiment we performed a venous and arterial cuff occlusion on the forearm of subject 1. 250 mmHg pressure was applied to occlude arterial and venous blood flow [41] for ten minutes. We started OA recording 20 seconds before opening the cuff occlusion and expected a return to baseline for arterial blood after two minutes. All image processing was similar to section 4.2; the resulting time series are shown in Fig. 7. A pulse oximeter reference measurement could not be taken during this experiment due to the absence of a pulse in the cuffed arm.
Fig. 7. Reperfusion experiment. Blood flow was blocked for ten minutes by a blood pressure cuff inflated to 250 mmHg. OA imaging was started 20 seconds before reperfusion by deflation of the cuff. Similar to the baseline experiment, blood oxygen saturation (sO$_2$) estimates are shown for the radial artery, accompanying vein and a superficial vessel. Estimation is performed with linear unmixing (LU), multiple illumination learned spectral decoloring (MI-LSD) and fixed illumination learned spectral decoloring (LSD), for comparison.
This reperfusion experiment resulted in sO$_2$ estimates with MI-LSD that are consistent with the physiology of such hemodynamic challenges; with MI-LSD showing an estimated arterial sO$_2$ of $56\% \pm 10\%$ before reperfusion (averaged over the first 12 s) and a quick return to baseline with an arterial sO$_2$ of $91\% \pm 4\%$ (averaged over the last 12 s recorded). The arterial sO$_2$ estimates for LSD were $61\% \pm 19\%$ occluded and $82\% \pm 10\%$ reperfused; and for LU, $28\% \pm 4\%$ occluded and $82\% \pm 5\%$ reperfused. 100 seconds after reperfusion, radial artery and vein sO$_2$ returned to subject 1’s baseline, and so did the post-experiment pulse oximeter reading.
5. Discussion and conclusion
As expected, the MI-LSD estimations of sO$_2$ on the in silico test set were highly accurate, with median absolute estimation errors of 3.8 percentage points, slightly outperforming LSD with $|\Delta \textrm {sO}_2^\textrm {est}|$ = 5.2 pp and clearly outperforming LU with $|\Delta \textrm {sO}_2^\textrm {est}|$ = 13.1 pp. As can be seen in Fig. 5, MI-LSD (same as LSD and LU) performed worst at the boundaries of 0% and 100% sO$_2$. Accuracy of sO$_2$ estimations in the high regime is often medically most relevant as arterial sO$_2$ under 90% is already considered critical and usually chosen as the clinical threshold for hypoxemia [42]. The likely reason for this worse performance at the boundary is the non-existence of training samples outside the 0 – 100% range which punishes estimating values close to the boundary and leads to an overall bias toward 50% sO$_2$.
In silico accuracy was limited by simulation time, but less by the amount of training data and more by the signal-to-noise of the training data. Simulating more photon packets will likely improve the in silico performance but comes at a high computational cost. The MI-LSD model could further be improved by using machine learning algorithms that aim for more explainability [43] of their estimations. This may give us a realistic estimation of accuracy or add confidence estimates to the sO$_2$ estimates of our model [44]. This would of course also increase computational costs. Using lightGBM [32] is accurate, computationally cheap and greatly reduced memory usage compared to the random forests of our previous work [26].
Our in silico estimation accuracy of MI-LSD is high, even if it is only a slight improvement over LSD. However, this testing on in silico data generated with the same random distributions does not imply that MI-LSD or LSD generalizes well. The MI-LSD machine learning model was trained by generic and randomized transcutaneous tissue data provided by Monte Carlo simulations. These random data aim to simulate the full range of possible spectral coloring in and by blood and water and include spectral coloring by a wide range of melanin concentration in skin. The aim was to include the full ethnic variation in melanin. Previous quantitative OA oximetry work [19], when applied to in vivo or phantom test data, used in silico training data which was not generic but closely modeled after the relevant in vivo or phantom test data. In our in vivo experiments we estimated sO$_2$ with machine learning models trained with dissimilar and purposely generic in silico data. Still, our in vivo MI-LSD estimates of sO$_2$ yielded plausible and highly consistent results. The accuracy of these estimates can not be judged well, as these experiments on healthy human volunteers lack hard ground truth data for sO$_2$ which would only be available with invasive means (i.e. arterial blood gas analysis). We employed pulse oximetry to get a non-invasive reference measure for the true arterial sO$_2$. As an additional reference we used literature values for normal arterial sO$_2$ averaging around 97% [45].
All our arterial blood in vivo MI-LSD estimation results for radial arteries are consistent with both literature and pulse oximetry, except for subject 5’s. LU and LSD deviate from the arterial sO$_2$ reference values for all subjects (see Table 2). Arterial sO$_2$ baseline estimates for MI-LSD are systematically $\sim$2–5% lower than expected, which is in line with the deviation observed in silico. For venous blood we only have literature values as references. In central venous blood returning from peripheral tissues sO$_2$ averages around 75% [45], but with wide variation [46]. Normal venous sO$_2$ in the arm is usually found to average around 70% with a large standard deviation of 20% [47] or 12% [48]. All our in vivo MI-LSD results are consistent with this plausible range for venous blood. LU and LSD have several outliers in their venous sO$_2$ estimates (see Table 2 and Fig. 6).
It is unclear why subject 5 is an outlier in the arterial sO$_2$ estimation using MI-LSD. The falling trendline (see Fig. 6) could indicate an inadvertent occlusion (maybe partial arterial occlusion) due to pressure applied by the imaging probe. SpO$_2$ was not continuously measured and would have indicated if that were the case. In general, the unequal and uneven pressure applied by the probe to the forearm may also have resulted in a (partial) collapse of the radial veins in some experiments. This may have led to variation over time in sO$_2$ for some of the baselines and should be closely monitored in future work.
Generally, motion artifacts also have an additional (blurring) influence on the data. It should also be noted that the acoustic reconstruction with DAS and using a linear transducer array is an imperfect solution to the acoustic inverse problem. This means that we do not in fact estimate from a fully local OA signal spectrum but our OA signal spectra include confounding influences from the surrounding area. A major limitation of our pilot study is the light skin type in all our subjects. All subjects were light skin males, though with varying degrees of skin photoexposure. The group only included relatively light skin on the ventral side of their forearms (Type 1–3 on the Fitzpatrick scale). Therefore, this pilot study, while promising, yields only limited information about the method’s usefulness with the full range of ethnic variation in melanin. A follow-up study of MI-LSD OA oximetry should include continuous reference measurements with pulse oximetry, possibly a blood gas analysis of venous blood samples and definitely a larger, more diverse cohort.
The result of our hemodynamic challenge experiment showed reduced arterial sO$_2$ estimates for all methods after a ten minute cuff occlusion. Following reperfusion we could observe a characteristic overshoot (hyperemia) of arterial sO$_2$ followed by a return to the baseline of subject 1. Again the MI-LSD estimates followed the expected hemodynamics well. LU and LSD estimation performed worse. Maybe more importantly, we can see from this experiment that MI-LSD does not simply systematically estimate higher sO$_2$ values but indeed estimates plausible values in a wide range. In a follow-up study of MI-LSD OA oximetry such a hemodynamic challenge experiment could be extended by an oxygen challenge, which would allow a pulse oximetry reference measurement. Other experiments may include a local variation of perfusion e.g. with adrenaline injections [49].
In conclusion, MI-LSD proved highly accurate in silico and showed consistently plausible estimates of sO$_2$ in vivo. Our preliminary study shows that MI-LSD has the potential to be a robust tool enabling quantitative OA oximetry imaging in humans.
Funding
Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (205320-179038); Deutsche Forschungsgemeinschaft (471755457).
Acknowledgments
Calculations for the Monte Carlo simulations were performed on UBELIX, the HPC cluster at the University of Bern. We thank Adrian Jenk at the Institute of Applied Physics mechanical workshop. For their support of the open-source SIMPA framework we thank Janek Groehl and the Photoacoustics team at the Computer Assisted Medical Interventions division, German Cancer Research Center, Heidelberg.
Author contributions: Conceptualization, T.K. and M.F.; Software, T.K. and M.J.; Methodology, Investigation, Formal Analysis, Data curation, Validation, Visualization, Writing – original draft, T.K.; Writing – review & editing, T.K., M.J. and M.F.; Supervision, Funding acquisition, M.F.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are openly available at [50]. The code for the methods described in this paper is fully open-source; available at [51].
Supplemental document
See Supplement 1 for supporting content.
References
1. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 041101 (2006). [CrossRef]
2. J. Xia, J. Yao, and L. V. Wang, “Photoacoustic tomography: principles and advances,” Electromagnetic waves (Cambridge, Mass.) 147, 1–22 (2014).
3. J. Laufer, “Photoacoustic imaging: Principles and applications,” in Quantification of Biophysical Parameters in Medical Imaging, (Springer, Cham, 2018), pp. 303–324.
4. F. Knieling, C. Neufert, A. Hartmann, J. Claussen, A. Urich, C. Egger, M. Vetter, S. Fischer, L. Pfeifer, A. Hagel, C. Kielisch, R. S. Görtz, D. Wildner, M. Engel, J. Röther, W. Uter, J. Siebler, R. Atreya, W. Rascher, D. Strobel, M. F. Neurath, and M. J. Waldner, “Multispectral Optoacoustic Tomography for Assessment of Crohn’s Disease Activity,” N. Engl. J. Med. 376(13), 1292–1294 (2017). [CrossRef]
5. J. Jo, C. Tian, G. Xu, J. Sarazin, E. Schiopu, G. Gandikota, and X. Wang, “Photoacoustic tomography for human musculoskeletal imaging and inflammatory arthritis detection,” Photoacoustics 12, 82–89 (2018). [CrossRef]
6. R. I. Siphanto, K. K. Thumma, R. G. M. Kolkman, T. G. v. Leeuwen, F. F. M. d. Mul, J. W. v. Neck, L. N. A. v. Adrichem, and W. Steenbergen, “Serial noninvasive photoacoustic imaging of neovascularization in tumor angiogenesis,” Opt. Express 13(1), 89–95 (2005). [CrossRef]
7. A. Horiguchi, M. Shinchi, A. Nakamura, T. Wada, K. Ito, T. Asano, H. Shinmoto, H. Tsuda, and M. Ishihara, “Pilot study of prostate cancer angiogenesis imaging using a photoacoustic imaging system,” Urology 108, 212–219 (2017). [CrossRef]
8. M. Li, Y. Tang, and J. Yao, “Photoacoustic tomography of blood oxygenation: A mini review,” Photoacoustics 10, 65–73 (2018). [CrossRef]
9. T. Kirchner, J. Gröhl, M. A. Herrera, T. Adler, A. Hernández-Aguilera, E. Santos, and L. Maier-Hein, “Photoacoustics can image spreading depolarization deep in gyrencephalic brain,” Sci. Rep. 9(1), 8661 (2019). [CrossRef]
10. K. Aizawa, S. Sato, D. Saitoh, H. Ashida, and M. Obara, “Photoacoustic monitoring of burn healing process in rats,” J. Biomed. Opt. 13(6), 064020 (2008). [CrossRef]
11. S. Mallidi, G. P. Luke, and S. Emelianov, “Photoacoustic imaging in cancer detection, diagnosis, and treatment guidance,” Trends Biotechnol. 29(5), 213–221 (2011). [CrossRef]
12. J. L. Su, R. R. Bouchard, A. B. Karpiouk, J. D. Hazle, and S. Y. Emelianov, “Photoacoustic imaging of prostate brachytherapy seeds,” Biomed. Opt. Express 2(8), 2243–2254 (2011). [CrossRef]
13. B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, “Quantitative spectroscopic photoacoustic imaging: a review,” J. Biomed. Opt. 17(6), 061202 (2012). [CrossRef]
14. R. Hochuli, L. An, P. C. Beard, and B. T. Cox, “Estimating blood oxygenation from photoacoustic images: can a simple linear spectroscopic inversion ever work?” J. Biomed. Opt. 24(12), 1 (2019). [CrossRef]
15. K. Maslov, H. F. Zhang, and L. V. Wang, “Effects of wavelength-dependent fluence attenuation on the noninvasive photoacoustic imaging of hemoglobin oxygen saturation in subcutaneous vasculature in vivo,” Inverse Probl. 23(6), S113–S122 (2007). [CrossRef]
16. J. Gröhl, M. Schellenberg, K. Dreher, and L. Maier-Hein, “Deep learning for biomedical photoacoustic imaging: A review,” Photoacoustics 22, 100241 (2021). [CrossRef]
17. C. Bench, A. Hauptmann, and B. T. Cox, “Toward accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in three dimensions,” J. Biomed. Opt. 25(08), 085003 (2020). [CrossRef]
18. S. Tzoumas, A. Nunes, I. Olefir, S. Stangl, P. Symvoulidis, S. Glasl, C. Bayer, G. Multhoff, and V. Ntziachristos, “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nat. Commun. 7(1), 12121 (2016). [CrossRef]
19. J. Gröhl, T. Kirchner, T. J. Adler, L. Hacker, N. Holzwarth, A. Hernández-Aguilera, M. A. Herrera, E. Santos, S. E. Bohndiek, and L. Maier-Hein, “Learned spectral decoloring enables photoacoustic oximetry,” Sci. Rep. 11(1), 6565 (2021). [CrossRef]
20. A. Q. Bauer, R. E. Nothdurft, T. N. Erpelding, L. V. Wang, and J. P. Culver, “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16(9), 096016 (2011). [CrossRef]
21. L. Ulrich, L. Ahnen, H. G. Akarçay, S. S. Majos, M. Jaeger, K. G. Held, M. Wolf, and M. Frenz, “Spectral correction for handheld optoacoustic imaging by means of near-infrared optical tomography in reflection mode,” J. Biophotonics 12(1), e201800112 (2019). [CrossRef]
22. R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. 49(18), 3566–3572 (2010). [CrossRef]
23. P. Shao, B. Cox, and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt. 50(19), 3145–3154 (2011). [CrossRef]
24. K. G. Held, M. Jaeger, J. Ricka, M. Frenz, and H. G. Akarçay, “Multiple irradiation sensing of the optical effective attenuation coefficient for spectral correction in handheld OA imaging,” Photoacoustics 4(2), 70–80 (2016). [CrossRef]
25. M. Kim, G.-S. Jeng, M. O’Donnell, and I. Pelivanov, “Correction of wavelength-dependent laser fluence in swept-beam spectroscopic photoacoustic imaging with a hand-held probe,” Photoacoustics 19, 100192 (2020). [CrossRef]
26. T. Kirchner and M. Frenz, “Multiple illumination learned spectral decoloring for quantitative optoacoustic oximetry imaging,” J. Biomed. Opt. 26(08), 085001 (2021). [CrossRef]
27. M. B. Fonseca, L. An, and B. T. Cox, “Sulfates as chromophores for multiwavelength photoacoustic imaging phantoms,” J. Biomed. Opt. 22(12), 125007 (2017). [CrossRef]
28. ANSI, American National Standard for Safe Use of Lasers Z136.1 (Laser Institute of America, 2014).
29. L. Hacker, J. Joseph, A. M. Ivory, M. O. Saed, B. Zeqiri, S. Rajagopal, and S. E. Bohndiek, “A Copolymer-in-Oil Tissue-Mimicking Material With Tuneable Acoustic and Optical Characteristics for Photoacoustic Imaging Phantoms,” IEEE Trans. Med. Imaging 40(12), 3593–3603 (2021). [CrossRef]
30. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]
31. R. L. P. v. Veen, H. J. C. M. Sterenborg, A. Pifferi, A. Torricelli, and R. Cubeddu, “Determination of VIS- NIR absorption coefficients of mammalian fat, with time- and spatially resolved diffuse reflectance and transmission spectroscopy,” in Biomedical Topical Meeting, (Optical Society of America, 2004), p. SF4.
32. G. Ke, Q. Meng, T. Finley, T. Wang, W. Chen, W. Ma, Q. Ye, and T.-Y. Liu, “LightGBM: A Highly Efficient Gradient Boosting Decision Tree,” in Advances in Neural Information Processing Systems, vol. 30 (Curran Associates, Inc., 2017).
33. S. Tzoumas and V. Ntziachristos, “Spectral unmixing techniques for optoacoustic imaging of tissue pathophysiology,” Philos. Trans. R. Soc., A 375(2107), 20170262 (2017). [CrossRef]
34. T. Kirchner, J. Gröhl, F. Sattler, M. S. Bischoff, A. Laha, M. Nolden, and L. Maier-Hein, “An open-source software platform for translational photoacoustic research and its application to motion-corrected blood oxygenation estimation,” arXiv preprint arXiv:1901.09781 (2019).
35. J. Gröhl, K. K. Dreher, M. Schellenberg, A. Seitel, and L. Maier-Hein, “SIMPA: an open source toolkit for simulation and processing of photoacoustic images,” in Photons Plus Ultrasound: Imaging and Sensing 2021, vol. 11642 (SPIE, 2021), p. 116423C.
36. Q. Fang and D. A. Boas, “Monte Carlo Simulation of Photon Migration in 3D Turbid Media Accelerated by Graphics Processing Units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]
37. E. V. Salomatina, B. Jiang, J. Novak, and A. N. Yaroslavsky, “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11(6), 064026 (2006). [CrossRef]
38. P. Oltulu, B. Ince, N. Kokbudak, S. Findik, and F. Kilinc, “Measurement of epidermis, dermis, and total skin thicknesses from six different body regions with a new ethical histometric technique,” Turk J Plast Surg 26(2), 56 (2018). [CrossRef]
39. S. Alaluf, D. Atkins, K. Barrett, M. Blount, N. Carter, and A. Heath, “Ethnic Variation in Melanin Content and Composition in Photoexposed and Photoprotected Human Skin,” Pigm. Cell Res. 15(2), 112–118 (2002). [CrossRef]
40. T. Kirchner, F. Sattler, J. Gröhl, L. Maier-Hein, T. Kirchner, F. Sattler, J. Gröhl, and L. Maier-Hein, “Signed real-time delay multiply and sum beamforming for multispectral photoacoustic imaging,” J. Imaging 4(10), 121 (2018). [CrossRef]
41. A. Sudakou, F. Lange, H. Isler, P. Lanka, S. Wojtkiewicz, P. Sawosz, D. Ostojic, M. Wolf, A. Pifferi, I. Tachtsidis, A. Liebert, and A. Gerega, “Time-domain NIRS system based on supercontinuum light source and multi-wavelength detection: validation for tissue oxygenation studies,” Biomed. Opt. Express 12(10), 6629–6650 (2021). [CrossRef]
42. I. L. Tyler, B. Tantisira, P. M. Winter, and E. K. Motoyama, “Continuous monitoring of arterial oxygen saturation with pulse oximetry during transfer to the recovery room,” Anesth. Analg. (Philadelphia, PA, U. S.) 64(11), 1108–1112 (1985). [CrossRef]
43. A. Holzinger, “From Machine Learning to Explainable AI,” in 2018 World Symposium on Digital Intelligence for Systems and Machines (DISA), (2018), pp. 55–66.
44. J. Gröhl, T. Kirchner, T. Adler, and L. Maier-Hein, “Confidence estimation for machine learning-based quantitative photoacoustics,” J. Imaging 4(12), 147 (2018). [CrossRef]
45. J. E. Hall, Guyton and Hall Textbook of Medical Physiology (Elsevier, 2015).
46. K. Reinhart and F. Bloos, “The value of venous oximetry,” Current Opinion in Critical Care 11(3), 259–263 (2005). [CrossRef]
47. A. M. Belhaj, J. P. Phillips, P. A. Kyriacou, and R. M. Langford, “Comparison of non-invasive peripheral venous saturations with venous blood co-oximetry,” J. Clin. Monit. Comput. 31(6), 1213–1220 (2017). [CrossRef]
48. A. Keys, “The oxygen saturation of the venous blood in normal human subjects,” Am. J. Physiol. 124(1), 13–21 (1938). [CrossRef]
49. J. Bunke, A. Merdasa, A. Merdasa, R. Sheikh, J. Albinsson, T. Erlöv, B. Gesslein, M. Cinthio, N. Reistad, and M. Malmsjö, “Photoacoustic imaging for the monitoring of local changes in oxygen saturation following an adrenaline injection in human forearm skin,” Biomed. Opt. Express 12(7), 4084–4096 (2021). [CrossRef]
50. T. Kirchner, M. Jaeger, and M. Frenz, “Machine learning enabled multiple illumination quantitative optoacoustic oximetry imaging in humans - Data,” Zenodo, 2022, https://doi.org/10.5281/zenodo.5929161.
51. T. Kirchner, “MI-LSD-in-vivo,” Github, 2022, https://github.com/thkirchner/MI-LSD-in-vivo.git.
