2.1.

Computational Models

To assess the effects of skin tone on PAI in a controlled setting, we developed a computational forearm model imitating the structure and physical properties of the skin. The model was constructed in Python using absorption and reduced scattering spectra obtained from the literature.³⁶ Three variations of the model were used: a layer-only model for calibration with adding-doubling reflectance simulations, a varying blood oxygenation model for optical modeling, and a realistic model for full photoacoustic simulation. These distinctions were made to allow thorough characterization of the models under computational constraints. To bridge the gap between simulation and experiment, we also simulated optical absorption in the agarose phantom model (described in Sec. 2.5).

The first model consisted of three parallel layers. The first layer represented the epidermis, modelled as $60\mu \mathrm{m}$ thick, with optical absorption from a varying concentration of melanosomes. The melanosome volume fraction was chosen by first simulating diffuse reflectance (Sec. 2.2) over a range previously identified in the literature (1.3% to 43%³⁷) and adjusting the simulation range to align as closely as possible with individual typology angle (ITA) values associated with different Fitzpatrick types.³⁸ The resulting parameters were then qualitatively verified using a publicly available integrating sphere forearm skin reflectance dataset (Fig. S1 of the Supplementary Material).³⁹ We settled on six logarithmically-spaced values of melanosome concentration between 2% v/v and 40% v/v, representing the six Fitzpatrick skin types. The second layer, representing the dermis, was 1 mm thick, containing blood (0.2% v/v, 70% oxygenation), water (65% v/v) and lipid (20% v/v).³⁶ The third layer, representing generic tissue, filled the rest of the simulation volume. It consisted of blood (2.5% v/v, 70% oxygenation), water (65% v/v), and lipid (20% v/v). The absorption spectra for these three layers are shown in Fig. 1(a). Scattering properties were defined by a combination of Rayleigh and Mie scattering, assuming a constant scattering anisotropy of 0.9. These values are representative of typical values in each region.³⁶ This simplified layer model was used for validation purposes, using adding-doubling simulations of diffuse reflectance. The structure formed the basis of all simulations, with blood vessels added for Monte-Carlo simulations.

Fig. 1

Diffuse reflectance simulations validate Monte–Carlo simulation parameters and enable the assignment of simulations to equivalent Fitzpatrick skin types. (a) Optical absorption spectra in the tissue model were assigned according to the desired melanosome concentration (2% to 40% v/v, gradations of brown) in the skin layer, dermis below (gray) and background absorber (blue). (b) Diffuse reflectance spectra were calculated using the adding-doubling method. (c) Reflectance values simulated at 685 nm cover the physiological range of skin reflectance. (d) Red, green, and blue images were simulated from the reflectance spectra and the six models were assigned to the six Fitzpatrick skin types (shown by inset label I to VI).

To model the effects of changing blood oxygenation in a blood vessel on the photoacoustic initial pressure, Monte-Carlo simulations were run with the varying blood oxygenation model, giving the optical absorption distribution. This model was based on the layer structure described above with the addition of a cylinder of radius 1.25 mm, at a depth of 2.5 mm below the epidermis, representing a blood vessel of varying oxygenation (between 0% and 100%).

To evaluate the effects of a more realistic tissue model, a third model, the realistic forearm model, was subjected to both optical and acoustic simulations, giving full photoacoustic data that was processed using normal reconstruction and analysis procedures. Two cylinders were added to the layer model, one representing an artery (100%-oxygenated blood, radius 1.25 mm at a depth of 2.5 mm below the epidermis) and one representing a vein (70% oxygenated blood, radius 0.625 mm at a depth of 2.5 mm).

We also ran optical simulations of the cylindrical phantoms (Sec. 2.5), with approximately identical optical absorption and scattering properties as the agarose phantoms. We defined the model in Python, with a grid size of $100\mu \mathrm{m}$ in an array of size ($220\times 220\times 220$). The tube in the center of the phantom was modeled as a cylindrical absorber of radius 0.75 mm containing blood. Concentric absorbers were added representing the inner part of the phantom (23 mm diameter, water absorption and intralipid scattering) and the outer layer of the phantom (1 mm thick, varying melanosome volume fraction, intralipid scatterer). The absorption coefficients of the melanin layer were defined by finding the melanosome volume fraction that gave the same optical absorption coefficient at 700 nm as measured in the double integrating sphere system (Table S1 in the Supplementary Material, Sec. 2.6).

