figure table

# Distribution-informed and wavelength-flexible data-driven photoacoustic
oximetry

Janek Gröhl  Cancer Research UK Cambridge Institute, University of Cambridge,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK
Kylie Yeung  Cancer Research UK Cambridge Institute, University of Cambridge,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK
Kevin Gu  Cancer Research UK Cambridge Institute, University of Cambridge,
Cambridge, UK  Department of Physics, University of Cambridge, Cambridge, UK
Thomas R. Else  Cancer Research UK Cambridge Institute, University of
Cambridge, Cambridge, UK  Department of Physics, University of Cambridge,
Cambridge, UK  Monika Golinska  Cancer Research UK Cambridge Institute,
University of Cambridge, Cambridge, UK  Department of Physics, University of
Cambridge, Cambridge, UK  Ellie V. Bunce  Cancer Research UK Cambridge
Institute, University of Cambridge, Cambridge, UK  Department of Physics,
University of Cambridge, Cambridge, UK  Lina Hacker  Department of Oncology,
University of Oxford, Oxford, U.K.  Sarah E. Bohndiek  Cancer Research UK
Cambridge Institute, University of Cambridge, Cambridge, UK  Department of
Physics, University of Cambridge, Cambridge, UK

###### Abstract

Significance: Photoacoustic imaging (PAI) promises to measure spatially-
resolved blood oxygen saturation, but suffers from a lack of accurate and
robust spectral unmixing methods to deliver on this promise. Accurate blood
oxygenation estimation could have important clinical applications, from cancer
detection to quantifying inflammation.

Aim: This study addresses the inflexibility of existing data-driven methods
for estimating blood oxygenation in PAI by introducing a recurrent neural
network architecture.

Approach: We created 25 simulated training dataset variations to assess neural
network performance. We used a long short-term memory network to implement a
wavelength-flexible network architecture and proposed the Jensen-Shannon
divergence to predict the most suitable training dataset.

Results: The network architecture can handle arbitrary input wavelengths and
outperforms linear unmixing and the previously proposed learned spectral
decolouring method. Small changes in the training data significantly affect
the accuracy of our method, but we find that the Jensen-Shannon divergence
correlates with the estimation error and is thus suitable for predicting the
most appropriate training datasets for any given application.

Conclusions: A flexible data-driven network architecture combined with the
Jensen-Shannon Divergence to predict the best training data set provides a
promising direction that might enable robust data-driven photoacoustic
oximetry for clinical use cases.

###### keywords:

quantitative imaging, oximetry, deep learning, image processing, simulation

✉ Send correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Blood oxygen saturation (sO2) is an important indicator of individual health
status used routinely in patient management [1]. Photoacoustic imaging (PAI)
is a promising medical imaging modality for real-time non-invasive spatially-
resolved measurement of sO2 [2], with early clinical applications [3] shown
in, for example, inflammatory bowel disease [4] and cardiovascular diseases
[5]. In cancer, alterations in localized sO2 levels have been linked to
angiogenesis and hypoxia [6]; key hallmarks of cancer that are known to affect
treatment outcomes [7] and are measurable through PAI.

Unfortunately, it remains difficult to apply PAI to derive quantitative values
for sO2 from multispectral PA measurements [8, 9]. Linear unmixing (LU)
remains the de facto standard for sO2 estimation from PAI measurements [10]
because of its simplicity and flexibility, but has well-understood limitations
in its applicability and accuracy [11]. The limitations of LU are significant
in the context of artefacts arising from: the optical processes (e.g. non-
linear light fluence distribution [12] leading to spectral colouring),
acoustic processes (e.g. reflection artefacts [13]), or reconstruction
algorithms (e.g. model mismatch in sound speed [14]). Using LU is particularly
challenging in the presence of highly absorbing tissues, such as the
epidermis, [15] which can introduce reflection artefacts and a spectral bias
that leads to an overestimation of sO2, which increases with darker skin tone
[16].

Data-driven unmixing schemes have shown promise to alleviate some of the
shortcomings of LU [17, 18, 19] but suffer from three major drawbacks: (1)
inflexibility to receiving different input data after training [20]; (2)
performance determined by the composition of the training dataset [21]; and
(3) limited testing on diverse and representative use cases [22]. In
comparison to LU, data-driven methods are often inflexible regarding the input
data after training, impacting generalizability, making them difficult to use,
and requiring laborious tailoring to a specific application and imaging
system. Thus, data-driven sO2 estimation methods can struggle to translate
from promising findings in silico to in vivo data [23, 24]. The lack of high-
quality annotated training data and reliable validation data has made it
difficult to implement data-driven methods robustly. One way of tackling this
challenge lies in bridging the gap between simulated and actual PA data,
exploiting realistic phantoms [25], an approach that has recently started to
be explored [26, 27, 28]. Further, many data-driven approaches use single-
pixel input spectra for their inversion, as with LU, even though 3D inversions
would be preferred for realistic use cases [29]. Solving the full 3D problem
is typically computationally intensive, limiting its success in vivo [30]. To
make the inverse problem more tractable without full 3D information, some
approaches use priors in the inversion scheme [31], differential image
analysis [9], or multiple measurements with differences in illumination [32].

In this work, we set out to address the three aforementioned limitations. We
improve the flexibility of data-driven sO2 estimation using a Long Short Term
Memory (LSTM) network that enables input wavelength flexibility. We propose a
method to inform the choice of training data, which can either be used to
choose the best pre-trained model for a target application or to inform the
choice of simulation parameters when creating a training data set to underpin
a new model. We test these methods on diverse data sources across simulations,
phantoms, small animals and humans.

##  2 Materials and Methods

We begin by investigating sensitivity to changes in training dataset
simulation parameters, by defining a baseline dataset set (BASE) with typical
assumptions on the tissue geometry and functional parameter ranges, then adapt
it into 25 variations (Fig 1 A). We use the Jensen-Shannon divergence to
determine the ideal training dataset for a given use case. We propose a deep
learning network architecture based on a long short-term memory (LSTM) network
that is flexible regarding the input wavelengths (Fig 1 B) and use a testing
strategy that comprises computational studies in silico, phantoms in gello
[33], and in vivo data (Fig 1 C).

Figure 1: Overview of the methods used. A Photon transport was simulated with
a Monte Carlo light model for each of 25 distinct datasets adapted from a
baseline tissue assumption (BASE). Spectra with simulated amplitude for each
wavelength were extracted from these simulations. B A deep learning network
based on a long short-term memory (LSTM) network was introduced to enable
greater flexibility regarding input wavelengths for analysis. The hidden state
of the LSTM was passed to fully connected layers, which output the estimated
blood oxygenation sO2. C The performance of the LSTM-based method when trained
on datasets with different tissue simulation parameters was tested across
different datasets, ranging from in silico simulations, in gello phantom
measurements, to in vivo measurements. This figure was created with Inkscape
using BioRender assets.

###  2.1 In silico datasets

25 in silico datasets were simulated in Python using the SIMPA toolkit [34]
(Fig1 A). A Monte Carlo model (MCX [35]) was used for the optical forward
model with a 50 mJ Gaussian beam using $10^{7}$ photons and a 20 mm radius,
simulating at 41 wavelengths from 700 nm to 900 nm in 5 nm steps, and assuming
an anisotropy of 0.9. For the datasets that include acoustic modelling, a 2D
k-space pseudo-spectral time domain method implemented in the k-Wave toolbox
[36] was used. We assumed a uniform speed of sound of 1500 ms-1, a density of
1000 gcm-3, and disregarded acoustic absorption. For most datasets, a generic
linear ultrasound detector array was placed in the centre of the simulated
volume. The generic array consists of 100 detection elements with a pitch of
0.18 mm, a centre frequency of 4 MHz, a bandwidth of 55%, and a sampling rate
of 40 MHz. For the datasets mimicking the setup of commercial instruments, we
used the built-in device definitions provided in SIMPA.

