Clinically-translatable High-delity Photoacoustic Tomography Enhanced by Virtual Point Sources

Photoacoustic tomography (PAT), a hybrid imaging modality that acoustically detects the optical absorption contrast, is a promising technology for imaging hemodynamic functions in deep tissues. Particularly, PAT is capable of measuring the blood oxygenation level using hemoglobin as the endogenous contrast. However, the most clinically compatible PAT configuration usually employs a linear ultrasound transducer array and often suffers from the poor image fidelity, mostly due to the limited detection view of the transducer array. PAT can be improved by employing highly-absorbing contrast agents such as droplets and nanoparticles, which, however, have low clinical translation potential due to safety concerns and regulatory hurdles. Moreover, unlike hemoglobin, these exogenous contrast agents cannot report the functional hemodynamic information. In this work, we have developed a new methodology that can improve PAT’s image fidelity without hampering its functional capability or clinical translation potential. By using clinically-approved microbubbles as virtual point sources that strongly scatter the local pressure waves generated by surrounding hemoglobin, we can overcome the limited-detection-view problem and achieve high-fidelity functional PAT in deep tissues, a technology referred to as virtual-point-source PAT (VPS-PAT). We have thoroughly investigated the working principle of VPS-PAT by numerical simulations and phantom validations, showing the acoustic origin of signal enhancement and the superiority over traditional PAT. We have also demonstrated proof-of-concept applications of functional VPT-PAT for in vivo small-animal studies with physiological challenges. We expect that VPS-PAT can find broad applications in biomedical research and accelerated translation to clinical impact.

dependence on the vessel orientations. For an exact image reconstruction, PAT requires full record of the pressure wavefront travelling in all orientations, i.e., a solid angular detection aperture of 2π for a planar detection geometry or 4π for a spherical or a cylindrical detection geometry [24], [29]. When a linear transducer array with limited detection aperture is used, only partial wavefronts from the blood vessels can be recorded, which results in the limited view problem.
Current solutions for the limited view problem in PAT include increasing the detection aperture, employing exogenous highly-absorbing contrast agent, and reconstructing with deep learning approaches. Increasing the detection aperture is usually achieved by either rotating the linear transducer array to cover a large detection angle [30], [31], or using a ring-shaped or spherical transducer [32], [33]. The former method requires complex system design, and the later greatly increases the cost of the system. Both methods would substantially hamper the clinical translation. The limited view problem can also be addressed by introducing sparsely-distributed exogenous contrast agents such as gold nanoparticles, which have much stronger optical absorption than hemoglobin and can disrupt the homogeneity of the optical absorption inside the bloodstream. The exogenous contrast agents can generate spherically travelling pressure waves and be well detected by the linear transducer array. Blood vessels in arbitrary orientations can be imaged by accumulating sufficient number of frames [28], [34]. The major limitation of this approach is that the received PA signals are from the exogenous absorbers, and thus physiological information such as oxygen saturation of hemoglobin (sO2) cannot be estimated. In addition, none of these exogenous contrast agents are approved by FDA for clinical applications with unknown safety effects and difficult regulatory clearance. Alternatively, various deep-learning networks have been used for artifact removal in linear-array based PAT, which have shown improved image fidelity [35]- [41]. However, the deep-learning based methods may introduce overfitting artifacts for data with significant noise, and their in vivo performance are mostly unknown [42]. Therefore, there still lacks a clinicallytranslatable PAT technology with high imaging quality and functional imaging capability.
In this work, we have demonstrated a novel method to overcome the limited view problem and improve the image fidelity of clinically-translatable PAT, by using microbubbles as virtual point sources (VPS), referred to as VPS-PAT. In VPS-PAT, although the microbubbles themselves have close to zero optical absorption, they can redirect the local PA pressure waves travelling in any directions via Rayleigh scattering and emit spherical waves that can be well detected by the linear transducer array. By reconstructing the sparsely distributed microbubbles as VPS over a number of imaging frames, VPS-PAT can produce high-fidelity images of the vasculature with improved visibility of the vertical structures and reduced reconstruction artifacts. Importantly, since the detected signals from the VPS originate from the surrounding hemoglobin, VPS-PAT keeps the functional sensitivity to the blood oxygenation level. Moreover, since microbubbles have been routinely used in clinical ultrasound imaging, we expect that VPS-PAT should have accelerated clinical translation compared with other PAT technologies relying on highlyabsorbing contrast agents. We have thoroughly investigated the working principle of VPS-PAT by numerical simulations and phantom studies, showing its superior imaging performance over traditional PAT. We have applied VPS-PAT for small animal studies to demonstrate its functional imaging capability. Our experimental results have collectively illustrated that VPS-PAT can effectively improve the image fidelity of clinically-translatable PAT, and may find broad applications in preclinical research and clinical practice.

