1. Introduction

Blood oxygen saturation (sO₂) is a critical indicator in many clinical applications [1–3]. There are several imaging modalities that can assess sO₂. Functional magnetic resonance imaging (fMRI) is sensitive to the temporal changes of deoxyhemoglobin concentration, but it cannot measure the absolute value of sO₂ [4]. Positron emission tomography (PET) relies on exogenous tracers and poses a radiation risk [5]. Optical imaging methods like visible optical coherence tomography (OCT) and photoacoustic microscopy (PAM) can provide real-time high resolution sO₂ imaging, but they are limited by the low penetration depth of light [6]. Near infrared spectroscopy (NIRS) provides a practical means for imaging with centimeter-penetration, however it is of low spatial resolution [7].

Photoacoustic computed tomography (PACT) is ideal for sO₂ imaging thanks to its optical absorption contrast [8], ultrasound-based image reconstruction, and speckle-free nature [9]. Blood with varying levels of oxygen saturation exhibits distinct light absorption characteristics at wavelengths other than the isosbestic point [10]. As a result, PACT demonstrates sensitivity to the spatial variations of sO₂ within biological tissues. Furthermore, PACT can reach a penetration depth of several centimeters [11], enabling signal acquisition from depths much beyond the optical diffusion limit.

PACT has been widely used in both preclinical and clinical settings. In preclinical studies, PACT has provided insights into the mechanism of blood supply in the brain [12–14] and other body parts [15–17]. In clinical applications, PACT has been utilized for human organ imaging [18–21], brain functional imaging [21], and tumor oxygenation imaging [22–24]. For a more accurate quantification of sO₂, various methods have been involved for fluence compensation [25–28]. It was generally assumed that the sO₂ distribution within the cross-section of a single vessel is uniform. A previous animal study demonstrated a heterogenous distribution of sO₂ in a superficial vein using PAM [29]. In most cases, the imaging of sO₂ using PACT is accomplished through multispectral data acquisition. This process involves slow wavelength scanning, which hinders the assessment of the transient changes in sO₂.

PA imaging also possesses the capability to capture rapid hemodynamic processes occurring within blood vessels. This is facilitated by the high spatiotemporal resolution of PACT and its sensitivity to changes in hemoglobin. Recently, PACT has been demonstrated to provide a vector map of blood flow in human tissue. This new capability of flow measurement by PACT was based on the non-uniform spatial distribution of hemoglobin. However, there is a trade-off between the penetration depth and the range of velocity which can be estimated. Furthermore, it remained unclear how the inhomogeneity of sO₂ distribution contributed to the observed dynamic patterns [30].

In this study, we investigated the spatiotemporal distribution of blood oxygenation within a major vein and its branching vessels using PACT. All experiments were conducted on humans in vivo. We induced blood flow restriction through fist clenching and recorded variations in the PA signals inside the blood vessel. The analysis of multi-wavelength PA spectra within the vessel led us to the conclusion that the signal variations were attributed to the presence of deoxygenated blood flowing into the major blood vessels. To observe the blood flow pattern, we differentiated sequential image frames, allowing us to estimate blood flow velocity at various depths. Our findings suggest that PACT has the capability to visualize real-time distribution variations and has the potential to be applied in measuring hemodynamic parameters, such as flow speed and flow pattern, in regions with sO₂ heterogeneity.

2. Methods

2.1 Photoacoustic image contrast from sO₂ heterogeneity within a single vessel

Contrast of PA imaging comes jointly from the difference of optical absorption coefficient and changes of optical fluence. The PA signal ratio inside a vessel having two distinct sO₂ values can be expressed as

2.2 Numerical simulation

To verify that the experimentally observed pattern was due to sO₂ heterogeneity, we carried out numerical simulations in three dimensions on a digital phantom. A two-dimensional slice at the center of the phantom was shown in Fig. 1(a). Phantom vessel was placed inside a scattering media with the same optical scattering coefficients as human subcutaneous layer [31]. Only optical absorption of the phantom vessel was considered. The left half of the vessel exhibited an sO₂ level of 70%, consistent with typical values found in human large veins [32]. In contrast, the right half was assigned sO₂ values ranging from 10% to 50%. PA images of the digital phantom were reconstructed with a universal back projection method. We also calculated the theoretical PA signal ratio based on Eq. (1). The ratio is not unity, and is wavelength-dependent, due to the different sO2 levels between the left and right parts. The above calculation neglects the optical fluence variation inside the blood vessel. Furthermore, we considered the influence of the fluence variations on the PA signal ratio, based on a Monte-Carlo simulation.

 