Baseline dataset parameters. The parameters of the baseline dataset (BASE)
were chosen to reflect the typical parameter choices found in literature [20,
37, 38, 29, 39, 40]. A 19.2 x 19.2 x 19.2 mm cube was simulated with a 0.3 mm
voxel size, resulting in $64^{3}$ voxels per volume. The background started
from the second voxel from the top and was modelled as muscle tissue with a
blood volume fraction of 1% [20] with 70% oxygenation [41]. We added 0-9
cylinders with a diameter randomly chosen from U[0.3 mm, 2.0 mm], where U
denotes a uniform distribution. The tissue was divided into 3 x 3 equal-sized
compartments in the imaging plane and a vessel was included with a 50%
probability. Each vessel purely contained blood with a randomly drawn sO2
value from U[0%, 100%]. To prevent vessel overlap, they were located within
their respective compartments’ boundaries.

We defined 24 variations of BASE (Table 1). The variations encompassed a range
of tissue backgrounds, vessel sizes, illumination geometries and resolutions,
as well as the addition of a skin layer, performing acoustic simulation, and
using digital twins of two commercial instruments (MSOT Acuity Echo and MSOT
InVision 256-TF; iThera Medical GmbH, Munich, Germany). The MSOT Acuity Echo
is a handheld clinical PAI device with 256 detection elements and an angular
coverage of 120∘, whereas the MSOT InVision 256-TF is a tomographic pre-
clinical PAI device with 256 detection elements and a 270∘ view angle. The
full details for the digital twin parameters of these devices can be found in
prior work [34, 26] and are available within the SIMPA toolkit.

Dataset identifier |  Changes from BASE
---|---
BG_0-100 |  Background sO2 randomly chosen from U(0%, 100%).
BG_60-80 |  Background sO2 randomly chosen from U(60%, 80%).
BG_H2O |  Background modelled as water only.
HET_0-100 |  Background heterogeneously varied between 0% and 100%.
HET_60-80 |  Background heterogeneously varied between 60% and 80%.
RES_0.15 |  Simulation grid spacing: 0.15 mm.
RES_0.15_SMALL |  Simulation grid spacing: 0.15 mm; radii of vessels halved.
RES_0.6 |  Simulation grid spacing: 0.6 mm.
RES_1.2 |  Simulation grid spacing: 1.2 mm.
SKIN |  0.3 mm melanin layer; melanosome fraction: U(0.1%, 5%).
ILLUM_5mm |  Radius of incident beam: 5 mm.
ILLUM_POINT |  Radius of incident beam: 0 mm.
SMALL |  Radii of the vessels halved.
ACOUS |  Acoustic modelling with linear transducer array.
WATER_2cm |  2cm H2O layer added between illumination and tissue.
WATER_4cm |  4cm H2O layer added between illumination and tissue.
MSOT |  Volume extended: 75 x 19.2 x 19.2 mm; MSOT Acuity twin.
MSOT_SKIN |  MSOT + SKIN.
MSOT_ACOUS |  MSOT + ACOUS.
MSOT_ACOUS_SKIN |  MSOT + SKIN + ACOUS.
INVIS |  Volume extended: 90 x 25 x 90 mm; MSOT InVision 256-tf twin. Single vessel in a 10 mm radius tubular background.
INVIS_SKIN |  INVIS + SKIN
INVIS_ACOUS |  INVIS + ACOUS.
INVIS_SKIN_ACOUS |  INVIS + SKIN + ACOUS.
Table 1: Summary of datasets used in the study. Each row identifies the
changes in dataset creation parameters performed for each of the training
datasets relative to the BASE dataset (specified in Methods subsection 2.1).
The left column shows an identifier assigned to each dataset that is used
throughout the manuscript and the right column summarises the deviation from
BASE. U denotes a uniform distribution.

Data preprocessing. We extracted pixel-wise PA spectra from simulations of
initial pressure or reconstructed signal amplitude if acoustic simulations
were performed. We only extracted spectra from blood vessels where the signal
intensity at 800 nm was at least 10% of the maximum signal intensity in the
dataset. If less than 10% of voxels were chosen this way, we selected the 10%
of voxels with the highest signal amplitude from the dataset. We enforced this
selection criterion to effectively exclude voxels with low signal-to-noise
ratio (SNR) caused by optical attenuation, as previous works have shown that
idealised simulation can contain spectra where the spectral colouring is
unrealistically strong in regions with low signal amplitude [20, 32].

The number of extracted spectra from the different dataset variations ranged
from 25 thousand to 60 million, compared to BASE with 7 million spectra; the
mean over all datasets was just above 6 million. The large difference in
extractable spectra is primarily caused by two factors: (1) typical
simulations yield 3D p0 distributions, but adding acoustic forward modelling
leads to a 2D reconstructed image; and (2) some datasets only include a single
vessel in the centre of the phantom tube, mimicking a blood flow phantom [42]
(described below). To mitigate performance differences caused by this
discrepancy, we stratified the dataset sizes by randomly sampling 300,000
spectra with replacement per dataset, which represents a balanced compromise
between undersampling larger datasets and oversampling smaller ones [43].

We performed z-score normalization on each spectrum, setting the mean ($\mu$)
to 0 and variance ($\sigma^{2}$) to 1, which discards signal intensity
information and eliminates the need for quantitative simulation calibration.
Since we performed the same spectra-wise z-score normalisation on experimental
data, this normalisation allows us to apply sO2 estimation algorithms trained
in silico to experimental in gello and in vivo data.

###  2.2 Deep Learning Algorithm

To address the limited flexibility of data-driven oximetry methods [40, 20], a
custom unmixing network architecture was composed that contains a Long Short-
Term Memory (LSTM) network. Due to the recurrent nature of an LSTM, it can
process sparse spectra containing zeros at arbitrary positions (see Fig 1 B).

The network input size was fixed at 41, representing the maximum number of
wavelengths (between 700 nm and 900 nm in 5 nm steps) we consider during
inference. The number is a trade-off between maximising the spectral
resolution of the input features and simulation time efficiency. The network
started with an LSTM layer with a hidden size of 100. A masking layer was
introduced to identify zeroed values, creating a mask that instructs the layer
to ignore missing values. The LSTM output was then flattened into a fixed-
length encoding. Following the LSTM, a three-layer fully connected neural
network (FCNN) was used with input size 4200, hidden size 1000, and output
size 1. A leaky rectified linear unit was used after each layer and the
activation function after the final layer was a sigmoid function to constrain
sO2 predictions between 0 and 1 (Fig 1 B).

The deep learning networks were trained for 100 epochs, where each epoch
included the entire training set. The parameters were optimised with the Adam
optimizer from the Keras framework [44] using an initial learning rate of 10-3
and a mean absolute error loss. The learning rate was halved upon a 5-epoch
plateau of the validation loss.

###  2.3 In gello blood flow phantom imaging

Two variable oxygenation blood flow phantoms were imaged using a previously
described protocol [42]. Briefly, agar-based cylindrical phantoms with a
radius of 10 mm were created and a polyvinyl chloride tube (inner diameter 0.8
mm, outer diameter 0.9 mm) was embedded in the centre before the agar was
allowed to set at room temperature. The base mixture comprised 1.5 % (w/v)
agar (Sigma Aldrich 9012-36-6) in Milli-Q water and was heated until
dissolved. After cooling to approximately 40∘ C, 2.08% (v/v) of pre-warmed
intralipid (20% emulsion, Merck, 68890-65-3) and 0.74% (v/v) Nigrosin solution
(0.5 mg/mL in Milli-Q water) were added, mixed, and poured into the mold.
Imaging was performed using the MSOT InVision 256TF (iThera Medical GmbH,
Munich, Germany) according to a previously established standard operating
procedure [45]. PA images of the phantom were acquired in the range between
700 nm and 900 nm with 20 nm increments. The first phantom was imaged using
deuterium oxide (D2O, heavy water) as the coupling medium within the system,
while the second was immersed in normal water (H2O) during imaging.

