Ultrasound-Guided Breath-Compensation in Single-Element Photoacoustic Imaging for Three-Dimensional Whole-Body Images of Mice

Photoacoustic imaging can be used to extract functional information at the molecular level for imaging the whole body of small animals in vivo. However, the use of a mechanical scanner to image the whole body involves acquiring the data for a period of time which can lead to breathing-related artifacts. Thus, the acquired three-dimensional data cannot be analyzed and visualized using two-dimensional projection rendering techniques unless the respiratory distortions are minimized. In this study, an ultrasound-guided breath-compensation method for the three-dimensional data of photoacoustic macroscopy to facilitate visualization and analysis in the depth direction is presented. Photoacoustic and ultrasound data of the whole body of mice were acquired in vivo, then the distorted skin layers were compensated in ultrasound data. The extracted distortion parameters were then applied to the corresponding photoacoustic data to compensate for the motion. The compensation method was successfully applied to visualize three-dimensional hemoglobin oxygen saturation in the whole body of mice in vivo by using multispectral photoacoustic data. The resulting three-dimensional images show that the developed methodology can be used in various biomedical studies, including monitoring drug delivery, imaging of tumors, and analysis of vasculature networks.


INTRODUCTION
Photoacoustic imaging (PAI) is currently attracting much attention in biomedical studies due to its unique ability to visualize functional information in biological tissues at the molecular level at ultrasound (US) resolution. PAI is derived from the photoacoustic (PA) effect, which is energy transduction from light to sound [1]. When a short-pulsed laser stimulates biological tissues, the endogenous chromophores (e.g., oxy-hemoglobin, deoxy-hemoglobin, water, melanin, and lipid) absorb the light energy, then transfer it to heat energy. The heat energy produces volumetric changes through thermoelastic expansion, generating acoustic waves (called PA waves) that propagate in all directions. The generated PA waves are detected by the conventional US transducer and reconstructed to form the PA image [2][3][4][5][6]. In addition to the endogenous chromophores, exogenous contrast agents also can be visualized in PAI [7][8][9].
Although PA and US images have similar data acquisition and image generation procedures, they provide different but complementary information on the underlying biological structure. The US imaging (USI) is based on the transmission and reflection of acoustic waves providing structural information, while PAI is based on optical absorption, representing functional information at the molecular level [10][11][12]. The two imaging modalities are typically implemented in a single imaging platform, especially in clinical applications [13][14][15][16][17][18][19][20][21][22][23], to complement each other facilitating better clinical outcomes.
PAI has a unique advantage in that the resolution and imaging depth can be scalable according to its application [24]. The two major types of PAI systems for preclinical application are PA microscopy (PAM) and PA computed tomography (PACT). PAM systems generally equip a single-element focused US transducer to acquire highresolution images (lateral resolution of~5-50 μm) [25][26][27]. Those imaging systems have been utilized to visualize superficial areas in small animals including the ear [28], brain [29][30][31][32], and eye [33][34][35]. In addition to the shallow imaging depth (~1-3 mm), one of the key issues in PAM is relatively slow imaging speed because of scanning of both the light and transducer to achieve images. In contrast, PACT systems can achieve images with faster imaging speed by using multi-element detection from various geometry of array transducers including linear, arc, and ring arrays [36][37][38]. The biodistribution of deep tissue in small animals has been investigated using various PACT systems [39][40][41][42][43]. However, the mathematical image reconstruction algorithms in PACT may generate artifacts, which degrade the quality of the resulting images. In addition, the high cost of the multielement US array and corresponding multi-channel data acquisition (DAQ) module may burden the researchers.
The mechanical scanning-based PA macroscopic (PAMac) system is one of the initial configurations of PAI for small animal studies [44]. PAMac systems have been used for imaging the whole body of small animals at the organ level with imaging parameters in between PAM and PACT (depth, axial resolution, and lateral resolution of~20 mm, 150 μm, and 590 μm, respectively) [45][46][47]. The main applications were monitoring the biodistribution of various agents to investigate biological responses of drugs, nanomaterials, or contrast agents [48][49][50] in the whole body of living mice for drug delivery monitoring [51], therapeutic assessment [52], tumor imaging [53], gastrointestinal angiography [54], lymphography [55], and cystography [56].
For visualization purposes, the mechanical scanning-based PAMac usually suppresses down its three-dimensional (3D) imaging capability to two-dimensional (2D) projected images. The main reason for displaying 2D images is to avoid depth-related distortion that occurs due to respiratory movement during relatively long (typically 30-60 min for whole-body imaging [47,[57][58][59]) data acquisition time, which is due to the physical traveling time of a single-element US transducer in the mechanical scanning. For better 3D analysis and visualization, the respiratoryrelated motion artifacts along the depth direction should be compensated.
To compensate the motion distortion, several studies have been reported for mechanical scanning-based systems. Schwarz et al. compensated motion by detecting skin signals, which were generated by the melanin-containing layers in the skin [60]. Then, they corrected the distortions by smoothing the surface profile. Zhao, et al. proposed a respiratory compensation algorithm, which tracks enhanced vascular structures in high-resolution PA images [61]. The algorithm was successfully applied to rat iris imaging and mice skin imaging. Those methods showed feasibility for motion correction during the mechanical scan, but they require continuous PA signals for motion tracking. Since the PA signals are not generated from the optically transparent tissues, a method that does not rely on PA signals is required for robust motion compensation.
In this paper, the compensation of breathing distortion in 3D whole-body PA data of mice, acquired from the mechanical scanningbased PA and US macroscopy (PAUSMac) system, is demonstrated.
The key approach for compensating the motion is the detection of the skin layer in the acquired images, which is often not apparent in PA images. If the skin layer can be detected using US images, where it is easy to identify and motion-compensated, the corresponding compensation can be applied to corresponding PA images to remove motion artifacts. To accomplish this, the recently reported dual-modal PAUSMac system was proposed to be used for wholebody imaging of small animals [62]. This system simultaneously acquires both PA and US data without increasing data acquisition time. The results also showed the potential of the dual-modal imaging system for reinforcing structural information, which can be achieved in USI, and ported onto PAI.
In the proposed method, the skin layers were segmented by using cross-sectional US B-mode data, then used for the distortion compensation in the depth direction. The corresponding distortions in PA data were corrected to achieve 3D images. The resulting images from the compensated data showed better depth-resolving quality compared to the previous images. The 3D distributions of hemoglobin oxygen saturation level (sO 2 ) were also acquired by spectral unmixing of the multispectral data with correction of respiratory distortion. This technique will benefit biomedical research by providing more accurate information on both the anatomy and physiology of the underlying tissue, thus improving the outcome of the in vivo small animal studies.

Photoacoustic and Ultrasound Data Acquisition
The previously developed dual-modal PA and US imaging system, which can acquire both PA and US data simultaneously, was used for this study ( Figure 1). In brief, Nd:YAG pumped Q-switch laser (Surelite III-10, Continuum, United States) generated a laser beam, whose wavelength was tuned by the optical parametric oscillator (OPO) laser system (Surelite OPO PLUS, Continuum, United States). The optical components delivered the laser beam to the target animals, which were placed underneath the bottom-opened water tank. The generated PA waves were detected by a single-element US transducer (V308, Olympus NDT, United States), captured by a pulser/receiver (5072PR, Olympus NDT, United States), stored by a data acquisition module (MSO 5204, Tektronix, United States), and reconstructed as images using Matlab-based custom software. The conventional US gel and membrane-sealed water tank were used to match the acoustic impedance between the US transducer and the target. The US data were acquired between each PA data acquisition. More details of the system can be found in our previous publication [62]. This study used multispectral lasers of wavelengths 750, 800, and 850 nm, with average fluences of 9, 7, and 7 mJ/cm 2 , respectively. A 60 × 32 mm 2 scanning was performed to achieve whole-body images of mice with a step size of 0.2 mm in both x and y directions. The total data acquisition time was~100 min at each wavelength.

Ultrasound-Guided Breath-Compensation
A post-processing algorithm that automatically detects the distortion from respiratory movement during the scanning and compensates the distortion along the axial direction was developed ( Figure 2A). Prior to the compensation, the conventional signal processing procedures (i.e., frequency demodulation, envelope detection, and interpolation) were applied to each A-line signal for generating images (Supplementary Figure S1). The breath-compensation algorithm was then applied to each cross-sectional B-mode data. The first step in the breath-compensation algorithm was the segmentation of the skin layer. From the US B-mode images, it was easy to segment the membranes (blue lines in Figure 2A) and skin layers (yellow lines in Figure 2A) using a simple thresholding method. After removing the membrane signal, a low-pass filter in the spatial domain was applied to smoothen skin layers (i.e., compensated skin layers, green lines in Figure 2A). Next, at each lateral point, the difference between the original and compensated skin surface was measured and applied to shift the original data in the axial direction. Through this shifting, the spiking distortions caused by respiratory movement (yellow arrows in Figure 2A) were compensated (green arrows in Figure 2B). The same shifting intervals were also applied to the corresponding PA data to compensate for breathing. Finally, breath-compensation was achieved for 3D PA and US data after repeating this procedure for every cross-sectional B-mode data ( Figure 2B) using Matlab-based open-access software (3D PHOVIS, POSTECH, Republic of Korea) [63]. The accuracy of the resulting 3D images was verified by comparing them with anatomical positions from commercial software (3D Rat Anatomy, Biosphera, Brazil).

Three-Dimensional Distribution of Hemoglobin Oxygen Saturation Level
For obtaining the 3D whole-body sO 2 distribution in mice, multispectral PA and US data were acquired using excitation wavelengths of 750, 800, and 850 nm, which are deoxyhemoglobin dominant, isosbestic, and oxy-hemoglobin dominant points, respectively [64]. The same procedure described in Figure 2 was applied to the acquired volumetric data. The breath-compensated skin profiles were achieved at each wavelength. For accurate sO 2 , we calculated the intervals between the reference profile to each skin profile, shifted the data to match the axial position, then spectrally unmixed oxy-and deoxy-hemoglobin components at each pixel of volumetric data. The breathcompensated skin profile at 750 nm was used as the reference profile. The relative sO 2 level at each pixel was calculated by the conventional linear unmixing method [65,66]. After smoothing the acquired volume by a 3D 3 × 3×15 median filter, which corresponds to 0.6 × 0.6 × 0.4 mm 3 in x, y, and z directions, respectively, 3D distributions of sO 2 level in the whole-body mice were calculated.

Animal Preparation
All animal experiments followed the protocol approved by the institutional animal care and use committee (IACUC) of the Pohang University of Science and Technology (POSTECH-2021-0052). Healthy 6-week-old female Balb/c nude mice (~20 g) were used to acquire PA and US data. The fine hairs that generate distortion in data were removed using a conventional hair removal cream (Nair, Church & Dwight, United States). During the data acquisition, the mice were anesthetized using a vaporized isoflurane gas system with a mixture of 1 L/min oxygen and 0.75% isoflurane. The mice were placed under a membrane-sealed water tank. The acoustic impedance between the membrane and mice skin was coupled using a conventional ultrasound gel (Ecosonic, SANIPIA, South Korea). After in vivo imaging at the left sagittal view, the mice were sacrificed under anesthesia by using an overdose of carbon dioxide gas.

Breath-Compensated Whole-Body Images
After US-guided breath-compensation was applied to the 3D volume data, the rendered 3D images of both US and PA were generated (Supplementary Movies S1, S2). In the 3D rendered images, distortions are significantly reduced, and the accuracy of the depth-information is highly improved. In the conventional PAMac images, the volumetric data are projected onto the x-y plane using the maximum amplitude projection (MAP) method along with the whole range of depth of data to hide the respiratory distortions. The combined mean intensity projection (MIP) image of the US and MAP image of PA can visualize coexisting internal organs (e.g., intestine, spleen, and liver) and corresponding major blood vessels that match well with the anatomy (Figures 3A-C). However, the MAP or MIP projected images lose depth information of the signals.
Using the breath-compensation technique developed here, which utilizes the position of the skin surface to compensate for the motion, we can overcome this issue and can calculate the accurate depth of each signal. The sliced images can be obtained by sectioning the data at 1-mm intervals underneath the skin surface. In the sliced US MIP images, structural information is clearly visible at each depth ( Figure 3D). Signals from the skin are dominant at the surface, but the internal anatomy and pathology can be visualized at subsequent deeper slices. More specifically, the rib cage, spine, and thigh are identifiable due to the strong reflection of US waves arising from the impedance mismatch between bone and soft tissues. The lung surface also generates strong signals because of impedance mismatch between the lung surface and air. The continuous changes in the sectioned US MIP images from the surface to the depth of 5 mm are shown in the Supplementary Material (Supplementary Movie S3).
The sliced PA MAP images show the optical absorption at each level ( Figure 3E) of the volume along with the depth. The colormap used for displaying the blood flow in each image has the same dynamic range. Since the hemoglobin strongly absorbs light, the blood vessels are generally dominant in PA images. As the depth increases, optically absorbing organs, such as the spleen, intestine, and liver, also become visible. Since the PA waves propagate in all directions, the signals from the blood vessels are often reflected by the lung and appear as artifacts in the deep position (green arrows in Figure 3E). The continuous change in the sectioned PA MAP images is also demonstrated in the Supplementary Movie S4.
In addition to the x-y plane, we also successfully achieved PA MAP images in the x-z and y-z planes ( Figure 4A) after compensating the distortions. Compared to the conventional landscape visualization of the x-y plane, depth information can be achieved in the other planes. In the x-z and y-z MAP images, the breath-compensation can be noticed after the removal of the respiratory distortion (white arrows in Figures 4B,C) in the compensated images. The cross-sectional B-mode images clearly show the distortion correction in both US and PA images ( Figure 4D, the corresponding position in the 3D image is depicted as the green box in Figure 4A). The spiking respiratory motions (white arrows in Figure 4D) are successfully compensated. The results demonstrate the feasibility and effectiveness of our proposed US-guided breath-compensation technique for improving the quality of 3D information, especially for depth-related analysis.

Three-Dimensional Distribution of Hemoglobin Oxygen Saturation Level
Similar to other optical imaging techniques, PAI can also calculate the relative sO 2 level in biological tissue by unmixing multispectral data. In the conventional PAMac, sO 2 calculation is commonly performed using MAP images that are acquired at multiple excitation wavelengths ( Figure 5A). However, since the maximum signals may not originate from the same depth, we cannot assume that the multispectral data for spectral unmixing arise from the same position. Thus, in conventional techniques, accurate sO 2 calculation is limited due to the loss of depth information in the MAP images ( Figure 5B).
In this study, through breath-compensation at each excitation wavelength, we successfully obtained 3D sO 2 levels in the whole body of mice in vivo. Compared to the conventional PAMac studies, we performed spectral unmixing of oxy-and deoxyhemoglobins before the MAP procedure ( Figure 5C). We could then preserve the depth information of each data and match the axial position by shifting signals to the reference skin profile (described in the Method section). We achieved 3D distribution of whole-body sO 2 level at various view angles ( Figure 5D). The resulting 3D sO 2 distribution shows the depth-wise functional information, which otherwise cannot be achieved using conventional techniques. The 3D visualization of relative sO 2 distribution is also made available in the Supplementary Material (Supplementary Movie S5).

CONCLUSION
The mechanical scanning-based PAMac has been used in biomedical studies for analyzing the whole body of small animals in vivo. However, the mechanical scanning leads to a long data acquisition time, thus resulting in respiratory distortion in the depth direction of the acquired volume. With the distortions, it is difficult to visualize the 3D volumetric data in conventional PAMac. As an alternative, landscape MAP images have typically been used for the visualization of signal distribution. Although landscape MAP images are helpful in analyzing functional information distributed in the whole body of mice, they do not provide positions of signals in the depth direction, which are essential for accurate signal processing and visualization.
Here, we have successfully demonstrated a US-guided breathcompensation technique for the 3D visualization of mechanical scanning-based PAMac. We acquired PA and US data from the whole body of mice in vivo by using the previously developed dual-modal imaging system. Using the structural information in the cross-sectional US images, we successfully corrected respiratory distortions in the axial direction for both US and PA images. The successful 3D visualizations of the volumetric data showed the great potential of this technique for investigating the 3D distribution of optical contrast in vivo. The 3D wholebody sO 2 level was also calculated at each pixel after compensating respiratory movements. To the best of our knowledge, this is the first report that presents 3D whole-body sO 2 levels using PAMac.
The key procedure in the proposed breath-compensation method is detecting skin profiles with US data. The US signals are always generated from the skin surface due to the mismatching of acoustic impedance between skin and water medium. Therefore, the proposed method is not limited depending on the presence or absence of a PA signal. We expect that the US-guided breath-compensation method presented here will open new opportunities to visualize the comprehensive distribution of optical chromophores in various biomedical studies, including drug delivery monitoring, clearance assessment, lymphatic network mapping, and blood vessel imaging. In the near future, we plan to expand the use of the technique into many other areas as listed above.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of the Pohang University of Science and Technology.