Fig. 1. Numerical phantom of blood vessels containing blood with different sO₂ values. (a) 2D slice of the digital phantom for simulation described in 2.1. The phantom vessel had a heterogeneous sO₂ distribution. (b) Two-dimensional slice of the digital phantom of human forearm. The vessel (red oval) was assigned a heterogeneous distribution of sO₂. The upper part and lower part (illustrated with different colors) of the phantom vessel have different sO₂ value of 90% and 10% respectively, which was a rough estimation of the sO₂ distribution of the major vein we observed as shown in Fig. 6(a-c). Wavelength-dependent optical absorption coefficients in the 700∼900 nm range for all the layers are displayed. Scale bars: 2 mm.

Download Full Size | PDF

To analyze the PA spectra at multi-wavelength we got in the experiment, we adopted a digital phantom of human forearm as shown in Fig. 2(b) for simulation. Absorption coefficients of major absorbers including hemoglobin, water, melanin, and fat were set within the ranges as reported in [31,33,34]. The absorption coefficient of different skin layers varying from 700 nm to 900 nm were shown in Fig. 2(b) respectively. Wavelength-dependent scattering coefficients of epidermis, dermis, subcutaneous layer, muscle, and blood were set according to the values reported in [31,33]. We used the Monte Carlo method [35] to model light transport in the phantom. Two rectangular light sources symmetrically set on either side of the imaging plane were used for optical illumination. The phantom vessel was attributed to a heterogenous distribution of sO₂. A temporal increase in venous sO₂ after muscle contraction was previously reported in [36,37], rising up to an value of approximately 90%, promising that the sO₂ we set in the upper part of the phantom vessel was rational.

 

Fig. 2. (a) A schematic of the experiment system. DAQ: data acquisition unit, OPO: optical parametric oscillator, WT: water tank, TR: half-ring transducer array. (b) Photograph of the in vivo experiment setup. (c) 3D structure of the main blood vessels we observed in the experiment. Plane 1 refers to the imaging position in Fig. 4 and 5 using a single-wavelength acquisition mode. Plane 2 refers to the imaging position in Fig. 4 and 5 using a multi-wavelength acquisition mode.

Download Full Size | PDF

For the simulation of ultrasound propagation, we adopted a focused half-ring probe with 128 elements for ultrasound detection. The radii of the element curvature and the half-ring array were both set to 2.5 cm. All simulations assumed a homogenous medium without acoustic absorption or dispersion, and the speed of sound in the phantom was set to 1500 m/s. The electrical impulse response of our half ring PACT system was convolved with the sinogram generated by the k-Wave toolbox [38].

2.3 Imaging system and experiment procedure

We utilized a custom-built PACT system with a half-ring probe for the experiment as illustrated in Figs. 3(a) and (b). The system employed frequency-doubled Nd:YAG laser (Beamtech Opotronics Co. Ltd, Nimma-900) pumping an optical parametric oscillator (OPO, 50-80 mJ pulse energy in the 700-900 nm range, 10 ns pulse width, 10 Hz repetition rate) for PA signal excitation. A half ring transducer array of 128 channels (central frequency 2.5 MHz, 60% bandwidth; ULSO TECH) fixed on an adjustable positioner was used for signal detection. PA data was sampled at 80 MHz using a data acquisition unit (TsingPAI Technology Inc., Marsonix-128). For multi-wavelength excitation, the illumination wavelength was scanned per pulse from 700 nm to 900 nm at a 10 nm interval. The laser fluence applied to the skin remained below 12 mJ/cm² across the entire scanned wavelength range, significantly below the safety limit set by the American National Standards Institute (ANSI).

 