Imaging principle of VPS-PAT
FDA-approved microbubbles are usually made of a gas core (e.g., nitrogen or perfluorocarbon gases) encapsulated by a solid shell (e.g., lipid or protein). Microbubbles have been widely used as contrast agents in clinical ultrasound imaging, in which they can enhance the blood vessel signals by strongly scattering the transmitted ultrasound waves [43]- [46]. Typically, native microbubbles without being conjugated with other chromophores cannot be used as contrast agents for PAT, because microbubbles have negligible optical absorption compared to hemoglobin (Supplementary Fig. 1). Unlike the highly-absorbing nanoparticles and droplets, microbubbles may serve as negative contrast agents in the bloodstream and break the homogeneity of the optical absorption.
However, using numerical simulations, we have demonstrated that, due to the small size (~5 μm) and sparsity of the microbubbles (4×10 6 microbubbles/ml), the zero optical absorption of microbubbles has practically negligible impact on the generated PA signals and thus has little contribution to mitigating the limited view problem ( Supplementary   Fig. 2) [47], [48].
However, since microbubbles have large acoustic impedance mismatch with the surrounding bloodstream [49], we expect that the strong acoustic scattering of microbubbles can play a much more important role in addressing the limited view problem in PAT. The most common PAT system usually uses a linear array transducer (Fig. 1a).
The linear array transducer typically has a central frequency between 2-10 MHz, corresponding to an in vivo acoustic wavelength of 150-750 µm [50]. Microbubbles used for clinical applications typically have a diameter between 0.5 and 10 µm, which is much smaller than the PA pressure wavelengths and thus can be considered as point targets [49]. Similar to contrast-enhanced ultrasound imaging, sparsely distributed microbubbles can function as Rayleigh scatterers and redistribute the surrounding PA pressure waves travelling in all directions (Fig. 1b, Supplementary Fig. 3). Because Rayleigh scattering has low anisotropy (or directivity dependence), the redirected PA waves by microbubbles are close to spherical waves, and thus can be readily detected by the linear transducer array positioned at any orientation. By filtering out the static PA signals directly from the blood via single value decomposition (SVD), VPS-PAT detects only the scattered PA signals from the microbubbles, and then reconstructs the images of the microbubbles as virtual point sources (Fig. 1c). By accumulating the microbubble signals from multiple independent frames, VPS-PAT is able to reconstruct the otherwise invisible vessel structures and thus improve the image fidelity.

