2.1.

Endocavity Photoacoustic and Harmonic Ultrasound Imaging System Configuration

Our EPAUSI system comprises a clinical endocavity probe for light delivery [as depicted in Fig. 1(a)], which is connected to a combination of a research US system (Vantage 256, Verasonics, Washington, United States) and a cart-mounted tunable pulse laser (PhotoSonus, Ekspla, Lithuania) [as shown in Fig. 1(b)]. For both PA and US acquisitions, we employed a commercially available endocavity US probe (6EIX, Humanscan, Republic of Korea) with a center frequency of 6 MHz and a 95% bandwidth. This probe features a microconvex transducer array positioned along the periphery of a 10 mm-radius circular edge, offering a wide-angle 160 deg fan-shaped FOV beneath the application site. In light of the practical scenario involving direct contact with a canal wall, we implemented a dark-illumination scheme, ensuring that the transducer surface is not directly exposed to a laser beam reflected from the tissue surface. This dark-field illumination scheme presents several practical advantages. First, it eliminates the need for a transparent standoff layer filled with water or gel to clear the light path, resulting in the efficient reduction of the overall probe size by economizing the void volume. Moreover, the removal of the standoff layer helps maintain the signal intensity by reducing subsidiary attenuation and preserving the original FOV. Finally, this scheme prevents the direct exposure of the laser beam to the probe elements, mitigating strong photoinduced artifacts generated from the probe lens surface, which often interfere with PA contrast from microvessels at shallow depths. To further minimize PA artifacts arising from the transducer surface, we replaced the probe’s acoustic lens with a white one. To preserve the identical physical properties and geometrical form factors of the transducer, the replaced lens was crafted using the same elastomer and lens mold as employed in the original lens fabrication process, with the only variation being in color through adjustments in the pigment ratio.

Fig. 1

EPAUSI system configuration: (a) photograph and 3D rendering of EPAUSI probe assembled with optical casing. Optical and acoustic paths are depicted in red and blue, respectively. (b) System configuration of EPAUSI system on a cart. (c) Magnified rendering of EPAUSI probe head describing optical assemblies for line-beam shaping. (d) Photograph of a laser beam trajectory on the $x\text{-}z$ and $y\text{-}z$ planes imaged under a scattering suspension medium. PC, personal computer; OPO, optical parametric oscillator; DAQ, data acquisition unit; ED, engineered diffuser; LCF, large core optical fiber; ECP, endocavity ultrasound probe; and PC, probe casing.

For laser transport, we employed a large-core optical fiber, significantly reducing the size of the laser delivery module while achieving superior optical coupling efficiency (approximately 90%), in contrast to using a fiber bundle (approximately 50%). Specifically, a 1 mm core size multimodal fiber (Optibase, South Korea), equipped with air-gapped SMA-905 fiber connectors and bare FC/PC ferrules on each end, was coupled to a customized fiber coupling system within the portable OPO laser system (operating at a wavelength of 660 to 1064 nm, a pulse width of 3 to 5 ns, and a pulse repetition rate of 20 Hz). This fiber delivered the laser beam to the distal end of the endocavity probe, and the beam was directed into a line beam-shaping engineered diffuser (EDL-150, VIAVI Solutions, Arizona, United States) installed in front of the fiber outlet [as depicted in Fig. 1(c)]. The initially circular-shaped beam diverged into a single direction as it passed through the engineered diffuser, eventually transforming its pattern into a line shape [Fig. 1(d)]. These optical components were integrated within a biodegradable 3D printed casing, which was attached to the side of the probe head. To maximize the PA signal, the fiber channel within the casing was slightly tilted toward the imaging plane. Consequently, the total thickness of the probe head was 25 mm.

2.2.

Real-time, Video-rate PA and Harmonic US Image Processing