2.2.

Adding-Doubling Simulations

The layer-only skin model was validated using the adding-doubling model of diffuse reflectance⁴⁰ to allow for fast evaluation of the reflectance spectra at several wavelengths (450 nm to 900 nm in 20 nm steps) in the visible and near-infrared ranges. We additionally calculated the reflectance at 685 nm to enable direct comparison with literature values.⁴¹ The reflectance spectra were converted into an RGB image and an ITA metric using a D65 standard light source and a 10-deg standard observer in the Python color-science library.⁴² ITA has been proposed as a quantitative measure of skin tone and can be obtained from low-cost colorimeters. ITA was calculated from skin reflectance spectra using the following equation:

2.3.

Monte-Carlo Simulations

We simulated the PAI optical forward model using the varying blood oxygenation model described above using a custom version of MCX ( https://github.com/IMSY-DKFZ/mcx), giving initial pressure distributions that correspond to an idealized photoacoustic image. The model was defined on a grid with a pixel size of $60\mu \mathrm{m}$, which allowed the simulation to run in a sufficiently short time while ensuring that the simulations were sufficiently detailed. The light transport through the model was simulated with ${10}^7$ photons per simulation, at 21 equally spaced wavelengths between 700 and 900 nm inclusive.⁴³

We used an illumination geometry, which approximates a clinical PAI system (MSOT Acuity, iThera Medical GmbH, Munich, Germany), consisting of a single optical fiber with a beam divergence of 8.66 deg, 43.2 mm above the surface of the probe, directed at an angle of 22.4 deg from the imaging plane, intersecting the imaging plane 2.8 mm outside of the probe.⁴⁴ Optical simulations of the agarose phantoms were run using the same customized version of MCX with an illumination geometry approximating a pre-clinical photoacoustic system (MSOT inVision 256-TF, iThera Medical GmbH, Munich Germany), defined on a grid of size $100\mu \mathrm{m}$. The illumination geometry consisted of five radially distributed fiber bundle pairs.⁴⁴ The version of MCX used in these experiments is available on GitHub: https://github.com/BohndiekLab/melanin-phantom-simulation-paper. ${10}^7$ photons were simulated per fiber bundle, for a total of ${10}^8$ photons. All simulations were run on a computer with 2 NVIDIA Quadro 8000 GPUs (48 GB GPU RAM each), Intel Xeon Gold 6230 CPU at 2.10 GHz (40 cores), and 256 GB RAM.

2.4.

Acoustic Simulations

To determine the impact of acoustic modeling, we ran the forward acoustic model in j-Wave⁴⁵ to calculate the pressure time series detected by the ultrasound transducer array for the initial pressure distributions. Simulations were run in three dimensions, with a grid size of $240\mu \mathrm{m}$, and a speed of sound of $1500\mathrm{m}{\mathrm{s}}^{-1}$. Three-dimensional initial pressure distributions were taken from Monte-Carlo simulations and down-sampled by a factor of four, and their propagation was simulated. The resulting pressure distributions were evaluated using band-limited interpolation sensors in the position of the sensors in the photoacoustic system. The pressure time series were recorded and processed in the same way as the experimental data.

2.5.

Phantom Design

Phantoms were used to validate the simulations in a real photoacoustic system. A custom-made 3D-printed polylactic acid mold was constructed to fabricate a two-layer cylindrical phantom with an inner scattering compartment and an outer compartment containing melanin to model the epidermis. The mold consisted of two cylindrical compartments, one of radius 9 mm and one of radius 10 mm. A polyvinyl chloride tube (inner diameter 1.5 mm, outer diameter 2.1 mm; VWR 228-3857) was placed in the center of the phantom mold to enable blood flow through the phantom.

Agarose-based phantoms were fabricated using standard methods.⁴⁶ The inner phantom mixture consisted of a mixture of agarose for mechanical robustness and intralipid for scattering. 1.5% (w/v) agarose (Sigma-Aldrich 05039) was added to Milli-Q water and heated until it dissolved. Once the mixture had cooled to $\sim 40°\mathrm{C}$, 2.08% (v/v) of pre-warmed intralipid (20% emulsion, Sigma-Aldrich I141) was added and mixed, before being poured into the inner compartment of the phantom mold and allowed to set, constituting the inner 9 mm radius of the phantom. Synthetic melanin (Sigma, M8631, $2.5\mathrm{mg}{\mathrm{mL}}^{-1}$ stock in dimethyl sulfoxide) was added to a separate batch of the base mixture at the desired concentrations (Table S1 in the Supplementary Material). Once the inner compartment was set, the outer layer mixture was added and again allowed to set at room temperature (Fig. S2 in the Supplementary Material).

2.6.

Phantom Characterization

A custom double integrating sphere system was used to characterize the optical absorption and reduced scattering coefficients of the phantom material. Samples were placed between two integrating spheres (Avantes, AvaSphere-50, 50 mm internal diameter), both of which were connected to a spectrometer (Avantes, Starline Avaspec-2048) via an optical fiber. The reflectance integrating sphere was connected to a broadband light source (Avantes, Avalight-HAL-s-mini) via an optical fiber. Circular samples (approximate diameter: 5 cm) of the agarose mixture with and without melanin were prepared with a thickness of $\sim 3\mathrm{mm}$. Sample thickness was measured using vernier callipers. Each agarose sample was placed between the two spheres and the transmittance and reflectance were measured. Optical scattering and absorption coefficients were calculated using the inverse adding doubling method,⁴⁷ assuming a scattering anisotropy factor of 0.89 and refractive index of 1.34.⁴⁸ These measurements were made on three replicates of each sample, over a wavelength range of 450 to 900 nm at 1 nm steps.

2.7.

Blood-Flow Circuit

A blood-flow circuit was used to control the oxygenation of blood flowing through the tube centered in the two-layer cylindrical phantom, enabling photoacoustic measurements of blood with known ground-truth blood oxygenation as described previously.^49,50 Briefly, we measured the oxygen partial pressure (${\mathrm{pO}}_2$) and the photoacoustic spectrum of the blood simultaneously. We used fresh whole blood from Wistar rats (6 to 9 months, Charles River Laboratories). All blood samples were anti-coagulated in ethylenediaminetetraacetic acid (EDTA), stored at 4°C, and processed within 48 h. The blood was pumped around the circuit using a peristaltic pump (Fisher Scientific CTP100), and the oxygenation was controlled by initially oxygenating the blood using hydrogen peroxide (30% (w/w) in deionized water, Sigma-Aldrich 7722-84-1) and then slowly injecting 0.03% w/v sodium dithionite (ACROS Organics 7775-14-6, in PBS) using a syringe pump ($10\mu \mathrm{L}{\min }^{-1}$, Harvard, MKCB2159V). To measure the blood ${\mathrm{pO}}_2$, two oxygen fluorescence quenching needle probes (Oxford Optronix, NX-BF/O/E) were placed before and after the phantom; the average of these measurements was taken as the ground-truth ${\mathrm{pO}}_2$ in the phantom. ${\mathrm{pO}}_2$ was converted to blood oxygenation (${\mathrm{sO}}_2$) using the Severinghaus equation. ⁵¹ The values were extracted using custom code on an Arduino UNO and fed to a computer over USB where they were read and stored using custom code in MATLAB. Experiments were conducted at 37°C.

2.8.

Animal Procedures

All animal procedures were conducted under personal and project licences issued under the United Kingdom Animals (Scientific Procedures) Act, 1986 and reviewed by the Animal Welfare and Ethical Review Board at the CRUK Cambridge Institute. The protocols were approved locally under compliance forms CFSB2317 and CFSB2348.

To understand the effects of skin pigmentation in a model in vivo setting, we compared non-pigmented albino mice ($n=10$, C57BL/6 albino, Janvier Labs, France, 10 to 13 weeks old) to pigmented “B6” mice ($n=10$, C57BL/6J, Charles River, 10 to 13 weeks old). The B6 mice were further subdivided qualitatively based on the visual appearance of pigmentation patches in color photographs (Fig. S3 in the Supplementary Material) into a non-pigmented group ($n=8$ scans) and a pigmented group ($n=11$ scans). Albino mice were imaged once and B6 mice were imaged twice (2 weeks apart) to allow the development of pigmentation after hair removal. Repeat measurements of the same mouse were treated as independent measurements for statistical analysis given the time lapse between measurements.

2.9.

Photoacoustic Imaging

Phantoms and mice were imaged using a pre-clinical multispectral photoacoustic tomography system (MSOT inVision-256TF, iThera Medical GmbH).⁴⁶ Briefly, excitation pulses were provided by a tunable optical parametric oscillator, pumped by a nanosecond Nd:YAG laser (10 Hz repetition rate up to 7 ns pulse duration). A custom optical fiber assembly creates a diffuse ring of uniform illumination over the imaging plane within the sample. The sample is coupled to the transducers using a water bath. For ultrasound detection, an array of 256 toroidally focused ultrasound transducers covering a 270 deg circular arc and radius of 4 cm was used with a center frequency of 5 MHz and 60% bandwidth.

For phantom imaging, a single cross-sectional slice of the phantom was imaged continuously during the blood deoxygenation procedure until the blood was completely deoxygenated. Phantoms were placed in a custom animal holder (iThera Medical GmbH, Munich, Germany), wrapped in a thin polyethylene membrane, with water used to couple the phantom to the membrane. Imaging was performed at wavelengths from 700 to 900 nm inclusive with 10 nm steps and frame averaging was enabled (10 averages) to improve the signal-to-noise ratio.

Mice were prepared for imaging according to our standard operating procedure.⁴⁶ Briefly, the mice were shaved and depilatory cream was used to remove excess hair. The mice were anaesthetized using $<3\%$ isoflurane in a 50% pure oxygen/medical air mixture and placed in a custom animal holder (iThera Medical GmbH, Munich, Germany), wrapped in a thin polyethylene membrane, with ultrasound gel (Aquasonic Clear, Parker Labs) used to couple the skin to the membrane. The holder was then placed within the PAI system and immersed in deionized water maintained at 36°C for Hb and ${\mathrm{HbO}}_2$ imaging acquisition. The animal respiratory rate was maintained in the range 70 to 80 b.p.m. with 1.5% isoflurane concentration for the entire scan. A single cross-sectional slice through the kidneys and spleen was imaged. Frame averaging was enabled (four averages) to minimize the effects of breathing motion artifacts. 10 repeat scans were taken at each of the following wavelengths 700, 730, 750, 760, 770, 800, 820, 840, 850, and 880 nm.

2.10.

Image Reconstruction and Processing

Simulated, phantom, and in vivo photoacoustic data were processed using the Python photoacoustic tomography analysis toolkit (PATATO, https://github.com/bohndieklab/patato).

Acoustic simulation images were reconstructed by filtered backprojection with a grid spacing of $60\mu \mathrm{m}$ and a field of view of $18\mathrm{mm}\times 40\mathrm{mm}$. Phantom and in vivo images were reconstructed by filtered backprojection with a grid spacing of $75\mu \mathrm{m}$ and a field of view of 25 mm. A band-pass filter was applied to all of the time series data, with a bandpass range of 5 to 7 MHz.

Linear spectral unmixing was applied to all simulated, phantom, and in vivo datasets using the matrix pseudo-inverse function in NumPy⁵² and literature values for the optical absorption spectra of oxy- and deoxy-hemoglobin.³⁶ Photoacoustic estimates of the blood oxygenation (denoted with a superscript EST ${\mathrm{sO}}_2^{\mathrm{EST}}$ to differentiate it from the true blood oxygenation, ${\mathrm{sO}}_2$) were made pixel-wise from the unmixing components as the ratio of the oxyhemoglobin to the total hemoglobin. To evaluate the potential for a machine-learning approach to improve photoacoustic oximetry in the presence of melanin, a learned spectral decoloring ³¹ algorithm was trained on a numerical tissue model containing varying levels of melanin in the epidermis layer (0.1% to 5% melanosome v/v). We simulated 500 unique volumes $19.2{\mathrm{mm}}^3$ in volume with a resolution of 0.3 mm that contained up to nine blood vessels with a random radius and blood oxygenation. The tissue background contained a blood concentration of 1% (v/v) and an oxygenation of 70%. These design parameters were developed independently of this study and simulations were obtained using the SIMPA framework.⁴⁴ From the simulated initial pressure distributions, we used a histogram-based gradient boosting regression tree as the inversion algorithm with a maximum depth of 16 and otherwise default hyperparameters as implemented in the scikit-learn (version 1.2.1) Python package.

Unless otherwise specified, manual polygon regions of interest (ROIs) were drawn and average spectra or unmixed components were calculated by masking the defined regions. In phantoms, the ROIs were taken around the blood tube and the melanin compartment. The melanin compartment was defined by drawing an ROI around the whole phantom and eroding the polygon by 1 mm. The difference mask between the drawn and eroded polygons was applied to the photoacoustic images, giving the multispectral data within the melanin region. Mean photoacoustic and unmixed component values were extracted. In mice, ROIs were drawn around the whole mouse body. The body ROI was eroded by 1 mm to subdivide it into a “skin” region (outer 1 mm) and a “body” region (inner region). Negative pixel values were included in the median spectra calculation to illustrate image reconstruction artifacts. For unmixed ${\mathrm{sO}}_2^{\mathrm{EST}}$ calculation, pixels were excluded if either of the unmixed oxy- or deoxy-hemoglobin coefficients were negative; the overall ${\mathrm{sO}}_2^{\mathrm{EST}}$ was calculated as the ratio of the median oxyhemoglobin component to the median total hemoglobin. The median was chosen because of the high negative pixel percentage in the body of the pigmented mice. Pixels near the open region of the detector were manually excluded from the analysis as the image reconstruction is less well posed in that region, leading to artifacts.

Line plots of the photoacoustic signal through the mice were generated by constructing a line through the kidneys and vena cava of the mice. To ensure constant placement of this line, it was constructed to be 6.5 mm from the surface of the mouse closest to the spline, and perpendicular to an axis through the spine and aorta. This construction was observed to pass through the kidneys and vena cava in albino and non-pigmented mice. The same construction was applied to the heavily pigmented mice, where the kidneys and spleen were not visible, to give a direct comparison. The values of the reconstructed image, total hemoglobin map, and blood oxygenation were calculated along this line using the interpnd function of SciPy.⁵²

2.11.

Statistics

Statistical analyses were conducted in Python using the statsmodels library. The relationship between continuous variables was quantified by linear regression. A log transformation was applied where specified to ensure that the residual normality assumption is upheld in each linear model. The residual normality assumption was tested using the Shapiro-Wilk test ($p>0.1$ for all models). For the in vivo data, the relationship between each response variable and the mouse pigmentation was assessed using a linear model, with both mouse strain (Albino or B6) and pigmentation (not pigmented or pigmented) as independent variables. The model related the response variable (${\mathrm{sO}}_2$ or photoacoustic signal) to a sum of the effect due to mouse strain and pigmentation with no interaction term.

3.1.

Computational Model of Skin Shows Good Agreement with Literature Based on Adding-Doubling Reflectance Measurements

We first validated our computational model of human skin, with the assigned optical absorption spectra in the epidermis [Fig. 1(a)], by calculating diffuse reflectance spectra using the adding-doubling method. We chose six values of melanosome concentration to represent the typical physiological range observed in the population. The calculated spectra show decreasing reflectance with increasing melanosome concentration and broadly increasing reflectance with wavelength, as one would expect given the profile of optical absorption of melanin [Figs. 1(a) and 1(b)]. We then calculated the ITA for each reflectance spectrum, allowing direct comparison with in vivo measurements from skin colorimeters. These ITA values (67 deg to $-47\mathrm{deg}$) are consistent with previous observations made in vivo using the Fitzpatrick skin type (Table 1).⁵³ Reflectance at 685 nm was evaluated for each tissue model, giving a range of reflection coefficients from 0.20 to 0.72 over our simulated range of melanosome concentrations [Fig. 1(c), Table 1] again consistent with previous literature.⁴¹ To enable qualitative evaluation of our skin model, the reflection coefficients were converted into red, green, and blue (RGB) images [Fig. 1(d), Table 1]. Our results align closely with literature values made of skin reflectance using a single integrating sphere (Fig. S1 in the Supplementary Material).³⁹ For the remainder of the paper, results from different Fitzpatrick skin types are color coded according to the calculated RGB values for consistency [Fig. 1(d)].

Table 1

Parameters of the skin simulations (melanosome concentration, absorption coefficient) and quantities derived from adding-doubling simulations (reflectance, ITA, and HEX color value). The equivalent Fitzpatrick skin type was assigned based on literature values.53

3.2.

Photoacoustic Simulations with the Skin Model Show Spectral Coloring Due to Skin Pigmentation

Multispectral Monte-Carlo simulations using the computational skin model with an embedded blood vessel showed how the optical absorption values are affected by melanosome concentration and hence Fitzpatrick skin type [Fig. 2(a)]. In simulations of fully oxygenated blood, we observed an increase in the optical absorption within the skin at all wavelengths with increasing Fitzpatrick type [Fig. 2(b)]. The optical absorption in the epidermis increased more than 10-fold between Fitzpatrick I and Fitzpatrick VI at 700 nm [379 a.u. in FP I, 4128 a.u. in FP VI, Fig. 2(c)]. Spectral coloring is observed within the blood vessel region, which manifests as a decrease in the mean-normalized absorption at lower wavelengths and a corresponding increase at higher wavelengths [Fig. 2(d)]. Absorption measurements within the blood vessel region, however, decrease by a factor of 1.7 from Fitzpatrick I to Fitzpatrick VI [44.6 a.u. in FP I, 26.7 a.u. in FP VI, Fig. 2(e)]. Combined optical and acoustic simulations of a realistic tissue model revealed additional image reconstruction artifacts associated with backprojection, including negative pixels and streaking artifacts (Fig. S4 in the Supplementary Material). Quantitative results from the blood and skin regions in the acoustic simulations were consistent with those made with solely optical simulations, so the remainder of the simulation study was performed using the Monte-Carlo optical simulation alone for computational practicality.

Fig. 2

Monte–Carlo simulations of light transport through the skin model show decreased signal and increased spectral coloring in an underlying blood vessel with increasing Fitzpatrick skin type. (a) Photoacoustic initial pressure (upper), unmixed total hemoglobin (THb, middle), and unmixed blood oxygenation (${\mathrm{sO}}_2^{\mathrm{EST}}$, lower) in representative examples for fully oxygenated blood in the vessel. (b) Epidermis PA initial pressure as a function of wavelength for all Fitzpatrick types and (c) at 700 nm as a function of Fitzpatrick type. (d) Normalized blood PA initial pressure as a function of wavelength for all Fitzpatrick types and (e) at 700 nm as a function of Fitzpatrick (FP) type. Fitzpatrick types are denoted with graded brown coloration according to the RGB definitions in Fig. 2(d). Scale bar = 2 mm applies to all images in panel (a).

Spectral coloring observed in the blood vessel biased the linear spectral unmixing blood oxygenation estimate [Fig. 3(a)]. Encouragingly, the bias can be somewhat overcome when using learned spectral unmixing trained on independently generated numerical tissue models containing varying melanin concentrations [Fig. 3(a)]. Considering the spectral changes observed [Figs. 2(d) and 2(e)], it appears that learned unmixing can better cope with the spectral corruption that presents as a reduction in normalized optical absorption at low wavelengths and a relative increase at higher wavelengths for higher Fitzpatrick skin types. With fully oxygenated blood, linear unmixing estimates of photoacoustic ${\mathrm{sO}}_2^{\mathrm{EST}}$ reveal a non-linear relationship with Fitzpatrick skin type [Fig. 3(b)] and a linear relationship with melanosome concentration [$p<0.001$, gradient = 0.218, 95% CI (0.200, 0.235)]. Blood oxygenation estimated with learned unmixing also increases linearly with melanosome concentration but with a smaller gradient [$p<0.001$, gradient = 0.150, 95% CI (0.117, 0.183)] [Fig. 3(b)].

Fig. 3

Estimation of blood oxygenation in simulated data using linear and learned spectral unmixing. (a) ${\mathrm{sO}}_2^{\mathrm{EST}}$ estimates increase with increasing Fitzpatrick type when using linear spectral unmixing, whereas estimates are closer to the ground truth when using learned unmixing, particularly at high and low blood oxygenation levels. (b) ${\mathrm{sO}}_2^{\mathrm{EST}}$ estimates for fully oxygenated blood using both approaches are shown as a function of Fitzpatrick type and melanosome volume fraction. The shaded area shows 95% confidence intervals.

3.3.

Blood Flow Phantoms Confirm Melanin-Associated Bias in an Experimental Photoacoustic System

Cylindrical phantoms were enclosed in a melanin-containing outer layer with concentrations equivalent to Fitzpatrick types below I, II, and IV [Fig. 4(a)]. A blood flow circuit was connected to a tube through the center of the phantom, allowing continuous monitoring of the blood oxygenation. The optical properties of the phantoms were confirmed using double integrating sphere measurements (Fig. S5 and Table S1 in the Supplementary Material). By imaging the phantoms over time, as the blood was adjusted from fully oxygenated to fully deoxygenated, the photoacoustic measurements were related to the known blood oxygenation (Fig. S6 in the Supplementary Material). Images of these phantoms show strong absorption in the melanin outer layer and blood flow tube, with image reconstruction artifacts visible around the blood inclusion and near the melanin layer, particularly visible with the higher concentrations of melanin [Fig. 4(b)]. The photoacoustic signal in the outer layer increases with melanin concentration [Figs. 4(b) and 4(c)]. Wavelength-dependent spectral coloring is observed in the blood spectrum with increasing melanin concentration, with normalized photoacoustic signals relatively lower at shorter wavelengths and higher at longer wavelengths, consistent with the optical absorption properties of melanin [Fig. 4(d)] and simulation findings.

Fig. 4

Skin-mimicking blood-flow phantoms show spectral coloring. (a) Schematic diagram of the agarose phantom setup, showing a 1 mm-thick outer layer containing synthetic melanin and a blood inclusion with varying oxygenation. (b) Reconstructed photoacoustic (PA) images at 800 nm (upper), linear unmixed total hemoglobin (THb, middle), and linear unmixed blood oxygenation (${\mathrm{sO}}_2^{\mathrm{EST}}$, lower) maps. Image reconstruction artifacts can be observed near the blood inclusion and near the outer melanin layer (particularly clear in ${\mathrm{sO}}_2^{\mathrm{EST}}$). (c) The photoacoustic signal in the outer cylinder of the phantom shows an increase with increasing melanin concentration. (d) Normalized photoacoustic spectra from the blood region show the effect of spectral coloring over a range of blood oxygenation levels. Scale bar = 5 mm.

Furthermore, examining the unmixed ${\mathrm{sO}}_2^{\mathrm{EST}}$ as a function of the reference blood ${\mathrm{sO}}_2$ measurements reveals the effect of spectral coloring on PAI signal quantification [Figs. 5(a) and 5(b)], with higher melanin concentrations giving higher photoacoustic estimates of blood oxygenation when using linear spectral unmixing. Learned unmixing again reduces the disparity between low and high melanin concentrations, particularly at higher blood oxygenation levels [Fig. 5(b)]. Linear unmixing estimates of ${\mathrm{sO}}_2^{\mathrm{EST}}$ increase with melanin concentration [$p<0.001$; $\text{gradient}=0.162\mathrm{mL}{\mathrm{mg}}^{-1}$, 95% CI ($0.161\mathrm{mL}{\mathrm{mg}}^{-1}$, $0.163\mathrm{mL}{\mathrm{mg}}^{-1}$)], whereas learned estimates do not [$p=0.50$; $\text{gradient}=0.051\mathrm{mL}{\mathrm{mg}}^{-1}$, 95% CI ($-0.603\mathrm{mL}{\mathrm{mg}}^{-1}$, $0.706\mathrm{mL}{\mathrm{mg}}^{-1}$)]. Analyzing simulations of digital twins of these phantoms in the tomographic PAI system showed consistent trends (Figs. S7 and S8 in the Supplementary Material).

Fig. 5

Comparison of linear unmixing and learned unmixing of blood oxygenation in skin-mimicking phantoms. Linear unmixing and learned unmixing blood oxygenation estimates (${\mathrm{sO}}_2^{\mathrm{EST}}$) are calculated and compared for varying melanin concentrations (a) Note that discontinuities in this plot arose from the rapid decrease in blood oxygenation from around 0.5 in the $0.21\mathrm{mg}{\mathrm{mL}}^{-1}$. This is also shown as a function of melanin concentration at a high blood oxygenation, 95% (b). The line of best fit from a linear model is shown. Linear unmixing ${\mathrm{sO}}_2^{\mathrm{EST}}$ increases significantly ($p<0.001$) with melanin concentration, whereas learned unmixing increases but not significantly ($p=0.50$).

3.4.

Mouse Models Underscore the Severity of the Impact of Skin Pigmentation on Photoacoustics in vivo

Two closely related mouse models showing differential pigmentation were used to observe the effects of skin pigmentation on PAI in vivo (Figs. S3 and S9 in the Supplementary Material). Qualitative inspection of the image data clearly shows the increase in photoacoustic signal from the skin region in the presence of pigmentation and an associated decrease in photoacoustic signals measured from the body region in those pigmented mice [Figs. 6(a) and 6(b)]. The spectra extracted from the ROIs show spectral coloring is clearly present in the pigmented mice. Spectra also demonstrate an increased presence of artifacts, such as negative pixels [Fig. 6(b)]. Line profiles across the mice [Figs. 6(c)–6(e)] further emphasize these findings, highlighting the high level of optical absorption in the skin of the pigmented mice, artifactual increase in the total hemoglobin signals from the same region, and an increase in noise in the body of the mouse in the ${\mathrm{sO}}_2^{\mathrm{EST}}$ image.

Fig. 6

Photoacoustic images of mouse spleen and kidneys in albino and pigmented mice show substantial variations in signal intensity within the skin region and body regions. (a) Reconstructed photoacoustic images at 700 nm (upper), linear unmixed total hemoglobin (THb, middle), and linear unmixed blood oxygenation (${\mathrm{sO}}_2^{\mathrm{EST}}$, lower) maps. Regions of interest (body = central region of the blue dashed line, skin = outer region of the blue dashed line) and line profiles used for quantitative analysis are shown on the single wavelength image. (b) Median photoacoustic spectra from the skin and body of albino mice (blue, $n=10$), non-pigmented B6 mice (yellow, $n=8$) and pigmented B6 mice (orange, $n=11$). Line profiles through the kidneys, taken 6.5 mm from the surface of the mouse nearest the spine, of the photoacoustic signal at 700 nm (c), THb (d), and ${\mathrm{sO}}_2^{\mathrm{EST}}$ (e). The shaded area shows the standard error. Scale bar = 5 mm.

Considering data from a single wavelength [700 nm, Fig. S9(A) in the Supplementary Material], the photoacoustic signal in the skin region of pigmented mice is significantly higher than in non-pigmented and albino mice ($p=0.010$). Interestingly, there is also a significantly higher signal in the skin of non-pigmented B6 mice than in albino mice ($p<0.001$). There was no significant difference between the photoacoustic signal, taken at 700 nm, in the body of the albino mice compared to the non-pigmented B6 mice ($p=0.95$), however, there is a significantly lower signal in the body of the pigmented B6 mice ($p<0.001$) [Fig. S9(B) in the Supplementary Material]. The photoacoustic signal within the skin shows a negative correlation with the photoacoustic signal in the body of the mouse ($p<0.001$) [Fig. 7(a)]. Unmixed estimates of blood oxygen saturation are significantly higher [Fig. 7(b)] in pigmented mice compared to non-pigmented ($p=0.025$) and albino mice ($p=0.019$). There is a strong positive correlation between the skin photoacoustic signal and the ${\mathrm{sO}}_2^{\mathrm{EST}}$ in the body of the mouse [Fig. 7(c), $p<0.001$], demonstrating that a higher melanin pigmentation in the skin leads to an overestimation of the underlying ${\mathrm{sO}}_2$. Comparing these results to Fig. 6(e), on the other hand, where negative pixels were included, appears to show a decrease in ${\mathrm{sO}}_2^{\mathrm{EST}}$ in pigmented mice compared to albino and non-pigmented mice across the single line plot. These findings highlight the potential influence of different approaches to image reconstruction and quantification on the observed bias.

Fig. 7

Increased skin pigmentation reduces the photoacoustic signal in the body, and linear unmixing estimates of blood oxygenation are increased in pigmented mice. (a) The mean skin photoacoustic signal negatively correlates with the mean body photoacoustic signal ($p<0.001$). (b) ${\mathrm{sO}}_2^{\mathrm{EST}}$, averaged across the body of B6 mice is significantly higher than that of albino mice ($p=0.025$), and pigmentation is a significant factor in the linear model ($p=0.019$). (c) ${\mathrm{sO}}_2^{\mathrm{EST}}$ averaged across the body region, correlates with skin photoacoustic signal ($p<0.001$). Averages were taken across the regions defined in Fig. 6.

Learned spectral unmixing was applied to the mouse data, with fair performance in albino and non-pigmented B6 mice. In pigmented mice, no clear improvement was observed over linear unmixing, with reconstruction artifacts dominating the oxygenation maps (Fig. S10 in the Supplementary Material).