Simulation of VPS-PAT using microbubbles as virtual point sources
To demonstrate the concept of VPS-PAT, we performed 2D k-wave simulation of PA imaging of a vertical blood vessel inside a mouse liver model [51]. The simulation covered a field of view of 12 × 15 mm 2 region that included skin, muscle, fat, and liver, with the linear array transducer placed on top of the skin surface. A vertical blood vessel (5 mm length, 0.2 mm diameter) was placed 11 mm from the transducer surface (Fig. 2a). The acoustic properties were heterogeneous across different tissue types in speed of sound, tissue density, and acoustic attenuation coefficient (Supplementary Fig. 4) [52], [53]. We assumed that optical fluence was homogenous around the blood vessel and hemoglobin was the dominating light absorber. A time-reversal based method was used to reconstruct PA images from the simulated PA data [54]. The image reconstruction also assumed homogenous tissue density and sound speed. The microbubbles had a diameter of 5 µm, an optical absorption coefficient of zero, and an acoustic impedance of 500 kg/(m 2 s).
Without any microbubbles inside the vertical vessel (Figs. 2b), the forward simulation shows that the generated PA pressure waves propagated sideways along the vessel (Supplementary Video 1). Thus, only the top and bottom boundaries of the vessel were visible in the reconstructed PA image, while the middle portion between the two boundaries was missing, demonstrating the classic limited view problem (Fig. 2c). With a single static microbubble inside the vessel (Fig. 2d), the PA signals redirected by the bubble travelled as a spherical wave and was detected by the transducer array (Supplementary Video 2). The reconstructed image clearly shows the microbubble as a point source at the correct position (Fig. 2e). Similarly, with five static microbubbles sparsely distributed inside the vessel (Fig. 2f, Supplementary Video 3), their positions were all correctly reconstructed in the PA image (Fig. 2g), showing partially the vertical vessel lumen. We further analyzed the raw PA signals detected by the transducer in the three different scenarios, and the results show that the scattered signals by the microbubbles were ~3-4 times weaker than the boundary signals (Figs. 2h).
Consequently, in the reconstructed images, the microbubbles as the virtual points sources were also ~3-4 times weaker than the vessel boundaries (Fig. 2i). Nevertheless, the microbubbles enabled better visualization of the vessel lumen with a 3-fold increase in the contrast to noise ratio (CNR) with one microbubble, and a 5-fold increase with five microbubbles (Fig. 2j) [55], [56].
We then simulated VPS-PAT with flowing microbubbles in the bloodstream. The flowing microbubbles were sparsely and randomly distributed in blood vessels. When a blood vessel was imaged over multiple frames with a relatively low frame rate, each frame had a different set of microbubble locations. By accumulating a certain number of frames and merging the reconstructed locations of all virtual point sources, we can recover the full vertical vessel structure that is otherwise invisible in traditional PAT. To mimic microbubbles flowing in bloodstream, we simulated and reconstructed 20 different frames.
Each frame had five microbubbles defined at random locations within the vertical blood vessel (Fig. 2a), which corresponded to a concentration of 2.2×10 6 particles/ml or ~0.03% volume concentration. This concentration was slightly lower than the typical clinical dose in contrast-enhanced ultrasound imaging [47].
As the scattered signals from the virtual point sources were weaker than the nonscattered signals from the vessel boundaries, we investigated different signal processing techniques to suppress the non-scattered signals. We first applied singular value decomposition (SVD) to the reconstructed images and removed high singular values with an empirically optimized cutting threshold (Fig. 3a). The non-scattered signals from the top and bottom boundaries were effectively filtered out due to the high temporal coherence over different frames. Compared to the original B-mode images (Figs. 3c-d), the scattered signals from the VPSs were largely preserved, because of their low temporal coherence from frame to frame (Figs. 3e-f). By merging the SVD-filtered images and accumulating all the reconstructed VPSs as the final image, the visibility of the vertical vessel as well as the CNRs were greatly improved (Figs. 3i-k). We also explored another signal processing method by applying a high-pass filter on the RF data across consecutive frames (Fig. 3b), which highlighted the scattered signals from the VPSs (Fig.   3g). The high-pass-filtered RF data was then reconstructed to form the images of the VPSs, revealing similar CNR improvements as the SVD-filtered image (Figs. 3h-k).
However, the high-pass-filter-based image presented more reconstruction artifacts than the SVD-filtered image, likely because the SVD method was more tolerant to the multiscattering signals by considering the signal's spatial coherence as well.