###  2.4 In vivo human forearm imaging

Human forearm imaging was performed as part of the PAISKINTONE study, which
started in June 2023 following approval by the East of England - Cambridge
South Research Committee (Ref: 23/EE/0019). The study was conducted in
accordance with the Declaration of Helsinki and written informed consent was
obtained from all study participants. Participants were excluded if they could
not give consent, were under the age of 20 or over 80, or had a body mass
index (BMI) outside the range between 18.5 and 30. Imaging was performed using
the MSOT Acuity Echo (iThera Medical GmbH, Munich, Germany) using laser light
between 660 nm and 1300 nm, averaging over 10 scans each and analysis was
performed at five wavelengths (700 nm, 730 nm, 760 nm, 800 nm, 850 nm). One
forearm scan from N=7 randomly chosen subjects with Fitzpatrick Type 1 or 2
were selected for the purposes of testing the method proposed in this work.
The authors manually segmented the radial artery in each scan using the
medical imaging interaction toolkit (MITK) [46].

###  2.5 In vivo mouse imaging

All animal procedures were conducted under project and personal licences (PPL
no PE12C2B96, PIL no I53057080), issued under the United Kingdom Animals
(Scientific Procedures) Act, 1986 and compliance was approved locally by the
CRUK Cambridge Institute Biological Resources Unit. Nine healthy 9-week-old
female C57BL/6 albino mice were imaged using the MSOT InVision 256TF (iThera
Medical GmbH, Munich, Germany) according to a previously established standard
operating procedure [45]. Additionally, six 28-week-old healthy female BALB/c
nude mice were imaged while inhaling 100% CO2 as their terminal procedure. In
both cases, imaging was performed at 10 wavelengths equally spaced between 700
and 900 nm averaging over 10 scans each. The mouse body, kidneys, spleen,
spine, and aorta were manually segmented by the authors using MITK.

###  2.6 Performance evaluation

The performance of the LSTM-based method was evaluated using the median
absolute error ($\epsilon\text{sO}_{2}$) between the estimate
($\hat{\text{sO}_{2}}$) and the ground truth/reference sO2:

$\epsilon\text{sO}_{2}=\text{median}(|\text{sO}_{2}-\hat{\text{sO}_{2}}|)$ (1)

Ground truth values are available for the in silico and in gello datasets. For
the in vivo measurements, reference values were based on literature. We
assumed sO2 of mixed murine blood to be 60% - 70% [41] and of arterial murine
blood sO2 under anaesthesia to be 94% - 98% [47]. For the CO2 terminal
procedure, we assumed that CO2 binds to haemoglobin, forming
carbaminohaemoglobin, which leads to oxygen unloading [48] and has an
absorption spectrum similar to deoxyhaemoglobin [49, 50], thus continuously
decreasing the actual [51] and measured global blood sO2. In humans, arterial
blood sO2 of 95%-100% was assumed [52].

###  2.7 Simulation gap measure

To predict the best-fitting training data set for a target application, one
could use the sO2 estimates of a trained algorithm to calculate error metrics,
such as the absolute estimation error, but this is only possible when ground
truth or reference sO2 values for a representative dataset are available. For
in vivo applications this is typically not the case and additionally,
unsupervised methods for performance prediction in the context of PA oximetry
remain largely unexplored.

Our alternative solution to this problem uses the Jensen-Shannon divergence
[53] (D${}_{\text{JS}}$). To this end, we compute D${}_{\text{JS}}$ between
the training data distributions and the target data distribution and calculate
the Pearson correlation coefficient between the resulting D${}_{\text{JS}}$
and $\epsilon\text{sO}_{2}$ of the LSTM estimates. D${}_{\text{JS}}$ measures
the distance between unpaired samples drawn from two probability distributions
$P$ and $Q$. D${}_{\text{JS}}$ is a symmetric version of the Kullback-Leiber
divergence [54] (D${}_{\text{KL}}$) and is defined as

$\text{D}_{\text{JS}}(P||Q)=\dfrac{1}{2}\text{D}_{\text{KL}}(P||M)+\dfrac{1}{2}\text{D}_{\text{KL}}(Q||M),$
(2)

where $M=\dfrac{1}{2}(P+Q)$ and the discrete D${}_{\text{KL}}$ is defined as
the relative entropy between two probability distributions:

$\text{D}_{\text{KL}}(P||Q)=\sum\limits_{x\in
X}P(x)\log\left(\dfrac{P(x)}{Q(x)}\right).$ (3)

To apply these measures, it is important to consider (1) handling the
multidimensional probability distributions arising from multi-wavelength
measurements and (2) the transformation of two sample distributions into the
same sample space. We calculate an aggregate $\overline{D_{\text{JS}}}$ by
calculating the mean over the distance for each wavelength in the spectrum:

$\overline{\text{D}_{\text{JS}}}=\dfrac{1}{N_{\lambda}}\sum\limits_{\lambda\in\Lambda}\text{D}_{\text{JS}}(P_{\lambda}||Q_{\lambda}),$
(4)

where $N_{\lambda}$ is the number of all available wavelengths $\Lambda$. To
standardise the sample space, a z-score normalisation is performed for each
spectrum and a histogram with 100 bins ranging from -3$\sigma$ to 3$\sigma$ is
created. A Python implementation of the Jensen-Shannon distance, available in
the Scipy (v1.10.1) package [55], was used.

##  3 Results

###  3.1 The LSTM-based method achieves accurate results across a range of
available input wavelengths.

The accuracy of the LSTM-based network architecture was first tested on the
BASE data set, when varying the numbers of available wavelengths
($N_{\lambda}$) for training or inference. $\epsilon\text{sO}_{2}$ was
extracted when training and testing the network on a certain fixed
$N_{\lambda}$, ranging from 3 to 41 wavelengths (Fig 2 A). We found that
$\epsilon\text{sO}_{2}$ increased as more wavelengths were used for training.
Nevertheless, $\epsilon\text{sO}_{2}$ for $N_{\lambda}=10$ was only slightly
higher than for $N_{\lambda}=41$ wavelengths. For $N_{\lambda}<10$,
$\epsilon\text{sO}_{2}$ starts to rapidly increase, which aligns with prior
literature [31] and is potentially exacerbated due to the ’vanishing gradient
problem’ [56] in LSTMs, which arises when a substantial portion of the input
parameter space consists of zeroes.

Next, a network trained at a fixed $N_{\lambda}$ (in this case
$N_{\lambda}=20$) was tested on data with a different $N_{\lambda}$ (Fig 2 B).
Accuracy was found to decrease rapidly if fewer wavelengths were used for
testing, but the error remains low if slightly more wavelengths are used.
Nevertheless, the results show that the LSTM-based network performs best if
the $N_{\lambda}$ used during inference matches the $N_{\lambda}$ during
training.

Figure 2: The LSTM-based method shows wavelength flexibility. A LSTMs were
trained with varying numbers of wavelengths ($N_{\lambda}$) to show that with
an increasing number of wavelengths, the accuracy of the predictions
increases. B LSTM trained at a given $N_{\lambda}$ can be applied to data with
different $N_{\lambda}$ but yields best results when $N_{\lambda}$ of the test
spectra matches that of the training spectra (indicated by the green vertical
line).

###  3.2 In silico cross-validation reveals the effect of changing simulation
parameters.

For each of the 24 simulation parameter variations of BASE, 500 data points
with random spatial vessel distributions and a fixed random number generator
seed for reproducibility were generated. To investigate the sensitivity of
data-driven $\text{sO}_{2}$ estimates to changes in training dataset
parameters, we first trained LSTM-based networks using all 41 available
wavelengths on each simulated training dataset and one on a mixed dataset
(ALL). We then performed cross-validation on all datasets by applying every
trained network to each of the datasets and calculating the median estimation
error $\epsilon\text{sO}_{2}$ (Fig 3 A). The $\epsilon\text{sO}_{2}$ values
range from 0.5% to over 35%. As expected, the best performance occurs when
testing on the training set, however, it is important to note that
$\epsilon\text{sO}_{2}$ is not zero (instead ranging from 0.5% - 5%),
suggesting that the network has not overfitted the training data.