The signal acquisition, image reconstruction, and real-time display followed the event sequence timing outlined in Fig. 2. The acquisition of PA signals commenced with a laser trigger synchronized to a laser shot. This was swiftly followed by the acquisition of US signals, in which we implemented a pulse inversion harmonic US imaging technique to enhance the resolution and tissue contrast (Fig. S1 in the Supplementary Material). This method involves transmitting a pair of normal and inverted US pulses at each location, subsequently isolating the harmonic signal by canceling out the fundamental frequency signal. To match the probe’s bandwidth, we designed a pair of normal and inverted bipolar 4 MHz US transmit pulses, which ultimately resulted in an 8 MHz harmonic signal after the successful cancellation of the fundamental frequency signal by adding the radiofrequency (RF) signal pairs.

Fig. 2

The 20 Hz video-rate timing sequence of multispectral PA and harmonic US imaging. FWL, first wavelength laser trigger; PRF, pulse repetition frequency; DAQ, data acquisition unit; and RF, radiofrequency data.

For each US frame, 96 line-focused US beams swept across the entire FOV, with the transmit pairs being sequentially positioned to minimize tissue motion between the two acquisitions. In total, 1 PA receive event and 192 harmonic US receive events (comprising 96 pairs) were accommodated within the elapsed time between two laser emissions, amounting to 50 ms. Unlike US signal reception involving active pulse transmission, PA signal reception experiences greater attenuation in both optical and acoustical domains. Accordingly, an additional gain was allocated for the PA acquisition event, utilizing a programmable variable-gain preamplifier consisting of a time-gain compensator (TGC), programmable gain amplifier (PGA), and low noise amplifier (LNA) relay. Although the TGC was adjustable, it was maintained at the maximum value across all depths to ensure consistency and maximize the preservation of deep PA signals. Additionally, an extra 12 dB gain (PGA: 30 dB; LNA: 24 dB) was selectively applied in PA acquisition compared with harmonic US acquisition (PGA: 24 dB; LNA: 18 dB), with the preamplifier settings remaining constant throughout the experiments. The acquired PA and US RF data were transferred twice during a single frame acquisition to the host PC for online image reconstruction, utilizing a real-time Vantage reconstruction algorithm (Verasonics, WA, USA). Consequently, a pair of PA and harmonic US images was displayed on the screen at a real-time frame rate of 20 Hz. For data export and subsequent offline processing, such as energy density compensation and spectral unmixing (as detailed in Sec. 2.3), a conditional loop was activated. This loop involved an additional trigger event and a data save event. Operating in a multiwavelength mode, the laser system emitted an extra trigger of 8 microseconds in advance of the emission of the first wavelength of the laser pulse, allowing the US system to organize the data in the order of the laser pulse relay. After repeating the acquisition event as commanded by the user, the multi-frame PA RF, PA image, and harmonic US image data were saved, and the system returned to the real-time mode.

2.3.

Offline Angular Energy Density Correction and Spectral Unmixing

We observed an angular deviation in laser intensity between the center and the sides of the wide-angle FOV. Consequently, we conducted post-processing on the PA images to achieve uniform PA intensity along the angular direction. Utilizing the original PA image of a tissue-mimicking phantom (as described in Sec. 3.1), we estimated the angular energy density distribution by analyzing the PA intensity of radial-patterned string targets. Initially, the PA image was transformed from Cartesian coordinates to polar coordinates. Subsequently, point targets spanning between $\pm 70\mathrm{deg}$ and 15 to 25 mm in depth, along with their neighboring pixels, were windowed [as illustrated in Fig. 3(a)]. We calculated to the average amplitudes of the top 10% of pixels from each window, sorted them by depth, and generated an angular distribution of PA amplitude, as shown in Fig. 3(b). We fitted Gaussian curves to the measured data points, resulting in a consistent standard deviation of approximately 60 at all depths. Based on this observation, we created an estimated angular energy density map that covered the entire FOV. This map was generated by transforming the rectangular Gaussian distribution block into Cartesian coordinates, as depicted in Fig. 3(c). Subsequently, all experimental PA images were homogenized using inverse compensatory weighting based on the estimated angular energy distribution map.

Fig. 3

Estimation of angular energy density distribution using a tissue-mimicking phantom: (a) PA imaging result of ${\mathrm{TiO}}_2$-gelatin phantom (45 dB). To build a fluence distribution map based on angles, the image was changed from Cartesian coordinates to polar coordinates. (b) PA amplitude distribution by angles. (c) Estimated angular fluence distribution map.