VPS-PAT for improving visibility of vertical structures
We validated the concept of VPS-PAT by imaging a loose knot made of a silicone capillary tube with a 300-μm inner diameter. The knot was embedded in 1% transparent agar and aligned with the linear transducer's imaging plane. We delivered pulsed laser excitation at 532 nm with a pulse energy of 7 mJ through the linear bifurcated light guide (Fig. 4a).
Whole bovine blood was mixed with microbubbles to yield a final concentration of 2×10 6 microbubbles/ml. The blood-microbubble mixture was perfused into the capillary knot with a flow speed of 20 mm/s, using a syringe pump to mimic blood flow (Fig. 4b,   Supplementary Video 4). A total of 500 consecutive frames were acquired with a frame rate of 10 Hz. As a control, we also imaged whole bovine blood without microbubbles (Fig.   4c, Supplementary Video 5). The same SVD filter was applied to both experiments. The US Doppler imaging was also performed on the same structure as a reference (Supplementary Video 6).
As the knot phantom had both horizontal and vertical segments, it was an ideal target for demonstrating the effectiveness of VPS-PAT in visualizing arbitrary target orientation. In the B-mode PA image, the vertical segments of the knot were poorly visible, as scattered signals from the microbubbles were overshadowed by non-scattered signals (Fig. 4b).
As expected, the SVD filter was able to suppress the non-scattered signals directly from hemoglobin and recover the correct shape of the knot (  (Fig. 4f). We also extracted the image intensity profiles in both lateral and axial directions of the SVD-filtered images. With microbubbles, vertical sections of the knot were clearly distinguished from the background (Figs. 4g-h). CNR was calculated at four representative regions, and microbubbles helped improve the CNR for both horizontal and vertical regions of the knot (Fig. 4i). This phantom experiment shows that microbubbles as VPS were able to improve the visibility of the vertical structures, and thus mitigate the limited view problem in PAT.

VPS-PAT for measuring oxygenation of hemoglobin (sO2)
In addition to improving the visibility of vertical structures, the scattered PA signals from the microbubbles preserve the functional information of the surrounding hemoglobin and sO2 can be estimated with multispectral PA imaging, which is fundamentally different from other highly-absorbing exogenous contrast agents. To demonstrate the functional sensitivity of VPS-PAT, we prepared fully oxygenated and deoxygenated bovine blood by adding sodium bicarbonate (37 mg/ml of blood, followed by 5-minute 100% O2 infusion) and sodium dithionite (2.5 mg/ml of blood, followed by 5-minute 100% CO2 infusion), respectively (Fig. 5a) [57], [58]. The oxygenated and deoxygenated blood were expected to have nearly 100% and 0% sO2 respectively. The oxygenated and deoxygenated blood samples were then mixed with microbubbles and perfused into a capillary tube with a constant flow speed of 20 mm/s (Fig. 5b). The capillary tube was imaged by VPS-PAT with dual-wavelength excitation at 532 and 1064 nm (Figs. 5c-f). The sO2 image was calculated from both the B-mode and the SVD-filtered images using the linear spectral unmixing method [59]- [61]. Here, we did not consider the spectral coloring effect in the transparent agar phantom [62]. The results clearly show that the B-mode PA images were not able to recover the sO2 levels in the vertical segments of the phantom (Figs. 5c and   5e). By contrast, VPS-PAT was able to accurately quantify the sO2 levels of the blood phantoms (Figs. 5d and 5f). It is worth noting that the sO2 estimation by VPS-PAT was consistent at all segments along the capillary tube.

VPS-PAT for in vivo rat liver imaging under hypoxia challenge.
To validate the in vivo performance of VPS-PAT, we imaged the liver vasculature of a rat under normoxia and hypoxia. We placed the rat in the prone position on the imaging window, with the laser light illumination (1064 nm) and the linear transducer array immersed in coupling medium underneath the imaging window (Fig. 6a). The same tissue plane around the liver region was imaged before and after microbubble injection. A total of 1000 frames were acquired with a frame rate of 50 Hz and a data acquisition time of 20 seconds. The reconstructed PA images were first stabilized to eliminate the respiration motion artifacts, and then processed with the SVD filter to extract the scattered signals from the microbubbles. Here, we used a dual speed-of-sound (SOS) image reconstruction method, because of the relatively large difference in SOS between the coupling medium and the rat tissue (Supplementary Fig. 5). The dual SOS image reconstruction was able to better converge the relatively weak PA signals scattered from microbubbles.
The results have clearly demonstrated that VPS-PAT was able to image more blood vessels with improved contrast and reduced artifacts, when compared with the B-mode PA images (Figs. 6b-c). The VPS-PAT image showed similar anatomical structure as the Doppler ultrasound image (Fig. 6d) and revealed clutters of small vessels as well as vertical vessels inside the liver region, with a penetration depth of ~1 cm (Figs. 6e-f). In order to highlight the necessity of microbubbles, we also tested the imaging process on the same animal without the microbubbles. In this case, the flowing blood may still provide some PA signal fluctuations over multiple frames, due to the rearrangement of red blood cells from frame to frame [63]. Without microbubbles, the SVD process had reduced effect on improving the image quality (Supplementary Fig. 6), with weaker CNR and no major vessels beyond 5 mm depth.
To demonstrate the capability of functional imaging with VPS-PAT, we conducted hypoxia challenge on the rat and monitored the hemodynamic response around the liver region.
Here, we used a single wavelength at 1064 nm for the relative oxygenation measurement, because deoxy-hemoglobin has a 20 times lower absorption coefficient at 1064 nm than that of oxygen-hemoglobin. The PA signal change at 1064 nm should mostly reflect the relative change in the oxy-hemoglobin concentration. We alternated the normoxia and hypoxia states for three cycles, with 1000 frames acquired for each condition in each cycle (Fig. 7a). The hypoxia challenge induced a systemic decrease in blood oxygenation level, and therefore, we observed a clear decrease in total PA signal and the vessel visibility, when compared to the normoxia conditions (Figs. 7b-c, Supplementary Fig. 7,   Supplementary Videos 10). Again, the results show that the SVD-filtered images could visualize more vessels at larger depths under both conditions than the B-mode images, highlighting the fact that the SVD filtered signals were redirected from the surrounding hemoglobin by the microbubbles. We further analyzed regional signal intensities for vessel clutters (Fig. 7d) and signal profiles for individual vessels (Fig. 7e). The results showed that hypoxia led to a >65% drop in the overall PA signal intensity at 1064 nm, reflecting a decrease in the blood oxygenation level. We also quantified the relative signal changes from normoxia to hypoxia in the B-mode images (Fig. 7f) or SVD-filtered images (Fig. 7g). The different signal changes under hypoxia challenge were indicative of the original blood oxygenation level under normoxia condition [64].

Discussion
We have developed a novel imaging technology VPS-PAT that uses acoustically- We applied SVD filtering to remove the strong non-scattering signals to improve the CNR of virtual point sources. The SVD-filtered images show improved structure visibility and reduced limited-view artifacts, both of which are due to the spherical waves scattered by microbubbles as virtual point sources.
One important limitation of the current VPS-PAT system is that multiple frames are needed to accumulate enough scattered signals from microbubbles, which leads to a relatively long data acquisition time and thus respiration-induced motion artifacts for in vivo studies. In our current setup, the imaging speed is further limited by the low pulse repetition rate of 10-50 Hz. We have shown that, with as few as 200 stabilized frames and a total acquisition time around 4-6 seconds, the visibility of major vessels and small vessel clusters were already improved compared to the B-mode images ( Supplementary   Fig. 8). Incorporating more frames could help improve the SNR of SVD-filtered image with more averaging, but also introduced respiration-induced artifacts that led to blurring of small features. One improvement of VPS-PAT is the incorporation of laser sources with a higher PRF to improve the imaging speed. High imaging speed will also significantly reduce the motion susceptibility, and may even allow tracking VPS to measure the blood flow speed [28], [65]. As our method relies on the relatively weak PA signal scattered by microbubbles, high SNR would help separate the VPS from the background. In deep tissues with low optical fluence and thus low SNR, advanced optical acoustic detection method with improved detection sensitivity may be applied [66]- [68].
For image reconstruction, we adapted a dual-SOS based method that was superior over the traditional single-SOS based method (Supplementary Figure 5). Although the dual-SOS method did not significantly improve the B-mode PA images, it was critical for improving the SVD-filtered VPS images that were sensitive to the signal convergence of point targets. Nevertheless, we still assumed all types of tissues had the same speed of sound, which was not accurate especially for tissues with high fat content and fibers [52].
A more accurate reconstruction method with tissue-dependent speed of sound, such as fast marching method, can help improve the VPS signal convergence and thus further improve the detection sensitivity and SVD-filtered imaging quality [69]. Moreover, instead of global SVD filtering, block-wise SVD filtering with adaptive cutoff values will improve the visibility of small vessels as well as noise suppression [70].
Compared with previous methods using ring-shaped or spherical transducer arrays, our method using microbubbles as virtual point sources has much lower cost and computational burden. Compared to exogenous contrast agent for improving structure visibility, microbubbles have negligible optical absorption. Microbubbles only redirect the signals from the surrounding blood and thus preserve the functional information, which allows estimation of physiological parameters such as sO2 and blood flow speed.
Microbubbles act as virtual point sources, which potentially allow super-resolution imaging via high-precision localization. Most importantly, as microbubbles are FDA approved contrast agents for US imaging, our approach is more clinically translatable than other preclinical contrast agents.
In summary, we expect that our new method VPS-PAT can provide high-fidelity hemodynamic information for preclinical and clinical applications such as disease monitoring, drug assessment, and surgical guidance.

Initiative Grant on Deep Tissue Imaging 2020-226178 by Silicon Valley Community
Foundation. The authors would like to thank Dr. Pengfei Song for the useful discussion and the engineering team at Sonovol for the technical support.

Competing interests
The authors declare no competing interests.

Additional information
Supplementary information is available for this paper.

Photoacoustic tomography
In PAT, photons in pulsed laser are absorbed by chromophores in tissue, and their energy is completely or partially converted to heat. The transient temperature rise leads to increased pressure and subsequent sound wave generation. The amplitude of pressure increase, as well as the intensity of the sound wave, is directly proportional to the chromophore concentration and local optical fluence. The sound wave is then detected by an ultrasonic transducer and used to form an image to tomographic map the original energy deposition in the tissue. Using acoustic detection, the penetration depth of PAT is less limited by optical scattering and is comparable to that of ultrasound imaging.

VPS-PAT system
The VPS surface, which was way below the ANSI laser safety limit. The channel data from the array transducer was sampled at 60 MHz and saved for further processing.

Ultrasound Doppler imaging
The same commercial ultrasound scanner (Vantage 256, Verasonics, Washington) and the L7-4 linear array transducer were used to acquire in-phase/quadrature (IQ) data for both real-time B-mode and ultrafast Doppler imaging of microbubbles. US plane waves were transmitted at a center frequency of 5 MHz with a PRF of 5000 Hz. Coherent spatial compounding was implemented with 7 steering angles (-3.5 degree to 3.5 degree with an increment of 1.16 degree), which led to a post-compounding frame rate of 500 Hz. For a single frame, each transmission event was separated by 80 μs to cover an imaging depth of 35 mm. An ensemble of 4000 compounded frames at the region of interest was acquired in 8 seconds and the IQ data were saved. The image data were stabilized using phase correlation-based sub-pixel registration [71]. Global singular value thresholding (SVT) was applied to the motion-corrected data to separate the microbubble signals from the tissue signals. As tissue signals typically resided in the large spectral components with high singular values, the SVT cutoff was determined by searching for the turning point which the singular value curve began to flatten [72]. The singular values larger than the cutoff were set to zero, and the ensemble after thresholding were averaged for further SNR improvement.

k-Wave simulation
The k-Wave toolbox was used for simulation [51]. In 2D k-Wave simulation, the field of view covered a 12 × 15 mm 2 with a grid size of 5 × 5 µm 2 . The microbubbles' location was randomly generated within the blood vessel region and was defined as a single pixel with acoustic impedance of 500 kg/(m 2 s). As the initial pressure in PAT is always smaller than 10 kPa, the microbubble oscillation is negligible, and we assumed no harmonics were generated. The sensor defined in the simulation had the same configuration as the L7-4 linear array, which contained 42 elements and 1/3 of the original aperture. The acquired PA channel signal was filtered using a 5 th -order Butterworth filter with a 4-7 MHz cutoff.
As the pressure generated during the PA imaging is usually less than 10 kPa, the oscillation and the diameter change of the microbubble was negligible [73].

Microbubble preparation
Perfluorocarbon-filled phospholipid microbubbles (VesselVue, SonoVol) with mean diameter of 4.37 µm were used for PAT. The vial had an initial concentration of 2.21 × 10 9 particles/ml and was mixed with PBS to yield a final concentration of 2 × 10 6 and 2 × 10 7 particles/ml for phantom and animal study respectively. A 27-gauge needle was used for injection via the tail vein to avoid microbubble bursting.

Singular value decomposition filter
Similar to slow-time signal in the US color Doppler image, the slow-time signal from PAT with microbubbles is comprised of three constituent components: 1) non-scattered signals, which are PA signals transmitted directly from blood and received within the detection aperture; 2) scattered signals, which are PA signals isotropically scattered from the microbubbles and received within the detection aperture; 3) random white noise, which originates from thermal and electronic noises [74]. For an image pixel that does not contain microbubbles, the slow-time signal should have only non-scattered signals and random white noise. A slow-time signal ensemble x can be expressed in following two forms for a give ensemble size ND: with microbubbles without microbubbles . (1) In (1), x(n) represents the n th slow-time data sample, and d, s, n are vectors of length ND for non-scattered PA signals, scattered PA signals, and white noise, respectively. Based on our simulation results, the non-scattered signal, usually from horizontal structures, is stronger than other signal components, including scattered PA signals from both vertical and horizontal structures. Therefore, we need to suppress the non-scattered component.
It is known that eigenvalues contain information about the signal energy represented by the individual eigen-components. In particular, a large eigenvalue signifies that its respective eigen-component accounts for a substantial energy portion of the slow-time ensemble. Such a representation matches the characteristics of non-scattered PA signals.
In order to use eigenvalue-based filtering, we first decomposed the slow-time signal contents into a sum of mutually orthogonal components that have the minimum meansquared modeling error. For a given slow-time ensemble x [75]- [77]: is the k th eigenvector and is the corresponding expansion weight that satisfies the following orthogonality relation: * 0 .
By identifying the eigen-components of non-scattered signals as the ones whose eigenvalues are above a predefined threshold λd, we can remove these eigencomponents and suppress the non-scattered signals. The resultant eigen-component subset for scattered siganls is For our application, the threshold λd was empirically set relative to the largest eigenvalue as λd ≈ 0.1λmax.