Figure 3: Dimensionality reduction and cross-validation reveal systematic
differences between training datasets. A An LSTM-based network trained on each
dataset is then applied to every other dataset and all $\epsilon\text{sO}_{2}$
(median absolute error in percentage points) can be visualised as a
performance matrix. Dataset names are shortened for visibility but are
detailed in Table 1. B Uniform Manifold Approximation and Projection (UMAP)
projections of the four representative example datasets onto an embedding of
all training data. C Mapping the ground truth sO2 onto the same projection
reveals a correlation along the first UMAP axis.

We calculated a Uniform Manifold Approximation and Projection (UMAP) [57]
embedding of 200,000 randomly chosen spectra from all datasets. With UMAP, we
visualised the location of the spectra of representative datasets on this
embedding (Fig 3 B). Examining this visualisation indicates that some changes
in the dataset parameters result in highly different spectra (e.g.: adding a
layer of water on top of the tissue), while others lead to only minor
variations (e.g.: changing the background oxygenation). Labelling the UMAP
embedding with the corresponding ground truth sO2 values reveals a correlation
from low to high oxygenation along the horizontal UMAP axis (Fig 3 C).

From the in silico cross-validation heatmap (Fig. 3 A), we can derive several
key observations concerning the design of simulated data:

  1. 1.

Variation in background sO2 has minimal effect with the used 1% blood volume
fraction used, however, this could become more significant at higher blood
volume fractions.

  2. 2.

Resolution matters: Performance improves with higher spatial resolution
simulations in the training data, suggesting fine details in the spectral data
are important for accurate sO2 estimation.

  3. 3.

Illumination matters: When changing from a Gaussian to a point source
illumination the error increased.

  4. 4.

Chromophore inclusion: When the test dataset contained melanin, but not the
training data, the estimation error increased by an average of 5.8 percentage
points. When designing a training dataset, all chromophores relevant to the
target application should thus be included.

  5. 5.

Acoustic modelling causes systematic changes: Using acoustic modelling and
image reconstruction introduced systematic spectral changes that increase
$\epsilon\text{sO}_{2}$. This can have detrimental effects on the estimation
accuracy and should be considered during training data simulation.

  6. 6.

Training on a combined dataset is better: Including random samples from all
training datasets yielded more accurate estimates for all test datasets in
silico. It should already be noted that this finding was not reproducible on
the experimental datasets, suggesting that the LSTM-based method was not able
to generalise better by training on a combined dataset.

###  3.3 The Jensen-Shannon divergence correlates with the estimation error
and can therefore be used to identify the best training data set.

Given the variance in performance introduced by the choice of training data,
it is desirable to automatically determine the best training dataset for a
given algorithm and target application. The Jensen-Shannon divergence
(D${}_{\text{JS}}$) allows one to quantify the distance between the data
distribution of each dataset and the target data. We calculated the
correlation between D${}_{\text{JS}}$ and the median absolute sO2 estimation
error ($\epsilon\text{sO}_{2}$). When applying all networks, each trained on a
distinct training dataset, to the BASE dataset, we found that
D${}_{\text{JS}}$ correlates strongly with $\epsilon\text{sO}_{2}$ (Pearson
correlation coefficient R=0.76). The RES_0.15 dataset, which is the dataset
simulated at the highest resolution, and not the BASE dataset, achieved the
best D${}_{\text{JS}}$ score. This may be because BASE dataset is a subset of
the RES_0.15 dataset, resulting in low $\epsilon\text{sO}_{2}$ (Fig 4 A).

Figure 4: The Jensen-Shannon divergence (D${}_{\text{JS}}$) can predict
estimation performance. A D${}_{\text{JS}}$ correlates with the median
absolute sO2 estimation error $\epsilon\text{sO}_{2}$ when applying all
networks, each trained on a distinct training dataset, to the baseline dataset
(BASE). B D${}_{\text{JS}}$ for the D2O flow phantom data, which shows a
similar correlation with $\epsilon\text{sO}_{2}$. C After removing two
outliers, D${}_{\text{JS}}$ shows the same degree of correlation with
$\epsilon\text{sO}_{2}$ for the H2O flow phantom data.

Extending the analysis to experimentally acquired in gello D2O flow phantom
data showed a similar correlation (R=0.74). The network trained on SMALL
achieved the best score with D${}_{\text{JS}}=0.35$ and $\epsilon
sO_{2}=4.2\%$ (Fig 4 B). The application of D${}_{\text{JS}}$ to the H2O flow
phantom experiment initially revealed no correlation (R=-0.1), however,
networks trained on ILLUM_POINT and MSOT_ACOUS_SKIN were outliers and after
removing these, the correlation was comparable to other datasets (R=0.72, Fig
4 C). The network trained on WATER_4cm achieved the best score with
D${}_{\text{JS}}=0.47$ and $\epsilon sO_{2}=12.6\%$ (Fig 4 B). The presence of
outliers emphasises the importance of expert oversight when applying summary
metrics such as D${}_{\text{JS}}$.

D${}_{\text{JS}}$ correlates with $\epsilon sO_{2}$ across multiple simulated
and experimental data sets, providing evidence that the Jensen-Shannon
divergence can predict algorithm performance. This is particularly relevant
for previously unseen datasets where the true sO2 is unknown. For each
training dataset, a D${}_{\text{JS}}$ value can be computed by drawing random
samples from both the training and unseen dataset as outlined in section 2.7.
An LSTM-based network pre-trained on the dataset with the lowest corresponding
D${}_{\text{JS}}$ would then be chosen for data analysis since lower
D${}_{\text{JS}}$ correlates with a lower $\epsilon sO_{2}$. The same strategy
could also be used to guide an optimisation process to tailor the simulation
parameters to create a new training dataset that matches the target
application.

###  3.4 In gello testing shows that the LSTM method outperforms learned
spectral decolouring.

Algorithm performance on the oxygenation flow phantom was compared with LU as
the de facto state of the art and with a previously proposed learned spectral
decolouring method [20]. We show three example PA intensity images at 700 nm
of the D2O flow phantom at three time points t=[0 mins, 44 mins, 70 mins] (Fig
5 A), annotated with reference oxygenation (sO${}_{2}^{\text{ref}}$)
calculated from $pO_{2}$ reference measurements using the Severinghaus
equation [58].

Comparing the estimates of all methods trained on the SMALL training dataset
by plotting the estimated sO2 over time (Fig 5 B) reveals that the LSTM-based
method is, on average, more than twice as accurate as the learned spectral
decolouring (LSD) method and four times as accurate compared to LU. We show
the LU estimates for the three example images (Fig 5 C), demonstrating the
restricted dynamic range of LU estimates from t=0 min to t=44 min, ranging
from 80% to 32%, compared to a ground truth of 99% to 1%. Example images for
LSD (Fig 5 D) and the LSTM-based method (Fig 5 E) are shown as well, where the
latter can recover the widest dynamic range, extending from 97% to 5%. For
both methods, we chose SMALL as the training data set, as it was assigned the
lowest D${}_{\text{JS}}$ score.

Figure 5: Estimation of flow phantom data highlights performance dependence on
training dataset. Three example images of the D2O flow phantom are shown at
different time points (0, 22, 44 min) displaying A the photoacoustic signal
intensity at 700 nm with a red contour marking the blood-carrying tube and B
the sO2 estimations from different methods. We visually compare the
performance of linear unmixing (LU, C), learned spectral decoloring (LSD, D)
and the lstm-based method (LSTM, E) by plotting the sO2 estimations over the
same image section and time points shown in A.