We conducted linear spectral unmixing on the multispectral PA images to distinguish between multiple contrasts originating from distinct chromophores, such as clinical dyes and oxy and deoxyhemoglobin (as detailed in Secs. 3.2 and 3.5). Multispectral PA amplitudes at a specific location represent a combination of PA contrasts from various chromophores with distinct spectral properties, which are expressed using the following equation:

Considering the clinical scenario in which the probe is manually maneuvered by hand, we reduced the number of wavelengths ($n$) to maximize the spatial correlation among the PA images. To enhance the accuracy of ${\mathrm{sO}}_2$ measurements, we estimated an optical fluence for each original spectral PA image, and we performed preliminary compensation prior to its utilization as input for spectral unmixing. We divided the provided PA image into 1 cm depth increments, and within each step, we determined the mean value of the upper 50% of PA signals. These values were fitted with an exponential curve, assuming that optical fluence attenuation follows Beer–Lambert’s law. Distinct attenuation coefficients were obtained for each wavelength, and the respective PA images were inversely compensated by their corresponding exponential curves. Data post-processing was conducted using a custom script developed with MATLAB (MATLAB R2021a, MathWorks, Massachusetts, United States).

3.1.

Spatial Resolution Benchmarking in ${\mathrm{TiO}}_2$-gelatin Phantom

The measurement of the image resolution and angular energy density distribution was conducted using an in vitro experiment with a tissue-mimicking phantom. To create the phantom mold, transparent acrylic boards were designed, and holes were evenly punched on both sides of the side face in a radial pattern, spaced at 10 deg intervals in the angular direction and 5 mm intervals in the depth direction. A $100\mu \mathrm{m}$ black nylon string was threaded through these holes to serve as the imaging target. For the medium, a gelatin solution was prepared by adding 16.7 w% gelatin to water, which was thoroughly dissolved through heating at 60°C. The medium’s optical turbidity was achieved by adding $1\mathrm{g}/\mathrm{L}$ of ${\mathrm{TiO}}_2$ while constantly stirring.⁷⁰ This prepared medium was subsequently poured into the threaded mold and allowed to cool to room temperature, solidifying into a gel phantom. To ensure complete contact with the cylindrical surface of the probe, a groove with a 12.5 mm radius was created by dipping a cylindrical mold during fixation. For imaging, the assembled endocavity probe was acoustically matched with the phantom using a droplet of water and securely fixed upright onto the groove, with adjustments guided by a real-time harmonic US view. PA images were acquired at 700 nm with an energy density of $20\mathrm{mJ}/\mathrm{c}{\mathrm{m}}^2$, adhering to the American National Standards Institute (ANSI) criteria for maximum permissible exposure (MPE) at NIR wavelengths. As a result of this experiment, the simultaneous acquisition of PA and harmonic US B-mode images with a strong spatial correlation of the point targets was achieved, as illustrated in Fig. 4(a) and Fig. S2 in the Supplementary Material. The PA contrast appeared to be evenly distributed across the full 160 deg fan-shaped FOV scan range, similar to the harmonic US images, which were reconstructed with geometric homogeneity. However, contrary to the US image contrast, which extended down to 50 mm, the PA contrast logarithmically decreased with depth and reached the noise-equivalent depth at approximately 35 mm due to optical attenuation in the phantom media.

Fig. 4

EPAUSI spatial resolution measured from a tissue-mimicking phantom. (a) overlaid PA and harmonic US images highlight the spatially correlated contrasts from points in a radial pattern. Quantification of (b) radial and (c) angular resolution from line profiles of five 20-mm depth point targets, marked with dashed white box in panel (a). Mean ± standard error.