Fig. 3. The results of the theoretical analysis and simulation as described in Section 2.2. (a) Cross-sectional image at 700 nm of the vessel with sO₂ of 70% (left) and 10% (right). Averaged signal strengths in the rectangles were selected for calculating the PA signal ratio. Scale bar: 2 mm. (b) PA signal ratios obtained by simulation and theoretical estimation. The simulation results were obtained from the ratio between the averaged signal strengths of the right and left regions in the boxes in (a). The theoretical calculations disregarded the influence of optical fluence variations.

Download Full Size | PDF

One healthy volunteer (21 years old, Asian female) participated in the imaging experiments. We have obtained informed consent from the participant. We instructed the participant to horizontally immerse her arm in the water tank while maintaining an extended position. The transducer was also immersed in water for acoustic signal coupling. The room temperature was kept at 25°C and the water temperature was maintained at 32°C using a water temperature control system. The water temperature was kept relatively close to the body temperature in order to avoid vessel constriction induced by loss of heat [39]. We displayed a 3D modeling image of the main blood vessels observed in our experiment in Fig. 2(c). The model was built based on the 3D photoacoustic scanning results to reflect the exact structure of the blood vessels. We centered the field of view (FOV) on the major vein along with one of its branches near the elbow (as shown in Plane 1 in Fig. 2(c)), as well as a position adjacent to the initial location, approximately 1 mm away along the direction of venous blood flow (as shown in Plane 2 in Fig. 2(c)).

In the imaging experiment, we employed two conditions: fist clenching for durations of 20 seconds and 5 seconds. Fist clenching can cause several hemodynamic responses. It can lead to increased muscle activity and consequently a higher oxygen uptake of the muscles involved, which can generate venous blood with a lower sO₂ level during this process. The decrease in sO₂ is largely dictated by the exercise level of forearm muscle [37,40]. Previous study have also shown that contraction can accelerate blood flow in the forearm during the relaxation period compared to a pure resting state, and a longer contraction duration can induce a greater blood flow [36]. It is expected that a longer clenching time will contribute to more oxygen consumption and thus more deoxygenated blood will be generated. A more rapid blood flow is expected as well for a longer clenching period. PA data were collected at the elbow location of the large vein and its branching vessel following muscle relaxation.

During the experiment, the participant was asked to keep her fist clenched for 20 seconds and subsequently release it quickly. Single-wavelength PA data was acquired at 700 nm at Plane 1 in Fig. 2(c) during the recovery period. After the blood flow condition returned to normal, the participant was asked to clench her fist for 5 s and PA signal from the same location was acquired. We adopted a single-wavelength (700 nm) acquisition mode to acquire PA images for the aforementioned experiments, through which we can get a sampling rate of 10 Hz to capture the dynamics of blood flow alterations during the relaxation processes. We then moved the transducer 1 mm towards the proximal end (referred to Plane 2 in Fig. 2(c)) to acquire multi-wavelength data after 20 seconds of fist clenching. Since our experiment result exhibited a more stable blood flow in this case, we utilized the multi-wavelength (700 to 900 nm at 10 nm interval) imaging mode to acquire PA spectra for the subsequent analyses.

2.4 Image reconstruction

Before the reconstruction, PA signals were processed by a 0.1∼ 4.5 MHz bandpass filter to remove the unwanted out-of-band frequency noise. We used a model-based method with total variation (TV) regularization in the iterative steps for image reconstruction [41] for both the in vivo and simulation data. This method considered the geometry of the transducer elements, thereby introducing the spatial impulse response (SIR) into the model. Compared to the conventional back-projection method in 2-D reconstruction, model-based reconstruction with non-negative constraints can help reduce background noise and retain the sharp vascular edges [42], making it more suitable for the spectral analysis in our study.

2.5 Velocity calculation

Laminar flow and turbulent flow coexist in the venous lumen. For laminar flow, the velocity on the cross section of the blood vessel is parallel to the central axis and the velocity profile follows a parabolic shape. There is no mixing of blood on different streamlines across the lumen apart from that caused by diffusion. While for turbulent flow, the velocity exhibits a more complicated distribution and usually there is a thorough mixing of the fluid. In this study, we estimated the flow velocity in the branching vessel under the assumption of laminar flow since it manifested a profile stemming from the inflowing blood, which is characteristic of the laminar flow pattern. To estimate the flow velocity in the branching blood vessel (Fig. 5), corresponding points on the same streamlines of differential images were manually selected and flow velocity was calculated using the following equation:

 