###  3.5 Static in vivo testing shows the applicability of the proposed method
to different use cases.

Figure 6: Jensen-Shannon divergence (D${}_{\text{JS}}$) proves valuable for in
vivo data. LU and LSTM applied to measurements of the human forearm (A-F) and
mice abdomens (G-L). Panels A and G show the photoacoustic signal at 800 nm
and panels B and H show the spread of D${}_{\text{JS}}$ estimates for the
training datasets. Panels C and I show boxplots of the highlighted regions of
interest over all N=7 subjects. The horizontal lines show expected sO2 values
for arterial blood (red) and mixed blood (blue). sO2 images are shown for
models trained on BASE (D, J) and the best (E, K) and worst (F, L) dataset as
predicted by D${}_{\text{JS}}$. On the bottom right of these images, the value
distribution is shown as a grey histogram with the mean values of the regions
of interest highlighted in their respective colour.

Having examined performance in the phantom system with known ground truth, we
next apply the data-driven methods to static photoacoustic measurements of
seven human forearms (Fig 6 A-F) and seven mouse abdomens (Fig 6 G-L). For
each, we: show an example PA image, calculate D${}_{\text{JS}}$ on the data
distribution, and compare the estimated blood oxygenation in a region of
interest (human forearm: radial artery; mouse: aorta and spine) with
literature references.

For the forearm data (Fig 6 A), the MSOT_ACOUS_SKIN dataset is objectively the
best fit and was assigned the second highest D${}_{\text{JS}}$ score, whereas
ILLUM_POINT dataset was the worst fit (Fig 6 B). Some estimated sO2 values
were close to the expected radial artery sO2 value of 95%-100%. Notably, the
network trained on the MSOT_ACOUS_SKIN dataset results in $sO_{2}\approx$ 90%,
while the network trained on ILLUM_POINT dataset produces $sO_{2}\approx$ 95%
(Fig 6 C).

The sO2 estimates of the network trained on the MSOT_ACOUS_SKIN dataset (Fig 6
D) seem to have three primary modes: high values $>$ 80% from the vessel
structures, values in the 60% to 80% range in the surrounding tissue, and low
values from 10% to 50% in the skin and deep tissue. The sO2 estimates of the
network trained on the ILLUM_POINT dataset (Fig 6 E), on the other hand, are
concentrated on high sO2 values in all superficial tissue and only seem to be
below 85% in the skin and in deep tissue. The network trained on the BASE
dataset (Fig 6 F) estimates low sO2 values throughout the entire tissue and
does not exceed 80%. The ILLUM_POINT dataset, while seemingly successful if
only considering values from the radial artery, was assigned the lowest
D${}_{\text{JS}}$ value. The estimates and marginal histograms show that many
estimates are mapped to $<$20% and $>$90%, explaining the good score in the
radial artery. This finding demonstrates a common limitation of data-driven
oximetry methods, where the estimated value distributions do not agree with
expectations based on human physiology. The combination of all datasets (ALL)
results in an extremely low sO2 estimate in the radial artery (median sO2
$<50\%$), which contradicts the in silico cross-validation results and
indicates overfitting of the method to the training datasets.

For mouse images, the aorta and the area around the spinal cord are examined
(Fig 6 G), assuming from literature a physiological arterial sO2 of 92% - 98%
and for the spinal cord, a mixed arterial and venous blood with sO2 of 60% -
70%. WATER_4cm is objectively the best matching dataset and was assigned the
highest score according to D${}_{\text{JS}}$ (Fig 6 H). All data-driven
methods significantly increase the sO2 estimate in the aorta, and lie within
the desired bounds in the spinal cord (Fig 6 I). The BASE (Fig 6 L) and
WATER_4cm (Fig 6 J) datasets estimate a broad distribution and yield a higher
sO2 estimate in the aorta and a larger spread between sO2 in the aorta and
spinal cord. The limitations of the ILLUM_POINT (Fig 6 K) dataset are even
more evident in the mouse data, where even more pixels are either assigned 0%
or 100%.

###  3.6 Dynamic in vivo testing demonstrates that LSTM method can reveal
physiological processes.

To provide a quantifiable decrease in the in vivo sO2 levels in mice that
could test the capability of the LSTM-based method, we imaged N=6 mice when
experiencing asphyxiation breathing 100% CO2 [59]. sO2 estimates were
extracted from the major visible organs in the scan (spleen, kidneys, spinal
cord, and aorta) at 3 minutes before and 10 minutes after CO2 asphyxiation.

Figure 7: Data-driven methods estimate an increased sO2 dynamic range during
CO2 delivery compared to linear unmixing. A single representative mouse is
shown here. Panels A-D show the photoacoustic image (A) and sO2 estimation
results (B-D) 3 minutes before asphyxiation and panels E-H show the
photoacoustic image (E) and sO2 estimation results (F-H) 10 minutes after
asphyxiation. We show sO2 estimates for linear unmixing (B, F), the baseline
dataset (BASE, C, G), and the best dataset as predicted by the Jensen-Shannon
divergence (WATER_4cm D, H). Panels A and E show the outlines of the full-
organ segmentations, and all other panels B-D, F-H show outlines of the
segmented regions used in Table 2.

CO2 asphyxiation (before, Fig 7 A-D; after, Fig 7 E-F) increases the PA signal
amplitude at 800 nm (Fig 7 A,E) in the superficial organs up to a depth of
approximately 3 mm, while the centre of the mouse shows a decrease in signal.
sO2 values in the examined organs before CO2 asphyxiation are generally
consistent between LU and data-driven unmixing methods, with the network
trained on WATER_4cm estimating slightly lower sO2 values (5 to 8 percentage
points lower) (Table 2). Notably, the direction of predicted effects on sO2
levels aligns well between LU and data-driven unmixing methods. Still, the
effect sizes are up to three times greater when utilizing data-driven
approaches, demonstrating a wider dynamic range. Intriguingly, in the case of
the aorta, LU predicts an increase in sO2 levels despite the expected global
decrease in sO2 levels due to CO2 exposure. This may be caused by the
aforementioned increase in absorption coefficient in the periphery leading to
an increase in spectral colouring in depth.

Dataset $\rightarrow$ | LU | BASE | WATER_4cm
---|---|---|---
Organ $\downarrow$ | Before [%] |  $\Delta sO_{2}$ [%] | Before [%] |  $\Delta sO_{2}$ [%] | Before [%] |  $\Delta sO_{2}$ [%]
Body | 46$\pm$14 | -2 (n.s.) | 44$\pm$20 | -7 (**) | 37$\pm$17 | -8 (**)
Spleen | 45$\pm$10 | -10 (**) | 40$\pm$19 | -24 (**) | 35$\pm$15 | -26 (**)
Kidney | 52$\pm$6 | -14 (**) | 54$\pm$17 | -29 (**) | 47$\pm$11 | -32 (**)
Spine | 55$\pm$6 | -14 (**) | 60$\pm$11 | -35 (**) | 52$\pm$11 | -34 (**)
Aorta | 67$\pm$6 | +5 (n.s.) | 76$\pm$5 | -7 (*) | 68$\pm$17 | -18 (*)
Table 2: sO2 decreases during CO2 asphyxiation. Reported are the mean $\pm$
standard deviation of sO2 measurement before CO2 asphyxiation (Before) and the
change in the mean ($\Delta sO_{2}$) after the procedure. Values are reported
at a maximum depth of 3 mm into the mouse body for: linear unmixing (LU), an
LSTM-based network trained on the baseline training dataset (BASE), and on the
best fitting dataset according to the Jensen-Shannon divergence (WATER_4cm).
p-values for $\Delta sO_{2}$ are calculated using a non-parametric Mann-
Whitney U test and indicated by (n.s. = not significant; * = p$<0.05$; ** =
p$<0.01$). The segmentation masks are adjusted in Fig 7 B-D, F-H to show the
region of interest considered when calculating the results.