In terms of spatial resolution, we observed that the line profile of the point targets broadened with increasing the depth. The distinctive shape of the array dimensions, being convex, led us to measure the radial and angular resolutions in polar coordinates, in contrast to the lateral and axial resolution measurements in Cartesian coordinates. We obtained line profiles of the five central point targets at a depth of 20 mm, considering the elevational focal length of the probe. The maximum amplitude projections in the radial and angular directions were normalized, and the peak widths were measured following the full-width–half-maximum criteria [as depicted in Figs. 4(b) and 4(c)]. The radial resolution, which was narrower in harmonic US ($291\pm 22\mu \mathrm{m}$) compared with PA ($318\pm 16\mu \mathrm{m}$), can be attributed to the inherent disparity in carrier frequency between PA and US signals. PA signals are modulated by the fundamental frequency of the transducer (6.25 MHz), whereas harmonic US signals are modulated at the second harmonic frequency of the transmit pulse (8 MHz). The angular resolution of PA ($1.22\pm 0.02\mathrm{deg}$ and $639\pm 11\mu \mathrm{m}$ in lateral) was finer than that of harmonic US ($1.50\pm 0.09\mathrm{deg}$ and $785\pm 49\mu \mathrm{m}$ in lateral), as shown in Fig. 4(c). This observation can be explained by the greater limited view effect induced by the transducer’s geometry. Unlike the round trip travel of US signals, PA signals travel in one direction, and PA images are affected by the limited acceptance angle of the piezo elements. Additionally, the curved arrangement of piezoelements in microconvex probes results in a sparser acceptance angle overlay, contributing to a more pronounced limited view effect. Consequently, PA images may provide a better resolution in depicting microvessels, which can be particularly useful in identifying cancerous angiogenesis.

3.2.

Spectral Unmixing of Multiplexed Clinical Dyes In Vitro

To validate the molecular imaging capabilities of our EPAUSI system, we conducted a study in which we analyzed the constituent ratio of two Food and Drug Administration (FDA)-approved contrast dyes using multispectral PA imaging under in vitro conditions. Multispectral PA imaging, which leverages the distinctive optical spectral characteristics of chromophores in biological tissues, is a valuable technique for distinguishing the molecular-specific distribution of substances. This technique is widely employed in the current clinical translation of PA imaging. We prepared aqueous solutions of $65\mu \mathrm{M}$ methylene blue (MB) and $100\mu \mathrm{M}$ indocyanine green (ICG) and adjusted their concentrations to achieve absorbance levels comparable to that of hemoglobin in the NIR range, as documented on the Oregon Medical Laser Center’s spectral data repository.⁶⁹ To mimic the quantification of oxygen saturation, we mixed the MB and ICG solutions in volume ratios of 3:1, 2:2, and 1:3 [Fig. 5(a)]. Using a 27G syringe, these mixtures were subsequently injected into silicone microtubing threaded into the ${\mathrm{TiO}}_2$-gelatin phantom, as described in Sec. 2.1, with the MB content decreasing from left to right. To minimize variations in optical attenuation, all tubings were placed at the same depth (25 mm). For multispectral PA imaging, we selected 20 wavelengths at 2 nm intervals, ranging from 660 nm to 718 nm, as this range exhibited the most significant spectral absorbance variations for the two dyes.

Fig. 5

Differentiation of multiplexed clinical dye contrast with multispectral EPAUSI: (a) a photograph of five MB – ICG cocktail samples arranged in order of ICG constituent ratio. (b) An overlaid PA and US image of the microtubing phantom filled with five samples. Scale bar = 10 mm. (c) PA spectroscopy of five samples in the NIR range. (d) A series of PA images corresponding to four wavelengths (660, 680, 700, and 718 nm) involved in spectral unmixing. Multiplexed mapping of (e) a normalized MB concentration and (f) a normalized ICG concentration. Scale bar = 5 mm. (g) Quantified MB and ICG concentration from each corresponding tubing location. (h) Measured ICG constituent ratio and ground truth.