Fig. 4. Imaging results from the 20-seconds fist clenching test. (a) Snapshots captured immediately after fist release. Photoacoustic images were acquired at 700 nm, with timestamps indicated in the images. MV: major vein, BV: branching vessel, SK: skin, SPV: superficial vessel. (b) Inter-frame changes within the branching vessel are highlighted in pseudo-color, overlaid on the original images. Scale bar: 2 mm.

Download Full Size | PDF

 

Fig. 5. Laminar flow observed in the branching vessel after 5-seconds first clenching. (a) Four consecutive frames with differential values inside the vessel overlayed on the original images to expose the dynamic flow pattern. Arrows indicated the blood flow inside the branching vessel. (b) Propagating wavefront from 0.1 s to 0.4 s overlayed on the original image at 0.1 s. Differential values were highlighted in various pseudo-color in order for better manifestation. Scale bars: 2 mm.

Download Full Size | PDF

3. Results

The results of our simulation to verify the image contrast induced by sO₂ heterogeneity are shown in Figs. 3(a) and (b). The PA signal ratio between the areas highlighted with the rectangular boxes in Fig. 3(a), as a function of wavelength, is illustrated in Fig. 3(b). The theoretical results were estimated based on Eq. (1), neglecting the optical fluence differences. It can be observed that at the wavelength of 700 nm, the largest signal ratio inside the vessel was obtained, though ${\mathrm{\epsilon }_{\textrm{Hb}}}$ reaches a peak at 760 nm. The contrast in the simulation results were lower than the theoretical ones since the optical fluence in the deoxygenated region was smaller because of its stronger light attenuation. According to the simulation and theoretical results, we selected 700 nm for the observation of transient flow pattern inside blood vessels in our in vivo human imaging experiments.

Figure 4(a) depicts the original images at 700 nm capturing the entry of low-oxygen blood into the vessel following 20 seconds of fist clenching. These images were obtained under a single-wavelength (700 nm) acquisition mode at Plane 1 in Fig. 2(c) to record the rapid inflowing process. To enhance the visualization of blood distribution with varying sO₂ levels within the vessel, we presented the frame-to-frame changes in pseudo-color overlaid on the original frames in Fig. 4(b). Differential signal below a positive threshold was set to zero to suppress the influence of noise. The positive signals in the differential images originated from blood with a lower sO₂ value, since the light absorption coefficient of Hb at 700 nm is higher than that of HbO₂. It can be observed that deoxygenated blood entered the main vessel (depicted as cross-sections in the images) from a position close to the central axis of the vessel. Subsequently, it diverged toward the vessel wall before flowing out into the branching vessel. Worth noting that oxygen consumption may not happen right at our imaging plane. The deoxygenated generated during the clenching process mostly came from capillaries of the forearm where oxygen exchange in capillaries is active, and flew into the major vein we observed after several times of convergence, so the mixing condition of oxygenated and deoxygenated blood cannot be inferred in our experiment. The corresponding video are provided in Visualization 1 with interpframe changes overlayed on the original images, and Visualization 5 displaying and entire blood flow process in the vessels.

Figure 5 shows the time-lapse images after five seconds of fist clenching. As the former in-vivo experiment, we also adopted a wavelength of 700 nm at Plane 1 in Fig. 2(c) for imaging. The color-coded differential images are overlaid on the original images and are shown in Fig. 5(a). To better visualize the wavefront evolution along the branching vessel, the color-coded wavefronts are co-plotted on top of a common PA image background in Fig. 5(b). The differential images revealed a distinct laminar flow pattern within the branching vessel. Due to the higher blood flow velocity at the center of the vessel, the differential images exhibited a characteristic parabolic profile, indicative of laminar flow. This pattern gradually elongated as it propagated down the branching vessel, as recorded in corresponding video provided in Visualization 2. The observed flow pattern is characteristic of venous blood [46]. Blood flow velocities at the center of the branching vessel, as well as at two positions symmetrically located relative to the central line, were estimated and are listed in Table 1.

Table 1. Estimated blood flow velocities at different radial positions inside the branching vessel.