##  4 Discussion

We present an LSTM-based method for estimating sO2 from multispectral PA
images. We demonstrate that it can yield superior inference results compared
to the previously proposed learned spectral decolouring method while at the
same time being usable in a flexible manner, making it a promising candidate
to replace LU as the de facto state of the art. We also show that the
performance of the trained networks is highly dependent on the training data
and that changes in simulation parameters can lead to drastically different
data distributions. We thus propose to use the Jensen-Shannon divergence
(D${}_{\text{JS}}$) to complement the LSTM-based method. D${}_{\text{JS}}$
correlates with the median absolute sO2 estimation error and can thus be used
to select the best-fitting training dataset or to optimise the training data
distribution to fit the target application.

We highlight how the interplay of the LSTM-based method with D${}_{\text{JS}}$
can be used on a diversity of in vivo human and mouse data acquired with
different scanners and demonstrate that the LSTM-based method can reveal
significant dependencies in sO2 changes that conventional LU would fail to
identify. The LSTM-based method further consistently outperforms LU, with
estimated sO2 values aligning better with ground truth measurements in gello
and literature references in vivo. LU also shows significant outliers in
regions where imaging artefacts are present, which are either not present or
less pronounced when using the LSTM-based method.

Our in silico cross-validation reveals that acoustic modelling and image
reconstruction introduces systematic spectral changes not explained by the
initial pressure spectrum alone. Thus, an accurate digital model of the
clinically used device is crucial during data simulation to ensure the best
algorithm performance. While the combined dataset showed promising results in
silico, these were not replicated in the experimental datasets, which may
suggest that the network is overfitting and able to differentiate between the
different simulated datasets.

D${}_{\text{JS}}$ appears to be a valuable measure for determining the optimal
training dataset for the LSTM-based method, as it correlates with the median
absolute sO2 estimation error $\epsilon\text{sO}_{2}$ in silico (R=0.76) and
in gello (R=0.74/0.72). D${}_{\text{JS}}$ predicts a plausible training
dataset for all three in vivo applications tested in this study, where the
predicted value range was 0.4 - 0.7 on the mouse data, 0.5 - 0.8 on the
forearm data, and 0.3 - 0.7 on the CO2 data. With the development of fast and
auto-differentiable simulation pipelines [60], it should be possible to
optimise the simulation parameters for accurate sO2 estimates by iteratively
minimising D${}_{\text{JS}}$. When using differentiable implementations of
distribution distance measures, it might even be possible to integrate this
optimisation into an unsupervised training routine.

The in gello experiments with H2O as the coupling medium had high
$\epsilon\text{sO}_{2}$ errors for all sO2 estimation methods and
D${}_{\text{JS}}$ was consistently high. The predicted best training dataset
was WATER_4cm and the worst ILLUM_POINT, which is consistent with the in vivo
mouse experiments also having H2O as the coupling medium. Contrary to
D${}_{\text{JS}}$ prediction, ILLUM_POINT has the lowest
$\epsilon\text{sO}_{2}$, resulting in no correlation (R=-0.11); after removing
outliers, the correlation was on par with the other experiments (R=0.72). In
both cases, the estimation error was lower than predicted by
D${}_{\text{JS}}$; while the distributions were different, the trained
networks still managed to estimate accurate sO2 values from the training data.
Specifically for the ILLUM_POINT dataset, judging from the quasi-bimodal
distributions of the network’s estimates on in vivo data, appears to have
achieved a good agreement with high sO2 values purely by chance. Such outliers
are naturally obscured by summary measures such as D${}_{\text{JS}}$, thus
expert oversight is recommended.

The in vivo experiments with CO2 asphyxiation showed an increase in signal
amplitude at 800 nm in the periphery with a decrease in signal in the centre
of the mouse. These phenomena suggest an overall increase in the absorption
coefficient, which could be caused by a range of factors, including: blood
coagulation [61], erythrocyte aggregation [62], the presence of bicarbonate
ions (HCO${}_{3}^{-}$) in the blood formed by the dissociation of carbonic
acid into bicarbonate and hydrogen ions [48], or an increase in blood volume
due to blood stasis. Vasoconstriction of the capillaries leading to pallor
mortis could also play a role in the better visibility of the superficial
organs [63].

Having applied the LSTM-based method combined with D${}_{\text{JS}}$ to a
diverse range of data, we believe they provide a promising route to replacing
LU for pixel-based PA oximetry. Combining the flexibility in application of LU
with the increase in accuracy of deep learning-based unmixing methods is
attractive. The CO2 experiment suggests that LU can underestimate the effect
size of sO2 changes due to its compressed dynamic range and susceptibility to
artefacts. We thus recommend that LU should be complemented by a deep
learning-based estimation method. The codes and data of this study are
available open-source and the method will be integrated into the PATATO
toolkit [64], facilitating its widespread testing and future application.

Nonetheless, there remain limitations to this study that should be the subject
of further research. In this work, we investigated D${}_{\text{JS}}$
predictions on a predefined set of datasets. Based on the results, we believe
that D${}_{\text{JS}}$ might be suitable to be used within a minimisation
scheme to either manually or automatically determine the best choice of
simulation parameters for a given data set, but this remains to be
investigated. To account for spectral colouring artefacts more robustly, the
full 2D – or better 3D – tissue context should be taken into account within
the neural network, as quantitative photoacoustic imaging is only possible
with the full 3D context [12, 21]. Additionally, we have shown that good
training data are key for deep learning-based methods for PA oximetry, as
such, simulating data as realistically as possible is important. One direction
towards this is to make use of domain adaptation methods [28, 27, 65] that
adapt simulated training data to appear more realistic.

##  5 Conclusion

The presented LSTM-based approach for sO2 estimation from multispectral PA
images surpasses the performance of LU and a previously reported data-driven
oximetry method, making it a promising candidate to replace LU as the state of
the art. We address the impact of training data variations by introducing the
Jensen-Shannon divergence (D${}_{\text{JS}}$) as a valuable complement,
enabling selection of optimal datasets and fine-tuning for specific
applications. Our LSTM-based method consistently outperforms LU, aligning well
with ground truth measurements and literature references, while mitigating
outliers in regions prone to imaging artifacts. The combination of the
flexibility of the novel LSTM-based method with D${}_{\text{JS}}$ for training
data optimisation is a promising direction to make data-driven oximetry
methods robustly applicable for clinical use cases.

### Disclosures

The authors declare no conflict of interest regarding this work.

### Code, Data, and Materials Availability

The data and code to reproduce the findings of this study are openly
available.

The data is available under the CC-BY 4.0 license at:
https://doi.org/10.17863/CAM.105987

The code is available under the MIT license at:
https://github.com/BohndiekLab/LearnedSpectralUnmixing

### Acknowledgments

This work was funded by: Deutsche Forschungsgemeinschaft [DFG, German Research
Foundation] (JMG; GR 5824/1); Cancer Research UK (SB, TRE; C9545/A29580);
Cancer Research UK RadNet Cambridge (EVB; C17918/A28870); Against Breast
Cancer (LH); the Engineering and Physical Sciences Research Council (SB,
EP/R003599/1) The work is supported by the NVIDIA Academic Hardware Grant
Program and utilised two Quadro RTX 8000. The authors would like to thank Dr
Mariam-Eleni Oraiopoulou for the helpful discussions.

This work was supported by the International Alliance for Cancer Early
Detection, a partnership between Cancer Research UK (C14478/A27855), Canary
Center at Stanford University, the University of Cambridge, OHSU Knight Cancer
Institute, University College London and the University of Manchester.