In the resulting PA and harmonic US images obtained at each wavelength, the upper and lower boundaries of the tubing were clearly visible in the harmonic US images, whereas the PA contrast of the dye appeared as small dots in the center of the lumen in the PA images [see Fig. 5(b)]. As the wavelength changed from 660 to 720 nm, distinct spectral trends in the PA contrast were observed within the tubing samples containing MB and ICG. Specifically, MB exhibited a consistent diminishing pattern, whereas ICG displayed a consistent increase [see Fig. 5(c)]. For practicality, we performed spectral unmixing using only four out of the total 20 multispectral image ensembles. These four images were selected from the quartile range at wavelengths 660, 680, 700, and 718 nm [see Fig. 5(d)]. The resultant unmixed distribution images of MB and ICG demonstrated a consistent rise and decline, respectively, which strongly correlated with the actual fractional composition of each dye [see Figs. 5(e) and 5(f)]. We quantified the upper 10% pixel intensity within the designated region of interest (ROI) within each tube. This analysis allowed for clear discrimination between MB and ICG in the pure tubing samples, whereas an equal constituent ratio was observed in the 2:2 cocktail tubing [see Fig. 5(g)]. Furthermore, when we calculated the fractional ICG concentration from the unmixed MB and ICG images, a monotonic increase was highlighted; it closely correlating with the actual fractional ICG concentration [see Fig. 5(h)]. However, the values for B and D, which should exhibit a volumetric ratio of 25% and 75%, respectively, appeared to deviate from the most representative reference ratio among the measured values. This discrepancy may stem from variances between the reference molar absorption coefficient used in the pseudoinverse calculations and the actual absorption coefficients of ICG and MB.

3.3.

ICG Infusion Kinetics into Ex Vivo Porcine Ovary

To assess the feasibility of using EPAUSI for transvaginal applications, we conducted real-time imaging of the dynamic ICG infusion into ex vivo female porcine genital tissue. In gynecological conditions necessitating surgical intervention, such as genital cancer or endometriosis, ICG is frequently employed for intraoperative fluorescent staining of malignant lesions, such as tumor masses or sentinel lymph nodes. For this experiment, we utilized a freshly excised porcine uterus, which included the ovary, fallopian tube (FT), uterus, uterine cervix, and vaginal canal [refer to Fig. 6(a)]. To replicate the clinical scenario, we enveloped the lower and surrounding regions of the excised ovary with FT tissue and subsequently covered them with an unfolded vaginal tissue flap. The EPAUSI probe was positioned over the vaginal tissue, and PAUS imaging was performed to penetrate through these tissue layers. Following probe positioning, harmonic US images clearly visualized the ovary at a depth of 23 mm, along with the adjacent FT tissues and the upper vaginal tissue layer. Subsequently, a 27G syringe needle was guided into the center of the ovary using real-time US imaging. Both 700 nm PA and harmonic US images were acquired at a rate of 20 Hz for a duration of 4 s, starting simultaneously with the injection of 3 mL of a $200\mu \mathrm{M}$ ICG solution through the needle.

Fig. 6

Ex vivo imaging of a porcine ovary with ICG infusion: (a) a photograph of a porcine female genital. A red inset image shows a magnified appearance of a porcine ovary in the end of a long FT. Comparative overlaid PA and harmonic US images of a tissue-surrounded ovary (b) in the preinjection state and (c) after ICG injection highlights the PA contrast of infused ICG (Video 1). (d) ICG infusion dynamics under the conditions of increasing PA amplitude inside the ovary. (e) A photograph of an ovarian hemisection after ICG infusion with topical ICG stain inside the ovary. UT, uterus; UC, uterus cervix; and VC, vaginal canal. Scale bar = 10 mm.

Before the injection, the syringe needle piercing through the ovarian membrane and into the ovarian center was clearly visualized in both the 700 nm PA and US images [see Fig. 6(b)]. However, upon the initiation of ICG perfusion, the PA image exclusively displayed the formation of a strong PA contrast cloud at the distal tip of the needle. This PA contrast area continued to expand over a span of 2 s, whereas such an observation was not evident in the US image [see Fig. 6(c) and Video 1]. We extended our analysis to compare the temporal changes in PA and US amplitudes. This involved selecting an ROI along the ovarian boundary and calculating the upper 10% average of pixel amplitudes from both images. As a result, the PA signal exhibited a remarkable increase of up to 10-fold within 2 s from the point of injection, compared with the pre-injection levels [see Fig. 6(d)]. After acquiring the images, the ovary used in the experiment was dissected into halves. The localized ICG stain was distinctly observable in a location corresponding to where the ICG cloud was evident in the PA image [see Fig. 6(e)].