View Table

In addition to the data collected at a single wavelength of 700 nm, we gathered data at 10 nm intervals from 700 nm to 900 nm. Through multi-wavelength imaging at Plane 2 in Fig. 2(c), we conducted PA spectral analysis to verify the components that contributed to the patterns inside the PA images of blood vessels. This set of data was obtained at a distance of approximately 1 mm along the direction of venous blood flow from the preceding cross-section. At this new cross-section, a more stable flow pattern which lasted for over 20 seconds was observed, with deoxygenated blood predominantly at the bottom of the blood vessels. At 700 nm, a strong signal was clearly visible at the bottom of the vessel as shown in Fig. 6(a), while at 900 nm (Fig. 6(b)), the signal in this region was consistent with the surrounding tissue, and no clear boundary can be observed. A close-up view of the boxed region in Fig. 5(a) is displayed in Fig. 5(c). The corresponding video acquired at 700∼900 nm can be found in Visualization 3. The curve shown in Fig. 6(e) displayed the PA spectrum of point 1 in Fig. 6(c) after averaging 10 frames. It can be seen that in point 1, there was a peak around 760 nm, with an overall decreasing trend from 700 nm to 900 nm. This curve is characteristic of the absorption spectrum of deoxygenated blood. For the upper part of the vessel (point 2 in panel (c)), the PA spectrum displayed in Fig. 6(f) showed an upward trend between 700 nm and 900 nm, indicative of a higher sO₂ value according to its spectral shape.

 

Fig. 6. Spectral analysis of a major vein’s cross-section. (a, b) Still frames of the same location obtained at 700 nm and 900 nm. (c) Closed-up view of the selected area in image (a). (d) Simulation result of the phantom vessel at 700 nm. (e, f) PA spectra of points 1 and 2, respectively, as labeled in (c). These spectra were averaged 10 times. The shaded area represented the level of uncertainty in the measured PA spectrum. (g, h) PA spectra of points 3 and 4, respectively, as labeled in (d). Scale bar: 2 mm.

Download Full Size | PDF

We further conducted numerical simulations using a digital phantom with the structure shown in Fig. 2(b). The phantom vessel was assigned a heterogenous sO₂ distribution, with an sO₂ level of 10% at the bottom while the remaining part had an sO₂ of 90%. More details of the simulation can be found in section 2.2. The simulation result at 700 nm was shown in Fig. 6(d). PA spectra of points 3 and 4 are illustrated in Figs. 6(g) and (h), respectively. These simulated spectral profiles are similar to the experimental ones shown in Figs. 6(e) and (f).

4. Discussion

Utilizing PACT, we investigated the spatiotemporal patterns within human blood vessels that reflect transients in the distribution of sO₂. The imaging experiment was conducted at the subject's elbow position using both single-wavelength and multi-wavelength, and changes in sO₂ were induced through fist clenching. The inherent sensitivity to spatially-varying optical absorption in PACT enabled the comprehensive visualization of the dynamic process of venous blood flow through branching veins.

The enhancement of intravascular PA signals at the bottom part of the vessel at 700 nm as shown in Fig. 6(c) can be attributed to several factors: (a) an elevation in deoxyhemoglobin concentration, (b) an increase in total hemoglobin concentration, (c) an increase of other chromophores with high absorption at 700 nm, and (d) the influence of spectral coloring [28].

We conducted an analysis of the PA spectra at different positions to illustrate that the observed signal changes within the blood vessels predominantly stem from variations in the distribution of sO₂. As shown in Fig. 6(d), we scrutinized the PA spectrum of the bright spot near the bottom of the blood vessel, a feature consistently observed throughout the image acquisition process. The spectrum at this particular location manifested the characteristic profile of deoxygenated blood, marked by a peak at 760 nm and an overall declining trend. In stark contrast to the aforementioned spectra, the PA spectrum of point 2, situated at the top of the blood vessel, displayed a characteristic of oxygenated blood. In conclusion, Fig. 6 illustrates a scenario wherein disparate sO₂ levels coexist within a vein at the same cross-section. We further conducted numerical simulations to confirm our hypothesis of the sO₂ distribution inside the blood vessel. As illustrated in Figs. 6(g) and (h), the PA spectrum of the upper part in the simulated vessel with an sO₂ of 90% exhibited an increasing trend, while the region with an sO₂ of 10% had a peak at 760 nm and an overall declining trend. The simulation results showed a similar trend as the experiment data, supporting our hypothesis of heterogenous sO₂ distribution.