This research was supported by the NIHR Cambridge Biomedical Research Centre
(NIHR203312). The views expressed are those of the authors and not necessarily
those of the NIHR or the Department of Health and Social Care.

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## References

  * [1] A. Fawzy, T. D. Wu, K. Wang, et al., “Racial and ethnic discrepancy in pulse oximetry and delayed identification of treatment eligibility among patients with covid-19,” JAMA internal medicine 182(7), 730–738 (2022).
  * [2] P. Beard, “Biomedical photoacoustic imaging,” Interface focus 1(4), 602–631 (2011).
  * [3] M. Li, Y. Tang, and J. Yao, “Photoacoustic tomography of blood oxygenation: a mini review,” Photoacoustics 10, 65–73 (2018).
  * [4] F. Knieling, C. Neufert, A. Hartmann, et al., “Multispectral optoacoustic tomography for assessment of crohn’s disease activity,” New England Journal of Medicine 376(13), 1292–1294 (2017).
  * [5] A. Karlas, N.-A. Fasoula, K. Paul-Yuan, et al., “Cardiovascular optoacoustics: From mice to men–a review,” Photoacoustics 14, 19–30 (2019).
  * [6] M. W. Dewhirst, Y. Cao, and B. Moeller, “Cycling hypoxia and free radicals regulate angiogenesis and radiotherapy response,” Nature Reviews Cancer 8(6), 425–437 (2008).
  * [7] D. Hanahan, “Hallmarks of cancer: new dimensions,” Cancer discovery 12(1), 31–46 (2022).
  * [8] J. Laufer, D. Delpy, C. Elwell, et al., “Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration,” Physics in Medicine & Biology 52(1), 141 (2006).
  * [9] C. Bench and B. Cox, “Quantitative photoacoustic estimates of intervascular blood oxygenation differences using linear unmixing,” in Journal of Physics: Conference Series, 1761(1), 012001, IOP Publishing (2021).
  * [10] S. Tzoumas and V. Ntziachristos, “Spectral unmixing techniques for optoacoustic imaging of tissue pathophysiology,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375(2107), 20170262 (2017).
  * [11] R. Hochuli, L. An, P. C. Beard, et al., “Estimating blood oxygenation from photoacoustic images: can a simple linear spectroscopic inversion ever work?,” Journal of biomedical optics 24(12), 121914–121914 (2019).
  * [12] B. Cox, J. G. Laufer, S. R. Arridge, et al., “Quantitative spectroscopic photoacoustic imaging: a review,” Journal of biomedical optics 17(6), 061202–061202 (2012).
  * [13] M. K. A. Singh, M. Jaeger, M. Frenz, et al., “In vivo demonstration of reflection artifact reduction in photoacoustic imaging using synthetic aperture photoacoustic-guided focused ultrasound (pafusion),” Biomedical optics express 7(8), 2955–2972 (2016).
  * [14] J. Jose, R. G. Willemink, W. Steenbergen, et al., “Speed-of-sound compensated photoacoustic tomography for accurate imaging,” Medical physics 39(12), 7262–7271 (2012).
  * [15] Y. Mantri and J. V. Jokerst, “Impact of skin tone on photoacoustic oximetry and tools to minimize bias,” Biomedical Optics Express 13(2), 875–887 (2022).
  * [16] T. R. Else, L. Hacker, J. Gröhl, et al., “The effects of skin tone on photoacoustic imaging and oximetry,” bioRxiv , 2023–08 (2023).
  * [17] G. P. Luke, K. Hoffer-Hawlik, A. C. Van Namen, et al., “O-net: a convolutional neural network for quantitative photoacoustic image segmentation and oximetry,” arXiv preprint arXiv:1911.01935 (2019).
  * [18] S. Agrawal, P. Gaddale, S. P. K. Karri, et al., “Learning optical scattering through symmetrical orthogonality enforced independent components for unmixing deep tissue photoacoustic signals,” IEEE Sensors Letters 5(5), 1–4 (2021).
  * [19] V. Grasso, R. Willumeit-Römer, and J. Jose, “Superpixel spectral unmixing framework for the volumetric assessment of tissue chromophores: A photoacoustic data-driven approach,” Photoacoustics 26, 100367 (2022).
  * [20] J. Gröhl, T. Kirchner, T. J. Adler, et al., “Learned spectral decoloring enables photoacoustic oximetry,” Scientific reports 11(1), 6565 (2021).
  * [21] J. Gröhl, M. Schellenberg, K. Dreher, et al., “Deep learning for biomedical photoacoustic imaging: A review,” Photoacoustics 22, 100241 (2021).
  * [22] H. Assi, R. Cao, M. Castelino, et al., “A review of a strategic roadmapping exercise to advance clinical translation of photoacoustic imaging: From current barriers to future adoption,” Photoacoustics , 100539 (2023).
  * [23] V. Sandfort, K. Yan, P. J. Pickhardt, et al., “Data augmentation using generative adversarial networks (cyclegan) to improve generalizability in ct segmentation tasks,” Scientific reports 9(1), 16884 (2019).
  * [24] L. Maier-Hein, M. Eisenmann, D. Sarikaya, et al., “Surgical data science–from concepts toward clinical translation,” Medical image analysis 76, 102306 (2022).
  * [25] L. Hacker, H. Wabnitz, A. Pifferi, et al., “Criteria for the design of tissue-mimicking phantoms for the standardization of biophotonic instrumentation,” Nature Biomedical Engineering 6(5), 541–558 (2022).
  * [26] J. Gröhl, T. R. Else, L. Hacker, et al., “Moving beyond simulation: data-driven quantitative photoacoustic imaging using tissue-mimicking phantoms,” arXiv preprint arXiv:2306.06748 (2023).
  * [27] A. K. Susmelj, B. Lafci, F. Ozdemir, et al., “Signal domain adaptation network for limited-view optoacoustic tomography,” Medical Image Analysis (2023).
  * [28] K. K. Dreher, L. Ayala, M. Schellenberg, et al., “Unsupervised domain transfer with conditional invertible neural networks,” arXiv preprint arXiv:2303.10191 (2023).
  * [29] C. Bench, A. Hauptmann, and B. Cox, “Toward accurate quantitative photoacoustic imaging: learning vascular blood oxygen saturation in three dimensions,” Journal of Biomedical Optics 25(8), 085003–085003 (2020).
  * [30] B. T. Cox, S. R. Arridge, and P. C. Beard, “Estimating chromophore distributions from multiwavelength photoacoustic images,” JOSA A 26(2), 443–455 (2009).
  * [31] S. Tzoumas, A. Nunes, I. Olefir, et al., “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nature communications 7(1), 12121 (2016).
  * [32] T. Kirchner and M. Frenz, “Multiple illumination learned spectral decoloring for quantitative optoacoustic oximetry imaging,” Journal of biomedical optics 26(8), 085001–085001 (2021).
  * [33] S. Manohar, “” in gello” imaging,” (2023).
  * [34] J. Gröhl, K. K. Dreher, M. Schellenberg, et al., “Simpa: an open-source toolkit for simulation and image processing for photonics and acoustics,” Journal of biomedical optics 27(8), 083010–083010 (2022).
  * [35] Q. Fang and D. A. Boas, “Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units,” Optics express 17(22), 20178–20190 (2009).
  * [36] B. E. Treeby and B. T. Cox, “k-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields,” Journal of biomedical optics 15(2), 021314–021314 (2010).
  * [37] C. Cai, K. Deng, C. Ma, et al., “End-to-end deep neural network for optical inversion in quantitative photoacoustic imaging,” Optics letters 43(12), 2752–2755 (2018).
  * [38] T. Chen, T. Lu, S. Song, et al., “A deep learning method based on u-net for quantitative photoacoustic imaging,” in Photons Plus Ultrasound: Imaging and Sensing 2020, 11240, 216–223, SPIE (2020).
  * [39] K. Hoffer-Hawlik and G. P. Luke, “Abso2luteu-net: tissue oxygenation calculation using photoacoustic imaging and convolutional neural networks,” Thesis (Bachelor’s) (2019).
  * [40] I. Olefir, S. Tzoumas, C. Restivo, et al., “Deep learning-based spectral unmixing for optoacoustic imaging of tissue oxygen saturation,” IEEE transactions on medical imaging 39(11), 3643–3654 (2020).
  * [41] D. G. Lyons, A. Parpaleix, M. Roche, et al., “Mapping oxygen concentration in the awake mouse brain,” Elife 5, e12024 (2016).
  * [42] M. Gehrung, S. E. Bohndiek, and J. Brunker, “Development of a blood oxygenation phantom for photoacoustic tomography combined with online po2 detection and flow spectrometry,” Journal of biomedical optics 24(12), 121908 (2019).
  * [43] A. Halevy, P. Norvig, and F. Pereira, “The unreasonable effectiveness of data,” IEEE intelligent systems 24(2), 8–12 (2009).
  * [44] F. Chollet et al., “Keras.” https://keras.io (2015).
  * [45] J. Joseph, M. R. Tomaszewski, I. Quiros-Gonzalez, et al., “Evaluation of precision in optoacoustic tomography for preclinical imaging in living subjects,” Journal of Nuclear Medicine 58(5), 807–814 (2017).
  * [46] I. Wolf, M. Vetter, I. Wegner, et al., “The medical imaging interaction toolkit,” Medical image analysis 9(6), 594–604 (2005).
  * [47] A. M. Loeven, C. N. Receno, C. M. Cunningham, et al., “Arterial blood sampling in male cd-1 and c57bl/6j mice with 1% isoflurane is similar to awake mice,” Journal of Applied Physiology 125(6), 1749–1759 (2018).
  * [48] H. Abramczyk, J. M. Surmacki, M. Kopeć, et al., “Hemoglobin and cytochrome c. reinterpreting the origins of oxygenation and oxidation in erythrocytes and in vivo cancer lung cells,” Scientific Reports 13(1), 14731 (2023).
  * [49] E. Dervieux, Q. Bodinier, W. Uhring, et al., “Measuring hemoglobin spectra: searching for carbamino-hemoglobin,” Journal of Biomedical Optics 25(10), 105001–105001 (2020).
  * [50] M. Taylor-Williams, G. Spicer, G. Bale, et al., “Noninvasive hemoglobin sensing and imaging: optical tools for disease diagnosis,” Journal of Biomedical Optics 27(8), 080901 (2022).
  * [51] P. R. Huber, M. E. Leimbach, W. L. Lewis, et al., “Co2 angiography,” Catheterization and cardiovascular interventions 55(3), 398–403 (2002).
  * [52] A. J. Williams, “Assessing and interpreting arterial blood gases and acid-base balance,” Bmj 317(7167), 1213–1216 (1998).
  * [53] J. Lin, “Divergence measures based on the shannon entropy,” IEEE Transactions on Information theory 37(1), 145–151 (1991).
  * [54] S. Kullback and R. A. Leibler, “On information and sufficiency,” The annals of mathematical statistics 22(1), 79–86 (1951).
  * [55] P. Virtanen, R. Gommers, T. E. Oliphant, et al., “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods 17, 261–272 (2020).
  * [56] R. M. Schmidt, “Recurrent neural networks (rnns): A gentle introduction and overview,” arXiv preprint arXiv:1912.05911 (2019).
  * [57] L. McInnes, J. Healy, and J. Melville, “Umap: Uniform manifold approximation and projection for dimension reduction,” arXiv preprint arXiv:1802.03426 (2018).
  * [58] J. W. Severinghaus, “Simple, accurate equations for human blood o2 dissociation computations,” Journal of Applied Physiology 46(3), 599–602 (1979).
  * [59] M. N. Fadhel, E. Hysi, H. Assi, et al., “Fluence-matching technique using photoacoustic radiofrequency spectra for improving estimates of oxygen saturation,” Photoacoustics 19, 100182 (2020).
  * [60] T. Rix, K. K. Dreher, J.-H. Nölke, et al., “Efficient photoacoustic image synthesis with deep learning,” Sensors 23(16), 7085 (2023).
  * [61] R. Su, S. A. Ermiliov, A. V. Liopo, et al., “Optoacoustic 3d visualization of changes in physiological properties of mouse tissues from live to postmortem,” in Photons Plus Ultrasound: Imaging and Sensing 2012, 8223, 105–111, SPIE (2012).
  * [62] R. K. Saha and M. C. Kolios, “A simulation study on photoacoustic signals from red blood cells,” The Journal of the Acoustical Society of America 129(5), 2935–2943 (2011).
  * [63] A. T. Schäfer, “Colour measurements of pallor mortis,” International journal of legal medicine 113, 81–83 (2000).
  * [64] T. Else, J. Gröhl, L. Hacker, et al., “Patato: a python photoacoustic tomography analysis toolkit,” in Photons Plus Ultrasound: Imaging and Sensing 2023, PC123790T, SPIE (2023).
  * [65] J. Li, C. Wang, T. Chen, et al., “Deep learning-based quantitative optoacoustic tomography of deep tissues in the absence of labeled experimental data,” Optica 9(1), 32–41 (2022).