3.4.

In Vivo Monitoring of Rat Prostate Vessel Pulsation with EPAUSI

We evaluated the in vivo imaging capabilities of EPAUSI by imaging the hemodynamics in the male reproductive organ of a rat. All experimental procedures involving animals were conducted in compliance with the protocol approved by the Institutional Animal Care and Use Committee of Pohang University of Science and Technology (Protocol POSTECH-2021-0052-C2, approved on February 17, 2022). A healthy 12-week-old male Sprague Dawley rat with a body weight of 347 g was initially anesthetized in an induction chamber with 3% to 5% isoflurane/${\mathrm{O}}_2$ at a flow rate of $1.5\mathrm{L}/\min$, using a preclinical anesthetic gas system (VIP 3000 Veterinary Vaporizer, Midmark, Ohio, United States). The rat was supplied with a continuous stream of anesthetic gas through a customized nasal cap. About 75% of its abdominal hair was removed using clippers and depilatory cream, and the rat was placed on a heating pad in a supine position to expose the lower abdomen. An EPAUSI probe was placed over the lower abdomen of the rat after applying ultrasonic gel for PAUS imaging. The probe was aligned on the sagittal plane with the left side directed toward the cranial direction [refer to Fig. 7(a)]. Below the probe, the male murine genitalia were located, as evidenced by the post-experiment autopsy image of the organ’s arrangement [refer to Fig. 7(b)]. The imaging plane encompassed the urinary bladder (UB), ventral prostate (VP), and urethra, with each organ being distinctly visible in the harmonic US image, showing unique echogenicity and structure [refer to Fig. 7(c)].

Fig. 7

Arterial pulsation in a murine prostate monitored in vivo with EPAUSI: (a) illustration of a transabdominal imaging setting of a male rat in a supine position. (b) Anatomical photograph of a murine male reproductive organ. The dotted line indicates the imaging plane of the probe. (c) Harmonic US B-mode image in the sagittal plane. (d) Overlaid PA and harmonic US image corresponding to (c). (e) Time-resolved normalized US and PA amplitude change in prostatic junction (yellow box in panels (c) and (d); Video 2). (f) US and PA temporal frequency spectra corresponding to (e). UB, urinary bladder; SV, seminal vesicle; VP, ventral prostate; DP, dorsal prostate; and USp, urethral sphincter. Scale bar = 10 mm.

Upon laser excitation, 700 nm PA and harmonic US images of the site were simultaneously displayed at a real-time frame rate of 20 Hz [refer to Fig. 7(d)]. In the center of both real-time views, we observed a repetitive pulsating motion at the boundary of the ventral and dorsal prostate junction at a depth of 8 mm (see Video 2). The prostate is a highly perfused organ with numerous interconnected blood vessels. Near its anatomical correspondence with the urethral axis, the identified vasculature appears to be a major branch from the central prostatic artery. Scattered PA contrast was also observed around the prostatic region, contributing to the subsidiary prostatic microvascular network. As we quantified the temporal amplitude changes in both PA and US images by comparing small regions of interest (ROI) at the pulsating site, we noticed synchronization between PA and US signals, evidenced by identical peaks [refer to Fig. 7(e)]. The signal changes were further analyzed by transforming them into temporal frequency spectra using fast Fourier transform (MATLAB R2021, Mathworks, Massachusetts, United States).⁷¹ These analyses revealed peaks at 4.7 and 9.4 Hz in both the US and PA spectra, representing the fundamental heartbeat frequency and its second harmonic [refer to Fig. 7(f)]. In summary, these findings confirm the practicality of our video-rate EPAUSI probe for capturing dynamic hemodynamics in live subjects.

3.5.

Hemoglobin Oxygen Saturation (${\mathrm{sO}}_2$) Mapping of a Human Blood Vessel In Vivo