The point where blood vessels intersect may harbor valve structures primarily composed of collagen, which has weak light absorption at the imaging wavelength. The streamlines of blood flow displayed distortion in the presence of valves, resulting in a displacement of the vertices of the parabolic wavefront onto a different streamline [30]. However, our analysis of the flowing pattern indicated that venous valves did not contribute to the observed transient pattern. As shown in Fig. 5(b), the vertices of the propagating wavefront align well within the same streamline, indicating an undisturbed flowing pattern and the absence of any nearby valves.

Through multispectral data and numerical simulations, we confirmed that the local signal enhancement and appearance of boundaries within the blood vessel were caused by deoxygenated blood. Previous studies have also shown that red blood cells at the bifurcation of blood vessels can have uneven distribution, leading to uneven photoacoustic signal strength [30]. Additionally, in the case of fist clenching, the concentration of total hemoglobin will also increase [36,47]. However, these findings are not contradictory to the conclusions we obtained regarding local changes in blood oxygen distribution through PA spectra analysis.

Using PACT imaging, we visualized the process of venous blood flowing inside blood vessels. Due to the sensitivity and high speed of the imaging system, we were able to record transient hemodynamic patterns immediately following blood occlusion with varying durations. Specifically, in Fig. 5(b), we depict a wave pattern characteristic of laminar flow, suggesting that, with appropriate signal processing, PACT can be valuable for analyzing flow speed and dynamic patterns. However, what we have demonstrated is the observation of a transient event involving a rapid influx of deoxygenated blood generated by fist clenching and sudden release. Therefore, we are cautious at this point about claiming a broader application in the general measurement of flow rates and flow patterns, and the method is currently more suitable for analyzing flow conditions with stable convergence of blood with different oxygen levels, as previously reported in [46]. Local changes in sO₂ are more likely to be caused by short-term changes in tissue metabolism and blood flow conditions [37,40]. Thus, PACT is potentially useful in visualizing local dynamics of metabolism [48].

Under normal blood flow conditions, the distribution of blood oxygen within a blood vessel is relatively homogeneous. However, we have observed dynamic patterns within major veins in humans under a resting state, presumably attributed to the temporary influx of venous blood with low sO₂ levels. The corresponding video is provided in Supplementary Visualization 4. We recorded the process of the evolution of an irregular pattern inside a vein at the middle of human forearm. No intentional exertion was conducted before or during this process, and the pattern inside the vein changed spontaneously. According to our observation, dynamic patterns usually appeared at venous junction, indicating a more complicated distribution of sO₂ and flowing mode.

In our current study, we utilized a pulse repetition rate of 10 Hz, which proved insufficient for capturing multispectral data during the transient flow process. And our sampling rate restricted us to use single-wavelength to record the blood flow variation shown in Fig. 4 and 5. However, by increasing the laser repetition rate, more reliable multispectral data can be obtained, potentially enhancing our analysis in multiple dimensions, including spatial, temporal, and spectral aspects simultaneously. Furthermore, in our study, we refrained from quantifying the sO₂ values because the PA spectrum inside blood vessels is distorted by the spectral coloring effect, making the conventional linear unmixing approach unreliable. For a more precise measurement of sO₂, advanced methods such as optical fluence compensation using deep-learning techniques can be employed [25–27]. However, even with these advanced methods, the heterogeneous sO₂ distribution encountered in this study poses significant challenges to the accuracy of sO₂ quantification, which is beyond the scope of our current investigation.

5. Conclusion

In this study, we documented the variations in PA signal strength inside a major vessel resulting from the transient flow of venous blood. We emphasize the capability of PACT to capture the time-varying, heterogeneous distribution of sO₂ within a vessel under circumstances such as blood occlusion and fast metabolic event. Further investigations could provide a more quantified estimation of the sO₂ throughout the entire vessel, linking sO₂ with blood flow velocity for more comprehensive understanding of intravascular hemodynamics.