Experimental animals
Female Wistar rats (Charles River Laboratories:150 -200g, 6-10 months old) were used for the in vivo studies. The laboratory animal protocol was approved by the Institutional Animal Care & Use Committee of Duke University. During the experiment, the rat's temperature was kept at 37 °C by a heating lamp. Inhalant isoflurane was used for anesthesia with 5% v/v for induction and 2% v/v for maintenance. Compressed air flow rate through the facemask was set to 1000 mL/min, and the breathing rate of rat was around 55-60 times/min. Rat's fur at the liver region was removed with shaving cream.
The rat was placed in prone position on the acoustically and optically transparent imaging window and secured with surgical tape. A dose of 100 μL microbubbles (2 × 10 7 particles/ml) was injected via the tail vein before each data acquisition session. Each data acquisition session was separated at least 10 minutes apart to ensure that injected bubbles had been cleared out. Each injection corresponded to approximately 170,000 microbubbles per mL of blood, which was lower than the typical dose in clinical practice.

Hypoxia challenge
The rats were housed and anesthetized under normoxia (21% O2). During the hypoxia challenge, isoflurane vaporizer was constantly supplied with a mix of compressed air and pure N2 (volumetric flow rate 1:1) to achieve an oxygen level of 10.5%. Each cycle started with a bolus injection of microbubbles. Immediately after the injection, 500 PA frames were acquired and saved under normoxia condition, which took around 80 seconds. The gas supply to the vaporizer was then switched from normoxia condition to hypoxia condition. PA data was acquired 30 seconds after the hypoxia started, allowing the blood oxygenation level to stabilize under hypoxia condition. Each hypoxia session last 1.5 minutes, after which the rat was allowed to recover with normoxia condition for 10 minutes.
The hypoxia-normoxia switching was repeated 3 times.