We further demonstrated the feasibility of EPAUSI by imaging representative blood vessels from various body parts of a healthy volunteer. These experiments were conducted in strict accordance with protocols approved by the Institutional Review Board of Pohang University of Science and Technology (Protocol POSTECH-PRIB-2023-R009, approved on April 3, 2023). Safety goggles were worn by both the examinee and the examiner to ensure laser safety, and the laser density was controlled to not exceed $25\mathrm{mJ}/{\mathrm{cm}}^2$ at 700 nm, complying with the ANSI MPE criterion. Various body parts of the examinee, including the neck, upper limb, and lower limb, were imaged by placing the EPAUSI probe in cutaneous contact with US gel. This allowed for the visualization of subcutaneous major vessels such as the carotid artery, jugular vein, radial artery, and popliteal vein. We conducted 20 Hz PAUS bimodal imaging, with PA images acquired at different laser wavelengths (690, 730, 756, and 796 nm) in series. Out of the acquired multispectral PA images, we decoupled oxyhemoglobin and deoxyhemoglobin via spectral unmixing, referring to the molecular extinction coefficient of hemoglobin (source: Oregon Medical Laser Center, Oregon, United States). We calculated pixelwise oxygen saturation (${\mathrm{sO}}_2$) by dividing oxyhemoglobin by the total hemoglobin, and the ${\mathrm{sO}}_2$ map was thresholded based on the spatially coherent 866 nm PA pixel intensity. This thresholded ${\mathrm{sO}}_2$ map was overlaid onto the ultrasound image, providing both structural and photochemical features of the vessels. To quantify $\mathrm{sO}2$, we manually labeled ${\mathrm{sO}}_2$ pixels neighboring the vascular contrasts, and their upper 75% means were measured as nominal ${\mathrm{sO}}_2$ values. For each blood vessel, five independent measurements of ${\mathrm{sO}}_2$ were taken from consecutive ${\mathrm{sO}}_2$ frames, and the statistical presentation included calculating their mean values and standard errors.

Among the numerous PA images, we identified the deepest vascular contrast at a depth of 17.2 mm, originating from the lower boundary of a popliteal vein [Fig. 8(a)]. Both the upper and lower boundaries were confirmed in the PA image, and their appearances were also presented in the subsequent ${\mathrm{sO}}_2$ images at the same depth [Fig. 8(b)]. The nominal ${\mathrm{sO}}_2$ was measured to be $76.7\%\pm 1.4\%$, a value well within the known range of oxygen saturation in venous flow. Similar strong PA signals and highly correlating ${\mathrm{sO}}_2$ values were observed for the rest of the vessels [Fig. 8(c)]. In the case of larger vessels, such as the jugular vein or carotid artery, PA contrasts were only visible at the upper boundary of the vessels, indicating both strong optical attenuation from the blood within the lumen and an enhanced limited view effect due to the large vessel diameter. Interestingly, it was possible to classify the imaged vessels as arteries or veins based solely on their nominal ${\mathrm{sO}}_2$ values [see Fig. 8(d)]. The ${\mathrm{sO}}_2$ levels were measured to be high at 89.9% and 88.7% for the carotid artery and radial artery, respectively. By contrast, the veins exhibited comparatively lower ${\mathrm{sO}}_2$ levels, with measurements of 74.9%, 75.1%, 76.7%, and 77.2% for the jugular vein, radial vein, and two popliteal veins, respectively. The high consistency in ${\mathrm{sO}}_2$ levels within the same vascular group, both across measurements and between different vessels, underscores the robustness of multispectral EPAUSI. This suggests its potential applicability in measuring ${\mathrm{sO}}_2$ in deep reproductive organs, which could efficiently capture hypoxic lesions commonly found in reproductive diseases, such as cancer or endometrioma.

Fig. 8

Human in vivo transcutaneous multispectral EPAUSI: (a) overlaid harmonic US B-mode and 796 nm PA images of a popliteal vein at a depth of 17.2 mm. (b) A corresponding oxygen saturation (sO2) map overlaid on the harmonic US B-mode image. (c) Series of sO2 maps identified from major vessels in neck, upper limb, and lower limb. (d) Quantitative statistics of nominal sO2 classified into arterial and venous group (two arteries, four veins, $n=5$). CA, carotid artery; RA, radial artery; JV, jugular vein; RV, radial vein; PV1, popliteal vein 1; and PV2, popliteal vein 2. Scale bar = 10 mm.