Janek Gröhl completed his M.Sc. degree in medical computer science in 2016 and
completed his PhD in April 2021 under Prof. Lena Maier-Hein at the German
Cancer Research Center. Janek is working as a postdoctoral fellow with Prof.
Sarah Bohndiek funded by the Walter Benjamin Programme of the German Research
Foundation focusing on using machine learning methods and physical modelling
to tackle the problem of quantitative photoacoustic imaging.

Kylie Yeung completed her MRes in Connected Electronic and Photonic Systems in
2022, jointly between University College London and the University of
Cambridge, during which she investigated the effect of different simulated
training datasets for quantitative photoacoustic imaging. She is currently
pursuing a PhD at the University of Oxford.

Kevin Gu completed his MSci degree in experimental and theoretical physics in
2022 at the University of Cambridge. Kevin’s dissertation was focussed on a
wavelength-flexible approach to data-driven photoacoustic oximetry. He now
works as a quantitative researcher in the Netherlands.

Thomas R. Else completed his MSci degree in experimental and theoretical
physics in 2019 and is currently pursuing a PhD in medical sciences, both at
the University of Cambridge. His research focuses on clinical translation of
photoacoustic imaging, with a focus on equitable application of the technology
through open access software and evaluation of skin tone biases.

Monika Golinska completed her PhD in oncology at the University of Cambridge
in 2012. During her postdoc at Prof Sarah Bohndiek’s lab, she studied the
relationship between sex hormones and tumour microenvironment in breast
cancer. Monika is currently a MCSA fellow at the Medical University of Lodz
and researches endometriosis and its link with ovarian cancer.

Ellie V. Bunce completed her MPharmacol degree (Hons) at the University of
Bath in 2021. She is currently pursuing a PhD in medical sciences at the
University of Cambridge, investigating the potential of vascular targeted
therapies to enhance cancer radiotherapy response using novel imaging
approaches.

Lina Hacker is a Junior Research Fellow at the Department of Oncology at the
University of Oxford, UK. Her research is focused on the medical and technical
validation of novel approaches for cancer imaging, specifically relating to
tumour hypoxia. She received her PhD degree in Medical Sciences at the
University of Cambridge, UK, and holds a Master’s and Bachelor’s degree in
Biomedical Engineering and Molecular Medicine, respectively.

Sarah E. Bohndiek completed her PhD in radiation physics at the University
College London in 2008 and then worked in both the United Kingdom (at
Cambridge) and the United States (at Stanford) as a postdoctoral fellow in
molecular imaging. Since 2013, she has been a Group Leader at the University
of Cambridge, where she is jointly appointed in the Department of Physics and
the Cancer Research UK Cambridge Institute. She was appointed as full
professor of biomedical physics in 2020. Sarah was recently awarded the CRUK
Future Leaders in Cancer Research Prize and SPIE Early Career Achievement
Award in recognition of her innovation in biomedical optics.