Image reconstruction based on dual speed of sound
It is critical to correctly reconstructed the scattered PA signals from microbubbles. With the setup of our in vivo imaging system, the delay-and-sum (DAS) reconstruction based on single speed of sound was not accurate. We assumed the sound speed distribution in rat tissues was homogenous, but the difference between the speed of sound in the coupling medium and rat tissue was not negligible. Therefore, dual speed of sound reconstruction was used for our in vivo system. In order to perform dual speed of sound reconstruction, we needed the information of the speed of sounds in the coupling medium and the rat tissue, as well as the location of the imaging window.
We first measured the sound speed in the coupling medium under in vivo experiment condition. Three single hair strands were placed orthogonal to the imaging plane and directly on top of the imaging window. The hair stands were imaged, and the acquired data was reconstructed by sweeping the sound speed from 1250 m/s to 1500 m/s. The best image convergence of three hair targets indicated that the coupling medium had a sound speed of 1330 m/s.
We then determined the imaging window location. As our acoustically-transparent imaging window had negligible optical absorption at 1064 nm but strong absorption at 532 nm, we used 532 nm to image the window. After completion of rat liver imaging, we kept the transducer at the same location and removed the rat. The illumination was switched from 1064 nm to 532 nm and the acquired PA data of the imaging window was reconstructed with the sound speed calibrated from the previous step to provide the actual location of the window. As the temperature of the rat was 37 °C, we assumed the sound speed in rat tissues was 1540 m/s. With the location of the imaging window, the dual speed of sound image reconstruction was performed.

Motion detection and data stabilization
Due to the limited pulse repetition rate of the laser, the acquired in vivo datasets were susceptible to breathing motion. Due to rat's breathing, the tissues were periodically moving in and out of the imaging plane (Supplementary Video 9). To extract the stable frames, we first picked a reference frame and calculated the normalized 2D crosscorrelation between the reference frame and all the other frames in the acquired dataset.
The normalized cross-correlation reflects the breathing pattern of the rat, and by setting a threshold, the moving frames can be identified and removed.

100%. (5)
Where and are molar concentrations of oxygenated (HbO2) and deoxygenated (HbR) hemoglobin respectively. We used the linear spectral unmixing method, which states that P (λi, x, y), the reconstructed PA image at wavelength λi at (x, Here Φ(λ) is local optical fluence. and are the known molar extinction coefficients (cm −1 M −1 ) of HbR and HbO2 at wavelength λi, respectively. We used an optical energy of 15 mJ/pulse for 532 nm and 90 mJ/pulse for 1064 nm. After correcting for optical fluence, , and , can be estimated by solving the following linear equations: , , .
where P is the SVD filtered PA images at 532 and 1064 nm, ε is the molar extinction coefficient matrix as , , , , . .