Article

- Open Access

- Published:

High-speed three-dimensional photoacoustic computed tomography for preclinical research and clinical translation

- Li Lin,

- Peng Hu,

- Xin Tong,

- Shuai Na,

- Rui Cao,

- Xiaoyun Yuan,

- David C. Garrett,

- Junhui Shi,

- Konstantin Maslov &

- Lihong V. Wang

Nature Communications

volume

12, Article number: 882 (2021) Cite this article

- 8978 Accesses

- 44 Citations

- 1 Altmetric

- Metrics

details

Abstract

Photoacoustic computed tomography (PACT) has generated increasing interest for uses in preclinical research and clinical translation. However, the imaging depth, speed, and quality of existing PACT systems have previously limited the potential applications of this technology. To overcome these issues, we developed a three-dimensional photoacoustic computed tomography (3D-PACT) system that features large imaging depth, scalable field of view with isotropic spatial resolution, high imaging speed, and superior image quality. 3D-PACT allows for multipurpose imaging to reveal detailed angiographic information in biological tissues ranging from the rodent brain to the human breast. In the rat brain, we visualize whole brain vasculatures and hemodynamics. In the human breast, an in vivo imaging depth of 4 cm is achieved by scanning the breast within a single breath hold of 10 s. Here, we introduce the 3D-PACT system to provide a unique tool for preclinical research and an appealing prototype for clinical translation.

Introduction

Photoacoustic (PA) computed tomography (PACT) is a noninvasive hybrid imaging modality that combines the functional optical contrast of diffuse optical tomography with the high spatial resolution of ultrasonography1. When a short-pulsed laser irradiates biological tissues, the light energy is partially absorbed by the tissue and converted into heat. Subsequently, a pressure rise induced by transient thermoelastic expansion propagates as a wideband ultrasonic wave (referred to as PA wave) through the biological tissue. The PA waves are then detected by ultrasonic transducers placed around the tissue surface to reconstruct the optical absorption distribution in the tissue2. The sensitivity to optical absorption in biological tissues provides PACT-rich contrast mechanisms relating to various intrinsic and extrinsic chromophores, enabling structural, functional, and molecular imaging3. In addition to the rich optical contrast, the conversion from light to ultrasound also provides this imaging technique with ultrasonically defined spatial resolution at depths beyond the ~1-mm optical diffusion regime, which limits the penetration capability of ballistic optical imaging modalities4.

This hybrid nature grants PACT unique advantages and makes it complementary to other mainstream imaging modalities in clinical practice: (1) compared with pure optical imaging techniques such as fluorescence imaging in humans5, PACT can sustain high spatial resolution in the deep tissue. (2) compared with ultrasonic imaging, PACT has rich intrinsic and extrinsic optical contrasts and is free of speckle artefacts4; (3) compared with X-ray computed tomography (X-ray CT) and positron emission tomography, PACT is free of radioactive material and uses nonionizing illumination; (4) compared with magnetic resonance imaging (MRI), laser-based PACT provides higher spatial resolution at a lower cost, though it faces a challenge to penetrate the entire human body. Specifically, the cutting-edge PACT system we are presenting in this manuscript can image the human breast in 10 s with full penetration (4 cm) and nearly isotropic spatial resolution (0.37–0.39 mm) that is finer than breast 3T-MRI (0.5 mm in-plane, 1.3 mm thick)6 or a typical 7T-functional MRI (fMRI) (isotropic 0.75 mm)7.

Employing different light illumination and detection schemes, several PACT systems with three-dimensional (3D) field-of-view (FOV) have been developed for preclinical studies8,9,10,11 and clinical practice12,13,14,15,16,17,18. Although these systems have advanced PACT performance, key limitations remain unaddressed. Specifically, the current systems’ limitations mainly arise from their shallow imaging depth, slow imaging speed, and/or compromised image quality owing to the limited noise-equivalent sensitivity19 and FOV. Here, we define high image quality as sufficient contrast-to-noise ratio (CNR) to reveal detailed structures on the order of the inherent spatial resolution (i.e., high clarity) within an FOV large and deep enough to cover the target region. A useful rule of thumb is that a CNR of 2.2 can be detected with a confidence level of 98%20.

To provide a large view aperture (i.e., view angle)21, some researchers utilized spherical detection matrices, which, however, were designed and implemented differently8,9,10,12,13,14,15,16. For example, Deán-Ben et al.9 recently reported a PACT system with 256 transducer elements integrated on a partial cap for small animal imaging. This design can generate a small volumetric image with a single laser shot, but its well-resolved FOV is theoretically limited to a few millimeters (e.g., diameter <4 mm) according to the spatial Nyquist sampling criterion22. Subsequently, the authors scanned the system helically around a mouse body to densely sample the surface tissue and generated an elegant superficial image. However, without scanning along the radial direction, the mouse’s internal organs were not imaged clearly. The limited-view aperture (i.e., <2π steradian solid angle) further compromised the image quality when a higher imaging speed was required for brain imaging8.

For human breast imaging, Matsumoto et al.12 developed a PACT system using a sparse hemispherical detector array that scanned in a spiral pattern on a plane. The dense sampling produced high-quality vascular images near the skin surface. Similar to the aforementioned system, however, the well-resolved FOV of the detector array could only recover a small breast volume and the system could not image the deep tissue well without an elevational scan. In addition, the illumination wavelengths they used suffer from higher scattering in the breast and a lower safety limit for maximum permissible exposure (1/5 of the limit for 1064-nm light)23. Accordingly, the imaging depth of the system was limited to the superficial region of the breast. More recent work published by Oraevsky et al.13 and Schoustra et al.14 employ scanned arc ultrasonic arrays for human breast imaging. Both studies showed neither cross-sectional breast images at different depths nor 3D-rendered views, presumably owing to the limited imaging depth. Similarly, the sparse spatial sampling limited the well-resolved FOV in breasts, thus compromising the image quality while also requiring a long scanning time of several minutes. In addition, the scanning of the laser beam varied the light energy distribution,14 which violates the assumption of 3D image reconstruction algorithms24,25 that light energy distribution remains consistent during scanning. Other than the systems that utilized spherical detection matrices, researchers also explored 3D imaging by detecting with a planar sensor17 or scanning a linear array18. Though the systems were simpler, the spatial resolution was anisotropic and the image quality was reduced owing to the limited-view aperture.

In this paper, we developed a PACT system, which produces volumetric images of biological tissues with superior imaging depth, image quality, and speed. Rather than being dedicatedly used for a single application, we have demonstrated its in vivo imaging versatility in the rat brain and human breast, representing two extreme imaging scales. In the rat brain, we clearly imaged angiographic anatomies and functions from the cortex to the Circle of Willis. In the human breast, the system can generate a volumetric image with a penetration up to 4 cm by scanning the breast within a single breath-hold of 10 s. Our in vivo imaging results distinguish 3D-PACT as a cutting-edge PA imaging system that meets all the following conditions: (1) deep penetration and scalable FOV to accommodate imaging objects from the rat brain to human breast, (2) isotropically high spatial resolution within the FOV in 3D space, (3) high noise-equivalent sensitivity to reveal small structures deep in tissues, (4) high imaging speed to minimize motion artefacts and enable functional studies, and (5) minimal limited-view artefacts. Accordingly, 3D-PACT could provide high-quality images with functional optical contrast and high imaging speed to directly benefit both preclinical research and clinical practice, which have been highlighted by increasing interest from clinical researchers and practitioners.

Results

Characteristics of the 3D-PACT

The 3D-PACT system is shown in Fig. 1a. An imaging platform is placed above the ultrasonic array housing to support the imaging subject. A disposable holding cup (plastic wrap) is clamped on the imaging aperture to separate the imaging system from the subject. The 3D-PACT system is composed of five modules: (1) four arc-shaped ultrasonic transducer arrays are integrated in a hemispherical housing (Fig. 1b, c). Each array has 256 transducer elements with a central frequency of 2.25 MHz and a one-way bandwidth of more than 98% (Supplementary Fig. 1). (2) 1024-channel pre-amplification circuits are attached to the base of the ultrasonic array housing (Fig. 1b, Supplementary Fig. 2). Four phase plates transmit amplified PA signals to four data acquisition (DAQ) modules. (3) Each DAQ module has 256-channel analog-to-digital converters with a tunable gain from 12 to 51 dB. (4) A mechanical scanning module rotates the aforementioned three modules coaxially for 90 degrees to form a densely-sampled hemispherical detection matrix. (5) The light beam is directed to the subject from an Nd:YAG laser through a one-inch lens tube glued on the bottom of the hemispherical housing (Fig. 1b, c). Two engineered diffusers are sealed in the lens tube to expand the beam for illumination with redundancy for safety (Methods and Supplementary Fig. 3); the laser beam would be expanded even if one of the diffusers was broken or misaligned during imaging. We take advantage of the low optical attenuation of 1064 nm light to achieve deep optical penetration in biological tissues26.

Fig. 1: Representations and photographs of the 3D-PACT system and the leaf skeleton images.

figure 1

a Perspective view of the system. DAQ, data acquisition module. b Cut-away view of the system with the imaging platform removed. c Photograph of the system with a close-up view of the ultrasonic arrays. d Photograph and maximum amplitude projection (MAP) images of a leaf skeleton reconstructed after geometric calibration (see Methods).

Full size image

Owing to manufacturing defects in the transducer arrays, which are slightly imperfect circular arcs, the reconstructed images of a leaf skeleton phantom are slightly blurred (Supplementary Fig. 4). To correct these defects, we measured a point source placed at multiple positions and applied Jacobi iteration to find the true geometric position of each transducer element (Methods and Supplementary Fig. 4). An improvement in image quality was achieved and the detailed leaf skeletons were clearly revealed (Fig. 1d). After calibrating the imaging system, we further imaged a 50-µm-diameter aluminum particle, allowing us to experimentally quantify the spatial resolution in 3D space (390 µm on X–Y plane, 370 µm along Z direction) (Supplementary Fig. 5). The scanning time to generate a volumetric image ranges from 2 to 10 s depending on the targeted FOV22.

3D-PACT of the rat brain anatomy and functions

Small animals are principal models for preclinical studies, and the imaging of their brains has an important role in neuroscience research. 3D-PACT is capable of providing high-quality images of the whole rat brain anatomy within 5 s or can continuously monitor the cerebral functions with a volumetric imaging rate of 0.5 Hz, which can be further increased by utilizing a laser with higher repetition rate (e.g., 5 Hz for a 100-Hz laser).

As shown in Supplementary Fig. 6, the rat head was mounted horizontally with the cortical region attached to the center of the membrane. Although the top part of the skull was thinned (see Methods) to acquire detailed angiograms of the brain by reducing the attenuation and distortion of the skull on the PA waves, the skull remained intact for functional imaging. The 2.25-MHz center frequency allows 3D-PACT to produce high-quality images with tolerable acoustic distortion from the rat skull27. Taking advantage of the deep penetration of 1064 nm light and the panoramic acoustic detection, a volumetric image of the rat whole brain (10 mm in depth) was acquired with detailed vasculature (Fig. 2a and Supplementary Movie 1). The image was acquired by averaging eight 5-second scans and subsequently reconstructed using a dual-speed-of sound 3D back-projection algorithm (see Methods). Fig. 2b, c show cross-sectional images of the brain on sagittal and coronal planes, respectively. Supplementary Fig. 7 shows the brain images with comparable clarity acquired in 5 s without averaging.

Fig. 2: 3D-PACT of a rat brain in vivo.

figure 2

a Perspective angiograms of the rat whole brain (Supplementary Movie 1). b Cross-sectional images of the brain on different transverse planes. c Cross-sectional images of the brain on different coronal planes. AACA azygos of the anterior cerebral artery, ACH anterior choroidal artery/vein, DTV dorsal thalamic vein, ICA internal carotid artery, ISS inferior sagittal sinus, LE left eye, MCA middle cerebral artery, MOA medial orbitofrontal artery, OA olfactory artery, PCA posterior cerebral artery, RE right eye, RRV rostral rhinal vein, SS sigmoid sinus, SSS superior sagittal sinus, SV supraorbital vein, TH transverse hippocampal artery/vein; THA thalamoperforating artery, TS transverse sinus.

Full size image

Serving as the center of the nervous system, the brain coordinates neural activities which, as in fMRI, can be related to cerebral hemodynamics28,29. Imaging at a 0.5-Hz volumetric frame rate, we used 3D-PACT to monitor the hemodynamics while modulating inhaled oxygen (Fig. 3a–c) and anesthesia content (Fig. 3d–f). Specifically, we measured the hemodynamics in four major vessels (Fig. 3c, f, and Supplementary Fig. 8) during the two processes to demonstrate the advantages of 3D-PACT’s spatiotemporal resolution. As oxyhemoglobin dominates the optical absorption at 1064 nm (i.e., μa(oxyhemoglobin)≈10 × μa(deoxyhemoglobin), μa: absorption coefficient)26, the changes in spatially averaged PA signals mainly reflect the variation of oxyhemoglobin content within the region of interest (ROI) over time. For instance, to measure hemodynamics in major blood vessels, the ROIs were selected slightly larger than the vessels’ cross-sections, thus both the oxygen saturation (sO2) and the blood volume (or vascular diameter) would affect PA signal amplitudes. Furthermore, we also measured the intrinsic functional connectivity (Fig. 3g) and the hemodynamics induced by the electrical stimulations to the front limbs (Fig. 3h).

Fig. 3: Functional rat brain imaging.

figure 3

a PA amplitude changes in different brain cross-sections after the rat inhaling pure nitrogen (N2) for 30 s. b PA amplitude changes in the circled regions in a. c PA amplitude changes in four major arteries and veins during hypoxic challenges. AACA azygos of the anterior cerebral artery, OA olfactory artery, SSS superior sagittal sinus, ISS inferior sagittal sinus. d PA amplitude changes in different brain cross-sections while awakening from deep anesthesia. e PA amplitude changes in the circled regions in d. f PA amplitude changes in the same four major vessels while awakening from deep anesthesia. g Correlation matrix of 20 functional regions across the whole brain. Notice the correlation between left and right hemispheres. h Brain cortical hemodynamics in response to the electrical stimulation to the front limbs. CP caudoputamen, GU gustatory area, HIP hippocampal region, MOp primary motor area, PRN pontine reticular nucleus, ORB orbital area, RSP retrosplenial area, SC superior colliculus, TH thalamus, VISp primary visual area.

Full size image

To observe the brain’s hemodynamic response to hypoxic challenges, we manipulated the oxygen concentration in the inhalation gas (see Methods) to during imaging. In Fig. 3a, we selected three cross-sections on the coronal plane from bregma to lambda and observed a global decrease in the PA signal in response to 30 s of nitrogen inhalation. Fig. 3b shows the spatially averaged PA amplitude changes in the circled regions in Fig. 3a.

During hypoxic challenges, we noticed that the PA signal variations in the four vessels were comparable (Fig. 3c and Supplementary Fig. 8a)30. We also observed larger fractional changes in PA signals from the four major vessels compared with brain tissue (Fig. 3b), which includes both blood vessels and brain cells. To explain this difference, we surmise that in Fig. 3b, the time-invariant PA signals from lipid, water, and proteins31 moderated the change in PA signals from hemoglobin. Another observation is that the slopes of the PA signals plateaued during the 40-second hypoxia in Fig. 3b, c, which is likely owing to a reduced oxygen consumption rate in response to extended hypoxia. In comparison, the slopes of the PA signals are more consistent during the 10-second hypoxia. This agrees with the expectation that oxygen consumption rate is more stable at the beginning of hypoxia32.

In addition to hypoxic challenges, we measured PA signal changes from deep to light anesthesia (Fig. 3d–f). Isoflurane, a well-known vasodilator33, can also reduce the cerebral oxygen extraction fraction34,35. Therefore, the decrease in PA signals within different brain regions from deep to light anesthesia (Fig. 3d, e) could be a combined effect of reduced blood volume and sO2.

Moreover, we also monitored hemodynamics in the same major vessels during the awakening process (Fig. 3f and Supplementary Fig. 8b). Similar hemodynamics were observed in the two arteries and two veins. To interpret this similarity, we surmise that in veins, both the decrease in sO2 and the contraction of vascular diameter35 contribute to the PA signals’ decrease. In arteries, although sO2 is expected to be more stable during the awakening process35, the vascular contraction may be more pronounced than that in veins35. Another observation is that the decrease in PA signals from major vessels (Fig. 3f) was comparable to that in brain tissues (Fig. 3e). To explain this, our hypothesis is that in Fig. 3e, although the time-invariant PA signals from lipid, water, and proteins could moderate the changes in PA signals from hemoglobin, small vessels in the brain tissue are expected to contract more than major vessels during the awakening process35.

Like fMRI36, we used 3D-PACT to assess the intrinsic functional connectivity across spatially separated brain regions through spontaneous fluctuations in PA signals (particularly during resting state) by measuring the spontaneous hemodynamic responses in cerebral regions (Supplementary Movie 2). We identified 20 functional regions across the whole brain and computed the correlation coefficients of every pair (see Methods). Fig. 3g shows clear correlation between the contralateral regions in left and right brain hemispheres. The high volumetric imaging rate also allows us to measure faster cerebral hemodynamics (Supplementary Movie 2). In Fig. 3h, we observed higher PA amplitude changes in the corresponding motor sensory regions when the rat’s front limbs were separately stimulated by periodic electrical pulses (see Methods). We further plotted the PA amplitude changes in the motor sensory regions in response to the electrical stimulation.

Our anatomical and functional brain imaging results demonstrate the potential of 3D-PACT as an imaging tool with deep penetration, isotropic resolution, and high speed for studying rat whole brain structures and functions, which were previously challenging to accomplish using optical contrast. For small animal functional imaging, whereas the imaging speed of 3D-PACT (0.5-Hz volumetric imaging rate) is not yet substantially faster than fMRI (9–11 slices in 2 s)37, the high spatial resolution of 3D-PACT enabled measurement of hemodynamics in individual vessels. Although the rat brain can be fixed in place during imaging using a nose cone and tooth bar, artefacts induced by periodic motions (e.g., breathing) can be largely removed by time-gated motion correction (see Methods). We have demonstrated this technique by imaging the abdomen of a pregnant rat, whose embryo was clearly revealed after correcting the respiratory motion (Supplementary Fig. 9). Time gating can also be used to reconstruct fast but periodic movements (e.g., heart beating). Accordingly, in addition to brain imaging, 3D-PACT can be used to reveal structural, functional, and molecular information in many other parts of small animals.

3D-PACT of human breasts in vivo

Multiple large prospective clinical trials have demonstrated the importance of early detection in improving breast cancer survival38,39. In addition, more recent studies have shown that tumor biology also has an important role in the prognosis of an individual cancer40,41. An imaging modality with high sensitivity and specificity which is cost effective and safe to use would provide an important unmet medical need42,43.

For breast imaging, PACT can offer high spatial and temporal resolutions with sufficiently deep nonionizing optical penetration44. In biological tissue, the optical absorption contrast of PACT is much higher than mammographic X-ray contrast45. However, the 1/e attenuation coefficient for 1064-nm light in an average breast is ~0.9 cm−1 26,46, which is slightly higher than that for X-rays (0.5–0.8 cm−1)47. In the near-infrared (NIR) region, hemoglobin is the principal optical absorber to provide PACT an endogenous contrast (i.e., blood vessel density and morphology, total hemoglobin concentration, oxygenation state, etc.), which is closely correlated with the tumor growth and metastasis48,49,50. The optical absorption contrast can also serve as an indicator of the breast cancer’s response to treatment51,52.

To demonstrate the 3D-PACT of breast imaging, we imaged the two breasts of a 31-year-old healthy female volunteer. Her breast size was 36 C, which is the average size in the United States53. The 3D-PACT system was placed underneath a dedicated imaging bed44 with minimal separation. With the human subject lying prone on the bed, the breast to be imaged was slightly compressed against the chest wall by the holding cup to reduce the optical tissue thickness from the nipple to the chest wall (Supplementary Fig. 6). By scanning the transducer arrays around the breast within a single breath hold of 10 s, we clearly revealed the blood vessels from the nipple to the chest wall (Fig. 4, Supplementary Movie 3 and 4). The cross-sectional images in Fig. 4b, d show detailed angiographic structures at different depths up to 4 cm, visualizing the vasculature down to an apparent vascular diameter of 400 µm. In our breast images, most large blood vessels are distributed near the skin surface (i.e., <2 cm deep, mainly derive from superior, lateral, and internal thoracic vessels), whereas blood vessels deeper in the breast (i.e., intercostal perforators) tend to be smaller and perpendicular to the chest wall. Such a distribution complies with the general angiographic anatomy of the human breast54,55. To further test the sensitivity of detecting small abnormalities in the breast, we imaged a breast-mimicking phantom with tumors of different sizes embedded at 2 cm depth (Supplementary Fig. 10). We made the tumor and breast phantoms with similar optical properties as in real human breasts56. The absorption coefficient of the tumors (μa = 0.105 cm−1) was 2.1 times that of the breast phantom (μa = 0.05 cm−1)56,57. The reduced scattering coefficient μ′s of the phantoms was 5 cm−1 at 1064 nm26,57. The system is sufficiently sensitive to detect a 1 mm-diameter tumor embedded in the breast-mimicking phantom.

Fig. 4: 3D-PACT of human breasts in vivo.

figure 4

a (top) Perspective angiogram of the right breast of a healthy human subject; (bottom) MAP image of the right breast view from the side. An imaging depth of 4 cm from the skin surface has been achieved (Supplementary Movie 3). b Cross-sectional images of the right breast on different coronal planes from the nipple to the chest wall. Each cross-sectional image is an MAP of a 1 cm-thick slice of the breast. c, d Left breast images of the same human subject (Supplementary Movie 4).

Full size image

Discussion

We have developed a 3D-PACT system that reveals detailed angiographic structures and functions with deep penetration, isotropic spatial resolution in scalable FOV, and high imaging speed. The massively parallel acoustic detection and one-to-one mapped low-noise amplifiers and DAQ circuits enable 3D imaging of an entire rat brain in 2 s and a human breast in 10 s. The high imaging speed allows 3D functional imaging of the whole rat brain, including the analysis of intrinsic functional connectivity and the measurement of hemodynamics while modulating inhaled oxygen, anesthesia content, and electrical stimulations. The large view aperture (i.e., ≈2π steradian solid angle) provides nearly isotropic spatial resolution and high image quality. In addition, by combining the 1064-nm light illumination with 2.25-MHz ultrasonic detection, 3D-PACT achieved up to 4 cm in vivo imaging depth in human breast, which is expected to be sufficient to image a D cup breast with a painless compression against the chest wall (over 99% of the US population have breast sizes of D cup or smaller53). Furthermore, 3D-PACT requires neither ionizing radiation nor an exogenous contrast agent, yielding zero risk to the imaging subject. The high imaging speed makes the 3D-PACT even more appealing.

The modular design also makes the 3D-PACT system easy to be maintained and upgraded. With the current hardware configuration, 3D-PACT has limitations in its single-wavelength illumination and mechanical scanning time (2–10 s). These limitations could be mitigated by replacing certain modules in the system: (1) using multi-wavelength illumination to provide additional imaging contrast58 (e.g., indocyanine green); (2) using lasers with higher repetition rates to image a larger FOV during the same scanning time or image the same FOV within shorter time; (3) replacing the stepper motor with a servomotor with higher torque to allow faster mechanical scanning; (4) using lasers with higher pulse energy to create a larger illumination area for breast imaging. For example, to improve the speed for small animal functional imaging, one could use a faster laser (e.g., LPY742-100, 100 Hz repetition rate, Litron Lasers Ltd.) and a servomotor with higher torque (e.g., EMG-50ASA, 3384 oz-in torque, Anaheim Automation Inc.). The system could then scan at 5 Hz for functional 3D imaging. One could also couple an OPO laser (e.g., SpitLight 1000 OPO, InnoLas Laser, Inc.) for multi-contrast imaging.

It is likely feasible to image awake rodents using 3D-PACT, but the system may need to employ ultrasonic arrays with smaller footprints (i.e., diameter) for easy operation. Inspired by PA microscopy of the awake mouse brain35, we could use a similar animal holder with the ultrasonic arrays placed above the animal. Considering the multiple applications this system could have, we designed the entire system to be rotatable (Supplementary Fig. 11) such that it could fit various scenarios with different imaging subjects and positions. The clinical values of the system are highlighted by the collaborations it has brought to us with the physicians from the City of Hope National Medical Center (breast cancer), Cedars-Sinai Medical Center (plastic and reconstructive surgery), Rancho Los Amigos National Rehabilitation Center (neurological surgery), and Ronald Reagan UCLA Medical Center (neonatology).

In this pilot study, we demonstrated the versatility and scalability of the 3D-PACT system that has a high reliability in reproducibly generating high-quality images. This system has advanced state-of-the-art PA tomography by providing a 3D imaging modality with superior performance for both preclinical research and clinical practice. To promote this technology, we share several design principles, which proved to be highly useful: (1) Since PACT reconstruction assumes light has a uniform or predicable distribution, light illumination needs to be carefully designed. High optical fluence within the safety limit is preferable for high signal-to-noise ratio; (2) a large view angle of the transducer array or detection aperture is helpful to minimize the loss of information and improve the image quality; (3) since PA signals generated from the deep tissue could suffer both light and acoustic attenuation, analog pre-amplifiers connected between the ultrasonic arrays and DAQs can amplify the weak PA signals before cable noise compromises the SNR; (4) in addition to the temporal Nyquist sampling law, one should also consider the spatial Nyquist sampling criterion22. Inadequate spatial sampling will result in limited FOV and undersampling artefacts; (5) high imaging speed is helpful to reduce the motion-induced artefacts. Although co-registration methods can partially mitigate the motion distortion, the non-rigidity of the biological tissue compromises the effectiveness of the co-registration; (6) proper grounding of the metallic housing of the transducer array, pre-amplifiers, and DAQs is important to reduce noise. The motors and their drivers are often well-shielded to eliminate electromagnetic interference; (7) finally, a reliable system with outstanding performance is always well-organized, where a concise design and neat arrangement facilitate problem debugging and reduce the likelihood of system failure.

Methods

System construction

The 3D-PACT system is mainly comprised of an ultrasonic detection module, pre-amplification and DAQ circuits, a rotation scanner, and an illumination laser. To achieve 3D panoramic ultrasonic detection, we employed four arc-shaped ultrasonic transducer arrays (Imasonic, Inc.; 2.25 MHz central frequency; 1024 elements) with a separation of 90 degrees that were mounted on a hemispherical surface (Fig. 1c). In each array, 256 wideband transducer elements were integrated on an 83-degree arc with a 130-mm radius. Each element had a flat-rectangular aperture (0.6 mm × 0.7 mm element size; 0.74 mm pitch). The hemispherical surface was constructed with white plastic such that back-scattered light from the tissue surface would be partially recycled and reflected to the tissue. The ultrasonic array housing was electrically grounded and was supported by eight metallic pillars that were fixed on an aluminum rotary board (Fig. 1b and Supplementary Fig. 2).

Above the rotary board, 32 sets of custom 32-channel pre-amplifiers (51 dB gain) were directly connected to the base of the ultrasonic arrays’ housing, amplifying the PA signals prior to experiencing cable noise. The amplified PA signals were then transmitted to four DAQ modules via four customized phase plates (PhotoSound, Inc.). Each DAQ module had 256-channel analog-to-digital converters (PhotoSound, Inc., ADC256; 40 MHz maximum sampling rate; 12 bit dynamic range), streaming the digitized data to a solid-state drive via USB 3.0. The DAQ modules were operated using National Instruments LabVIEW 2018 (64-bit).

To coaxially scan the ultrasonic transducer arrays for 90 degrees, we mounted two tapered-roller bearings with four-inch shaft diameter beneath the rotary board to provide coaxial alignment and steady rotation. The coaxial alignment was confirmed by aligning three one-inch holes on the ultrasonic array housing, rotary board, and driven spur gear with a one inch-diameter post. The rotation was motorized by a stepper motor (NEMA 34, 8 V) coupled with a set of two spur gears (Designatronics, Inc., KSS2-20J12 and KSS2-120, gear ratio = 1:6). The angular rotation speeds of the arrays were π4 rad/sec, π10 rad/sec, and π20 rad/sec for 0.5-Hz functional imaging, 5-second brain anatomy imaging, and single breath-hold breast imaging, respectively. According to the spatial Nyquist sampling criterion22, the FOV in brain anatomy imaging should have a diameter of ~2 cm, within which the spatial resolution is nearly isotropic (Supplementary Fig. 5). For rat abdomen and human breast imaging, we scanned the arrays for 10 s to provide a larger FOV (~4 cm). The lateral resolution outside the FOV was enlarged in relation to the increase in distance from the center (i.e., scanning axis)22. To generate a larger FOV within the same scanning time, one could use a faster laser (e.g., LPY7875-20, 20 Hz, Litron Lasers Ltd.) for illumination.

In 3D-PACT, an image reconstruction assumes a consistent light distribution in the imaging subject while scanning the ultrasonic arrays (i.e., spatial sampling)24. To achieve this, we guided the laser beam through the three coaxially aligned holes to a lens tube glued on the ultrasonic array housing. The lens tube was made of acrylic plastic to minimize the PA signal generated from it. The beam was further expanded by two engineered diffusers (EDC-40, EDC-15, RPC Photonics Inc.) sealed in the lens tube. The beam diameter was expanded to ~8 cm for human breast imaging (Supplementary Fig. 3). For rat brain and abdomen imaging, we added a condenser lens (ACL25416U-B, Thorlabs Inc.) after the diffusers (Supplementary Fig. 6a) to shrink the beam to a diameter of ~5 cm. As the light is scattered soon after propagating into biological tissues, we estimated that the energy distribution was relatively uniform across regions of the same depth. The optical fluence on the tissue surface was limited by the American National Standards Institutes safety standards23. Given the low optical attenuation in biological tissues, we utilized 1064-nm light for illumination, which is the fundamental wavelength of the Nd:YAG laser (PRO-350-10, Quanta-Ray, 10-Hz pulse repetition rate, 8–12-ns pulse width). To minimize the optical and acoustic energy loss, we used deuterium oxide (D2O, Isowater Corp.) to acoustically couple the holding cup (3–4 cm deep) with the ultrasonic arrays. The space between the biological tissue and holding cup was filled by a small amount of bathwater.

To synchronize the 3D-PACT system, the laser’s external trigger was used to trigger the DAQ modules, the rotation drive, and the linear power supply (IHE5-18/OVP, 5VDC, 18 A, International Power Corporation Ltd) of the pre-amplifiers (Supplementary Fig. 2). The power supply remained on during the experiment, whereas the power transmission to pre-amplifiers was controlled by a relay which functioned as a switch. Three seconds prior to the data acquisition, the relay was closed to transmit DC power to the pre-amplifiers. As the pre-amplification circuits were purely analog with the use of solid-state electronics, they amplified PA signals stably during data acquisition. Because the data acquisition time was limited, heat accumulation in the circuits was negligible. The four DAQ modules ran in parallel from a trigger supplied to the master DAQ module. The signal, trigger, and power cables were coiled around the rotation axis to reduce their required space. To prevent cable breaking, the system was programmed for round-trip scanning using a graphical user interface (GUI) generated by MATLAB 2018.

Geometric calibration of the ultrasonic arrays

Owing to manufacturing defects, the geometric position of each ultrasonic transducer element required calibration. We imaged a point source placed at multiple positions within the FOV and iteratively corrected the geometric location of each transducer element. To emit strong point source PA waves, a point absorber (~100 µm) was glued to an optical fiber tip such that most of the light coupled into the fiber was absorbed. Because the arrays’ housing was coaxially aligned and mounted, we assumed all the transducer elements share the same rotation axis.

The distance between the point source and each transducer element can be described by the following quadratic equations:

(xm−xn)2+(ym−yn)2+(zm−zn)2=(c⋅tm,n)2m=1,2,…,M;n=1,2,…,N

(1)

Here, M and N are the numbers of point sources and elements, respectively; (xm,ym,zm) denotes the location of the m-th point source; (xn,yn,zn) represents the location of the n-th transducer element; tm,n is the acoustic propagation time from the m-th point source to the n-th element; and c is the speed of sound in water. The center of the hemispherical detection matrix is chosen to be the origin.

We can measure tm,n precisely since the DAQ samples at 40 MHz. We have initial estimations of the unknown parameters c from its relationship with the water temperature, (xm,ym,zm) from the reconstructed images, and (xn,yn,zn) from the manufacture’s specifications. As our initial estimations should be close to the true values, Eq. (1) can be assumed to be locally linear and the Jacobi iteration59 was applied to find the true values of (xn,yn,zn). In each iteration, we updated the parameters as follows:

1. (1)

Speed of sound c:

c=1MN∑MM=1∑Mn=1(xm−xn)2+(ym−yn)2+(zm−zn)2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√/tm,n

(2)

2. (2)

Point source location (xn,yn,zn):

2(xn−xn−1)xm+2(yn−yn−1)ym+2(zn−zn−1)zm

=(c⋅tm,n−1)2−(c⋅tm,n)2+x2n+y2n+z2n−x2n−1−y2n−1−z2n−1

m=1,2,…,M;n=2,3,…,N

(3)

As M was much less than N, we applied the method of least squares method to solve the linear system and obtained (xm,ym,zm).

3. (3)

Element location (xn,yn,zn)

Because the ultrasonic arrays were manufactured by saw dicing60, the geometric displacement along the lateral direction is negligible. Accordingly, we only calibrated the geometric defects in the radial direction. We defined a radial calibration factor αn to confine the geometric defects along the radial axis and solved it from the quadratic equation:

∑Mm=1(xm−αnxn)2+(ym−αnyn)2+(zm−αnzn)2=∑Mm=1(c⋅tm,n)2

(4)

We iteratively updated the three sets of parameters until they converged (50 times iteration). Supplementary Fig. 4a shows the images of a leaf skeleton reconstructed before the geometric calibration. Supplementary Fig. 4b shows the radial calibration factor αn of each transducer element. The image quality is clearly improved after calibration (Fig. 1d).

3D image reconstruction

We used the universal back-projection (UBP) algorithm24 (Supplementary Methods) implemented in C++ to reconstruct all images in this work. The time-domain PA signals acquired at all scanning steps were back-projected into the 3D space using a GPU-accelerated implementation of the UBP. To mitigate aliasing artefacts caused by spatial undersampling in regions outside of the well-resolved FOV, we applied spatial interpolation and low-pass filtered PA signals with cutoff frequencies determined by the distance to the center of the imaging aperture22,61. A dual-speed-of-sound UBP algorithm61 was then used to mitigate the artefacts induced by the acoustic inhomogeneity between the biological tissue and the coupling fluid (D2O). We first reconstructed an initial image using the conventional UBP with a single speed sound24, then extracted the biological tissue’s surface (i.e., holding cup’s surface) from the initial image by moving detection windows towards the biological tissue and identifying the location of the first significant increase in PA amplitudes. The extracted surface was further fitted to a half ellipsoid. We then applied the dual-speed-of-sound UBP algorithm by assigning the speed of sound in tissue and D2O to regions inside and outside the ellipsoid, respectively.

Each volumetric image was reconstructed with a voxel size of 0.13 × 0.13 × 0.13 mm3. We further batch-processed all the reconstructed images to improve contrast. For image post-processing, we first applied a depth compensation (e0.81×depth(cm) for brain; e0.72×depth(cm). for breast) method to enhance the PA amplitude in the deep tissue. The compensated images were then denoised using sparse 4D transform-domain collaborative filtration62. We further applied Hessian-based Frangi vesselness filtration63 to the denoised images to enhance the contrast of blood vessels. Finally, we added the filtered images (self-normalized) with a weighting factor of 0.2 back to the original depth-compensated images (with a weighting factor of 0.8) and obtained the presented images. The 3D-rendered rat brain image (Supplementary Movie 1) was generated using VolView 3.4, and the 3D-rendered breast images (Supplementary Movies 3 and 4) were produced using the visvis library in Python 3.7.4.

Time-gated motion correction

Because the rat brain and human breast were well-fixed during imaging, motion correction was not applied to reconstruct the presented images. However, artefacts induced by periodic motions can be largely removed by time-gated motion correction (Supplementary Fig. 9). For example, when we scanned the abdominal region of a 3-month-old pregnant rat (Hsd:Sprague Dawley SD, Harlan Co.; Timed mated, 20 days) using 3D-PACT, respiration-induced motions appeared as spikes in the acquired PA signals (Supplementary Fig. 9a). To recover the signals acquired during motion, we detected the first-arrival times of the PA signals in different channels and applied a median filter with a kernel size of three channels to the first-arrival times. We then calculated the numerical difference of the signals’ first-arrival times. Signals with numerical differences >0.025 µs (temporal sampling step size) were regarded as motion-affected signals and were removed. The remaining signals formed a motionless partial-scan detection. Since the breathing cycles of the subject (i.e., rat abdomen) were spatially repeatable, we repeated the above procedures with four sequential 10-second scans and spliced the motionless partial-scan signals to form a full-scan detection (Supplementary Fig. 9a). The motionless full-scan signals were then used to reconstruct the image, in which the PA signals better converged (Supplementary Fig. 9b). Equipped with a faster laser, we could use this time-gated motion correction method to image faster phenomena (e.g., heartbeats).

Animal preparation

Adult, 1-month-old rats (Hsd:Sprague Dawley SD, Harlan Co.; 100–120 g body weight, male) were used for functional brain imaging. Before imaging, the rat was maintained under anesthesia with 1.5% vaporized isoflurane and its eyes were covered by pharmaceutical-grade ophthalmic ointment. The hair on the head was first removed using clippers and depilatory cream. We then opened the scalp, thinned the cranial surface with a handheld drill, and polished the surface with a green stone burr. For functional imaging, the skull was kept intact but the top part of the scalp was removed. Afterwards, we placed a 3D-printed nose cone and tooth bar on the rat and sealed the gap with silicone glue to prevent leakage.

The rat was then placed in a supine position and secured to a lab-made animal holder that was covered by a feedback-controlled heating pad to maintain the animal’s body temperature. To measure the brain’s response to hypoxic challenges, we switched the inhaled gas between air and pure nitrogen (Fig. 3b, c). To monitor the hemodynamics while awakening from deep anesthesia, we turned off the isoflurane vaporizer until the rat started to slightly move its limbs. During the measurement of intrinsic function connectivity and hemodynamics in response to the electrical stimulations, the isoflurane level was maintained at 0.6% so that the brain could be more active. The electrical stimulations were applied to the rat by binding two electrodes from an isolated pulse stimulator (A-M Systems, Inc., Model 2100) to one of the rat’s front limbs. The stimulator delivered 2-mA pulsed current to the electrodes with a pulse width of one millisecond. Each stimulation period lasted 12 s (4 Hz pulse frequency), followed by a 12-second rest (Fig. 3h).

Intrinsic functional connectivity of the rat brain

We adopted a data analysis method previously reported64 to measure the intrinsic functional connectivity of the rat brain. Each functional connectivity measurement took 8 min, with a 0.5 Hz volumetric imaging rate of the whole rat brain. After reconstructing all the 240 (0.5 × 60 × 8 = 240) time sequential images, a second-order Butterworth bandpass filter (0.008–0.09 Hz) was first applied to all temporal sequences, followed by a global signal regression based on the time sequences of all pixels within the brain. For functional connectivity analyses, we identified 20 functional regions throughout the brain, averaged the signals from pixels within each region, and computed correlation coefficients between regions to form a connectivity network matrix (Fig. 3g). Supplementary Movie 2 was made using a specialized 3D visualization software package, namely 3D PA Visualization Studio65.

Imaging protocols

All human and animal imaging experiments were performed in accordance with relevant guidelines and regulations. The human imaging experiment followed the protocol approved by the Institutional Review Board (IRB) of the California Institute of Technology. Human breast imaging was performed in a dedicated imaging room. After verbally agreeing to participate in the study, the subject signed the informed consent form, which was locked in a cabinet by researchers. The animal experiments followed the protocol approved by the Institutional Animal Care and Use Committee (IACUC) of the California Institute of Technology. IRB and IACUC were aware of both protocols before approving.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability

All data are available within the Article and Supplementary Files, or available from the corresponding author upon reasonable request. Imaging data from in vivo experiments and study approvals issued by the IRB and IACUC are available in the online folder: https://figshare.com/articles/dataset/3D-PACT_Data_and_Codes/13114544.

Code availability

The image reconstruction and processing algorithms are described in detail in Methods. The reconstruction code is not publicly available because the code is proprietary and is used in licensed technologies. Pseudo-code and description of the universal back-projection algorithm used in this study are provided in the Supplementary Methods. Imaging processing and analyzing programs are available in the corresponding data folders.

References

1. Wang, L. V. Prospects of photoacoustic tomography. Med. Phys. 35, 5758–5767 (2008).

Article PubMed PubMed Central Google Scholar

2. Zhou, Y., Yao, J. & Wang, L. V. Tutorial on photoacoustic tomography. J. Biomed. Opt. 21, 61007 (2016).

Article PubMed Google Scholar

3. Yao, J., Xia, J. & Wang, L. V. Multiscale functional and molecular photoacoustic tomography. Ultrason. Imaging 38, 44–62 (2016).

Article PubMed Google Scholar

4. Xia, J., Yao, J. & Wang, L. V. Photoacoustic tomography: principles and advances. Electromagn. Waves (Camb.) 147, 1–22 (2014).

Article Google Scholar

5. Marshall, M. V. et al. Near-infrared fluorescence imaging in humans with indocyanine green: a review and update. Open Surg. Oncol. J. 2, 12–25 (2010).

PubMed PubMed Central Google Scholar

6. Rahbar, H., Partridge, S. C., DeMartini, W. B., Thursten, B. & Lehman, C. D. Clinical and technical considerations for high quality breast mri at 3 tesla. J. Magn. Reson. Imaging 37, 778–790 (2013).

Article PubMed PubMed Central Google Scholar

7. Feinberg, D. A. & Yacoub, E. The rapid development of high speed, resolution and precision in fMRI. Neuroimage 62, 720–725 (2012).

Article PubMed Google Scholar

8. Gottschalk, S. et al. Rapid volumetric optoacoustic imaging of neural dynamics across the mouse brain. Nat. Biomed. Eng. 3, 392–401 (2019).

Article CAS PubMed PubMed Central Google Scholar

9. Deán-Ben, X. L., Fehm, T. F., Ford, S. J., Gottschalk, S. & Razansky, D. Spiral volumetric optoacoustic tomography visualizes multi-scale dynamics in mice. Light Sci. Appl. 6, e16247 (2017).

Article ADS PubMed PubMed Central CAS Google Scholar

10. Brecht, H. P. et al. Whole-body three-dimensional optoacoustic tomography system for small animals. J. Biomed. Opt. 14, 064007 (2009).

Article ADS PubMed PubMed Central Google Scholar

11. Jathoul, A. P. et al. Deep in vivo photoacoustic imaging of mammalian tissues using a tyrosinase-based genetic reporter. Nat. Photonics 9, 239 (2015).

Article ADS CAS Google Scholar

12. Matsumoto, Y. et al. Visualising peripheral arterioles and venules through high-resolution and large-area photoacoustic imaging. Sci. Rep. 8, 14930 (2018).

Article ADS PubMed PubMed Central CAS Google Scholar

13. Oraevsky, A. et al. Full-view 3D imaging system for functional and anatomical screening of the breast. Proc. SPIE 10494, Photons Plus Ultrasound: Imaging and Sensing (2018).

14. Schoustra, S. M. et al. Twente photoacoustic mammoscope 2: system overview and three-dimensional vascular network images in healthy breasts. J. Biomed. Opt. 24, 121909 (2019).

Article ADS PubMed Central Google Scholar

15. Kruger, R. A. et al. Dedicated 3d photoacoustic breast imaging. Med. Phys. 40, 113301 (2013).

Article PubMed PubMed Central Google Scholar

16. Dean-Ben, X. L. & Razansky, D. Functional optoacoustic human angiography with handheld video rate three dimensional scanner. Photoacoustics 1, 68–73 (2013).

Article PubMed PubMed Central Google Scholar

17. Heijblom, M., Steenbergen, W. & Manohar, S. Clinical photoacoustic breast imaging: the twente experience. IEEE Pulse 6, 42–46 (2015).

Article PubMed Google Scholar

18. Tan, Y., Xia, K., Ren, Q. & Li, C. Three-dimensional photoacoustic imaging via scanning a one dimensional linear unfocused ultrasound array. Opt. Express 25, 8022–8028 (2017).

Article ADS CAS PubMed Google Scholar

19. Winkler, A. M., Maslov, K. & Wang, L. V. Noise-equivalent sensitivity of photoacoustics. J. Biomed. Opt. 18, 097003 (2013).

Article PubMed PubMed Central Google Scholar

20. Contrast versus noise, eye versus machine. https://camera.hamamatsu.com/jp/en/technical_guides/contrast/index.html.

21. Xu, Y., Wang, L. V., Ambartsoumian, G. & Kuchment, P. Reconstructions in limited-view thermoacoustic tomography. Med. Phys. 31, 724–733 (2004).

Article PubMed Google Scholar

22. Hu, P., Li, L., Lin, L. & Wang, L. V. Spatiotemporal antialiasing in photoacoustic computed tomography. IEEE Trans. Med. Imaging vol. 39, 3535–3547 (2020).

23. American national standards institute. American National Standard for the Safe Use of Lasers ANSI, z136.131–2007 (Laser Institute of America, Orlando, FL, 2007).

24. Xu, M. & Wang, L. V. Universal back-projection algorithm for photoacoustic computed tomography. Phys. Rev. E 71, 016706 (2005).

Article ADS CAS Google Scholar

25. Pan, X., Zou, Y. & Anastasio, M. A. Data redundancy and reduced-scan reconstruction in reflectivity tomography. IEEE Trans. Image Process 12, 784–795 (2003).

Article ADS PubMed Google Scholar

26. Generic tissue optical properties. https://omlc.org/news/feb15/generic_optics/index.html.

27. O’Reilly, M. A., Muller, A. & Hynynen, K. Ultrasound insertion loss of rat parietal bone appears to be proportional to animal mass at submegahertz frequencies. Ultrasound Med. Biol. 37, 1930–1937 (2011).

Article PubMed PubMed Central Google Scholar

28. Chow, B. W. et al. Caveolae in cns arterioles mediate neurovascular coupling. Nature 579, 106–110 (2020).

Article ADS CAS PubMed PubMed Central Google Scholar

29. Lewis, S. A new kind of coupling. Nat. Rev. Neurosci. 21, 182 (2020).

PubMed Google Scholar

30. Choudhary, T. R. et al. Assessment of acute mild hypoxia on retinal oxygen saturation using snapshot retinal oximetry. Investig. Ophthalmol. Vis. Sci. 54, 7538–7543 (2013).

Article Google Scholar

31. Lei, L. et al. Label-free photoacoustic computed tomography of whole mouse brain structures ex vivo. Neurophotonics 3, 035001 (2016).

Google Scholar

32. Greenough, A. & Miller, A. D., Manual of Neonatal Respiratory Care. 2nd edn, 300–302 (2006).

33. Schwinn, D. A., Mcintyre, R. W. & Reves, J. G. Isoflurane-induced vasodilation: role of the alpha-adrenergic nervous system. Anesth. Analg. 71, 451–459 (1990).

Article CAS PubMed Google Scholar

34. Lyons, D. G., Parpaleix, A., Roche, M. & Charpak, S. Mapping oxygen concentration in the awake mouse brain. Elife 5, e12024 (2016).

35. Cao, R. et al. Functional and oxygen-metabolic photoacoustic microscopy of the awake mouse brain. Neuroimage 150, 77–87 (2017).

Article PubMed Google Scholar

36. Jonckers, E. & Van Audekerke, J. De Visscher, G., Van der Linden, A. & Verhoye, M. Functional connectivity fMRI of the rodent brain: comparison of functional connectivity networks in rat and mouse. PLoS ONE 6, e18876 (2011).

Article ADS CAS PubMed PubMed Central Google Scholar

37. Stenroos, P. et al. Awake rat brain functional magnetic resonance imaging using standard radio frequency coils and a 3d printed restraint kit. Front. Neurosci. 12, 548 (2018).

Article PubMed PubMed Central Google Scholar

38. Burton, R. & Bell, R. The global challenge of reducing breast cancer mortality. Oncologist 18, 1200 (2013).

Article PubMed PubMed Central Google Scholar

39. Kuhl, C. K. et al. Mammography, breast ultrasound, and magnetic resonance imaging for surveillance of women at high familial risk for breast cancer. J. Clin. Oncol. 23, 8469–8476 (2005).

Article PubMed Google Scholar

40. Hunter, C. P. Epidemiology, stage at diagnosis, and tumor biology of breast carcinoma in multiracial and multiethnic populations. Cancer 88, 1193–1202 (2000).

Article CAS PubMed Google Scholar

41. Klauber-DeMore, N. Tumor biology of breast cancer in young women. Breast Dis. 23, 9–15 (2006).

Article Google Scholar

42. Herranz, M. & Ruibal, A. Optical imaging in breast cancer diagnosis: The next evolution. J. Oncol. 2012, 863747 (2012).

Article PubMed PubMed Central Google Scholar

43. Rankin, S. Recent advances in imaging breast cancer. Cancer Imaging 2, 28 (2001).

PubMed Central Google Scholar

44. Lin, L. et al. Single-breath-hold photoacoustic computed tomography of the breast. Nat. Commun. 9, 1–9 (2018).

Article ADS CAS Google Scholar

45. Fang, Q. et al. Combined optical imaging and mammography of the healthy breast: Optical contrast derived from breast structure and compression. IEEE Trans. Med. Imaging 28, 30–42 (2008).

Article Google Scholar

46. Durduran, T. et al. Bulk optical properties of healthy female breast tissue. Phys. Med. Biol. 47, 2847 (2002).

Article CAS PubMed Google Scholar

47. Heine, J. J. & Thomas, J. A. Effective x-ray attenuation coefficient measurements from two full field digital mammography systems for data calibration applications. Biomed. Eng. Online 7, 13 (2008).

Article PubMed PubMed Central Google Scholar

48. Senchukova, M. A., Nikitenko, N. V., Tomchuk, O. N., Zaitsev, N. V. & Stadnikov, A. A. Different types of tumor vessels in breast cancer: Morphology and clinical value. Springerplus 4, 512 (2015).

Article PubMed PubMed Central CAS Google Scholar

49. Ueda, S. et al. Baseline tumor oxygen saturation correlates with a pathologic complete response in breast cancer patients undergoing neoadjuvant chemotherapy. Cancer Res. 72, 4318–4328 (2012).

Article CAS PubMed PubMed Central Google Scholar

50. Weidner, N., Semple, J. P., Welch, W. R. & Folkman, J. Tumor angiogenesis and metastasis—correlation in invasive breast carcinoma. N. Engl. J. Med. 324, 1–8 (1991).

Article CAS PubMed Google Scholar

51. Tromberg, B. J. et al. Imaging in breast cancer: diffuse optics in breast cancer: detecting tumors in pre-menopausal women and monitoring neoadjuvant chemotherapy. Breast Cancer Res. 7, 279–285 (2005).

Article PubMed PubMed Central Google Scholar

52. Zhou, C. et al. Diffuse optical monitoring of blood flow and oxygenation in human breast cancer during early stages of neoadjuvant chemotherapy. J. Biomed. Opt. 12, 051903 (2007).

Article ADS PubMed CAS Google Scholar

53. Average breast size (2019). https://www.beautytohealth.com/average-breast-size/.

54. Complete anatomy 2020. https://3d4medical.com/.

55. Jesinger, R. A., Lattin, G. E. Jr, Ballard, E. A., Zelasko, S. M. & Glassman, L. M. Vascular abnormalities of the breast: arterial and venous disorders, vascular masses, and mimic lesions with radiologic-pathologic correlation. Radiographics 31, E117–E136 (2011).

Article PubMed Google Scholar

56. Grosenick, D., Rinneberg, H., Cubeddu, R. & Taroni, P. Review of optical breast imaging and spectroscopy. J. Biomed. Opt. 21, 091311 (2016).

Article ADS PubMed Google Scholar

57. Durduran, T. et al. Bulk optical properties of healthy female breast tissue. Phys. Med. Biol. 47, 2847–2861 (2002).

Article CAS PubMed Google Scholar

58. Lei, L. Multi-contrast photoacoustic computed tomography. Dissertation (Ph.D.), California Institute of Technology, https://doi.org/10.7907/FYGX-7907M7929 (2019).

59. Datta, B. N. Numerical linear algebra and applications, (Siam, 2010).

60. Papadakis, E. P., Oakley, C. G., Selfridge, A. R. & Maxfield, B. Fabrication and characterization of transducers. Phys. Acoust. 24, 43–134 (1999).

Article Google Scholar

61. Li, L. et al. Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution. Nat. Biomed. Eng. 1, 0071 (2017).

62. Maggioni, M., Katkovnik, V., Egiazarian, K. & Foi, A. Nonlocal transform-domain filter for volumetric data denoising and reconstruction. IEEE Trans. Image Process 22, 119–133 (2012).

Article ADS MathSciNet PubMed MATH Google Scholar

63. Frangi, A. F., Niessen, W. J., Vincken, K. L. & Viergever, M. A. Multiscale vessel enhancement filtering. MICCAI 1496, 130–137 (Springer, 1998).

64. Nasiriavanaki, M. et al. High-resolution photoacoustic tomography of resting-state functional connectivity in the mouse brain. Proc. Natl. Acad. Sci. USA 111, 21–26 (2014).

Article ADS CAS PubMed Google Scholar

65. Cho, S., Baik, J., Managuli, R. & Kim, C. 3d phovis: 3d photoacoustic visualization studio. Photoacoustics 18, 100168 (2020).

Download references

Acknowledgements

We thank Dr. Vassili Ivanov for useful discussions. This work was sponsored by the United States National Institutes of Health (NIH) grants R01 CA186567 (NIH Director’s Transformative Research Award), R35 CA220436 (Outstanding Investigator Award), and U01 NS099717 (BRAIN2 Initiative).

Author information

1. Author notes

Xiaoyun Yuan

Present address: Department of Electronic Engineering, Tsinghua University, Haidian District, Beijing, China

2. Junhui Shi

Present address: Zhejiang Lab, China Artificial Intelligence Town, Hangzhou, Zhejiang, China

3. These authors contributed equally: Li Lin, Peng Hu, Xin Tong, Shuai Na.

Authors and Affiliations

1. Caltech Optical Imaging Laboratory, Andrew and Peggy Cherng Department of Medical Engineering, California Institute of Technology, Pasadena, CA, USA

Li Lin, Peng Hu, Xin Tong, Shuai Na, Rui Cao, Xiaoyun Yuan, David C. Garrett, Junhui Shi, Konstantin Maslov & Lihong V. Wang

Contributions

L.V.W. and L.L. conceived and designed the study. L.L., D.C.G., J.S. and K.M. constructed the hardware system. S.N., L.L. and X.Y. developed the control program. L.L. and X.T. performed the experiments. R.C. designed the functional rat brain imaging experiments and performed the skull-thinning surgery. P.H. developed the reconstruction algorithm. L.L., X.T. and P.H. analyzed the data. L.V.W. supervised the study. All authors wrote the manuscript.

Corresponding author

Correspondence to Lihong V. Wang.

Ethics declarations

Competing interests

All the authors declare no competing interests. K.M. has a financial interest in MicroPhotoacoustics, Inc. L.V.W. has a financial interest in Microphotoacoustics, Inc., CalPACT, LLC, and Union Photoacoustic Technologies, Ltd., which, however, did not support this work.

Article

- Open Access

- Published:

Visualising peripheral arterioles and venules through high-resolution and large-area photoacoustic imaging

- Yoshiaki Matsumoto,

- Yasufumi Asao,

- Hiroyuki Sekiguchi,

- Aya Yoshikawa,

- Tomoko Ishii,

- Ken-ichi Nagae,

- Shuichi Kobayashi,

- Itaru Tsuge,

- Susumu Saito,

- Masahiro Takada,

- Yoshihiro Ishida,

- Masako Kataoka,

- Takaki Sakurai,

- Takayuki Yagi,

- Kenji Kabashima,

- Shigehiko Suzuki,

- Kaori Togashi,

- Tsuyoshi Shiina &

- Masakazu Toi

Scientific Reports

volume

8, Article number: 14930 (2018) Cite this article

- 4415 Accesses

- 40 Citations

- 2 Altmetric

- Metrics

details

Abstract

Photoacoustic (PA) imaging (PAI) has been shown to be a promising tool for non-invasive blood vessel imaging. A PAI system comprising a hemispherical detector array (HDA) has been reported previously as a method providing high morphological reproducibility. However, further improvements in diagnostic capability will require improving the image quality of PAI and fusing functional and morphological imaging. Our newly developed PAI system prototype not only enhances the PA image resolution but also acquires ultrasonic (US) B-mode images at continuous positions in the same coordinate axes. In addition, the pulse-to-pulse alternating laser irradiation shortens the measurement time difference between two wavelengths. We scanned extremities and breasts in an imaging region 140 mm in diameter and obtained 3D-PA images of fine blood vessels, including arterioles and venules. We could estimate whether a vessel was an artery or a vein by using the S-factor obtained from the PA images at two wavelengths, which corresponds approximately to the haemoglobin oxygen saturation. Furthermore, we observed tumour-related blood vessels around breast tumours with unprecedented resolution. In the future, clinical studies with our new PAI system will help to elucidate various mechanisms of vascular-associated diseases and events.

Introduction

A blood vessel is a passageway for blood and plays a vital role in transporting various components, such as oxygen and nutrients, to each organ of the body. Ailments such as diabetes, cardiovascular disease and collagen diseases, such as rheumatism, involve injury to the peripheral blood vessels. Malignancies such as breast cancer instigate the formation of abnormal blood vessels inside and around tumours by inducing angiogenesis1,2. In addition, when performing plastic surgery to recover tissue that has been lost by surgical removal, as in cancer, it is imperative to transplant not only tissues but also blood vessels from the donor to nourish the transplanted tissue. Thus, it is essential to elucidate the form of blood vessels and arteriovenous veins for disease diagnosis and treatment.

The current vascular imaging method used in the medical field requires the use of a contrast agent, and obtaining an image requires X-ray exposure or an expensive contrast magnetic resonance imaging (MRI) examination. In contrast, blood flow imaging in the B-mode is feasible even with a Doppler imager mounted on an ultrasonic (US) diagnostic apparatus and requires no contrast medium nor invasive procedures; however, its resolution is limited.

Photoacoustic (PA) tomography (PAT) is an imaging method that utilises the light absorption of haemoglobin for non-invasive blood vessel imaging3. Previous clinical studies have actively investigated the impact of PAT on humans for application in the medical field4,5,6,7,8,9,10,11,12,13,14. Instrument configurations in which the linear probe used in a US diagnostic apparatus is provided with a light source have the advantage of enabling the rapid construction of a PA system14,15. However, obtaining high-quality images is challenging because of the ‘limited view problem’ of PA phenomena16. To resolve this problem, a PAT system comprising a hemispherical detector array (HDA) surrounding an object to be measured with a solid angle wider than that of a linear probe has been reported previously17,18,19. We developed a novel prototype (PAI-04) with enhanced specifications to realise large-area, high-resolution blood vessel analysis. Hence, this study aims to outline the device configuration of our novel PA imaging system prototype, presents examples of in vivo images obtained using this device and elucidates its possible clinical applications.

Results

The elementary design of PAI-04 inherited the design of the conventional prototype (PAI-03)9,10,19 equipped with an HDA receiving PA signals. We inserted a test object into the holding cup on the imaging table and performed a scan (Fig. 1; Fig. S1a–c). When imaging was required around the trunk, such as for the breast or thigh, subjects were asked to lie in the prone position on the scanning table. For scanning of the hand, subjects sat beside the table and placed their palm into the holding cup. Then, we scanned the laser light and the HDA helically around the holding cup in a horizontal plane and measured the subjects with a wide angle (Fig. 2a,c). The measurement time was approximately 2 min for an imaging region with a diameter of 140 mm.

Figure 1

figure 1

Schematic illustration showing the PAI-04 system configuration and the alternating irradiation sequence (pulse-to-pulse wavelength switching) adopted in this study. The PA controller controls the laser oscillation and PA wave reception. In this apparatus, laser beams of two different wavelengths are used for alternating irradiation. The generated PA wave is received by the hemispherical probe array (HDA) and sent to the data acquisition system (DAS). In addition to the PA controller, another unit controls ultrasonic (US) transmission and reception. Data received by the US transducer are sent to the US unit. The US transducer and the HDA are integrated as a transducer module, and the entire system is configured to move simultaneously during scanning. The space between the holding cup and the transducer module is filled with circulating water. Water for acoustic matching with the test object is poured into the holding cup.

Full size image

Figure 2

figure 2

Schematic illustration showing the positional correlation between the breast-holding cup and the HDA. (a) A cross-sectional view of the HDA and breast-holding cup cut at the centre during the PA measurement. (b) A cross-sectional view of the HDA and breast-holding cup cut at the centre during the US measurement. (c) A schematic illustration of the top view showing a scan sequence during the PA measurement. The HDA scans in a spiral shape with the central point of the holding cup as a rotation axis to acquire a PA image in a wide range. (d) A schematic illustration of the top view showing a scan sequence of the US measurement after the PA measurement. The US transducer moves under the holding cup and a transducer module and scans in the upward direction from the linear US transducer. When the scanning in the upward direction is completed, the next row is scanned by moving by one transducer width to acquire a wide range of US images.

Full size image

In this prototype, we changed the central frequency of the PA signal detector (Fig. S2a–d) from 2 MHz, used in our previous PAI-03, to 4 MHz. Owing to the phantom experiment, the spatial resolution improved from 0.5 mm to 0.27 mm. In our previous study9, when measuring images using two wavelengths of 756 nm and 797 nm, we performed the second wavelength scan by switching the wavelength after completing the scan using the first wavelength (Fig. S3a–d). In this study, we adopted alternating light irradiation20 with pulse-to-pulse wavelength switching in the PAI-04 (Fig. 1). Another key difference is the simultaneous acquisition of US echo images14,15,21,22. We placed a dedicated 256-channel linear transducer, equivalent to that used in a typical clinical setting, beside the HDA for US echo imaging. After performing a spiral scan of the HDA to obtain the PA signal (Fig. 2a,c), we moved the US probe under the subjects and scanned automatically (Fig. 2b,d). The three-dimensional US (3D-US) image acquisition time was approximately 1 min. Then, we reconstructed 3D-US images by stacking B-mode images acquired in a broad range. The 3D-US image acquisition time was approximately 1 min. Furthermore, the data acquisition area of the 3D-US in the x–y plane comprised a square inscribed in a circle with a diameter of 140 mm.

We evaluated the haemoglobin oxygen saturation (SO2) using an approximate parameter calculated from measurements obtained at two wavelengths, which we termed the S-factor. Background optical coefficients were measured beforehand for scanning mammary glands, and we attempted a numerical analysis of the S-factor. Notably, we did not measure the background optical coefficients for other tissues for anatomical reasons. In regions where the background was considered similar, such as for adjacent blood vessels, we projected whether the correlation between the S-factors along the two vessels was consistent with the correlation between their real SO2 values and compared the relative values. We acquired a 3D-US image to measure a breast tumour and extracted and stacked the hypoechoic area as a volume image. In addition, the tumour lesion was identified by a breast physician familiar with mammary US imaging. Then, segmentation processing was performed on the basis of the US image intensity. Furthermore, we applied the graph cut method23,24 to extract the tumour area.

Figure 3a–d present examples of image acquisition in the palm (Fig. 3a–c; Fig. S3a–d) and thigh (Fig. 3d) of a healthy subject. Figure 3e shows the scanned position in the living body. We confirmed the effect of alternating irradiation using the image of the palm (Fig. 3a–c; Fig. S3a–d). The S-factor colour of the common palmar digital arteries designated by A1–3 in Fig. 3c (by alternating laser irradiation) was uniform, and two accompanying vessels ran parallel to one another. Conversely, the colour was not uniform in Fig. S3d (conventional sequential irradiation), and the difference between the artery and the vein was not easily discerned. In addition, arteries and veins could be clearly distinguished by nuances in colour in the S-factor image obtained using alternating irradiation in the PAI-04. Figure 3d presents an image of an anterolateral thigh. The correlation between the locations of perforators (P) identified by PA imaging compared with conventional US exhibited discrepancies of <10 mm11. In addition, in PA images of bundles of blood vessels near the stem part, the image features comprising the vessel presumed to be the artery (A) and the vessel presumed to be the vein (V) from the image matched the anatomy of the perforators.

Figure 3

figure 3

Examples of palm and thigh PA images obtained from healthy subjects in vivo. (a–c) Examples of palm PA images obtained using an alternate irradiation sequence. No body motion correction was performed. (a) The maximum intensity projection (MIP) image of the whole palm. (b) An image after deletion of the subcutaneous veins from the whole palm image. (c) An enlarged image of the region of (b) indicated by the white dashed line shows the common palmar digital arteries. In Fig. 2c, the blood vessels anatomically determined to be the common palmar digital arteries are designated A1–3, and the veins accompanying them are V1–3. The assignment of these arteries and veins corresponded to the blood vessel colour of the S-factor image (i.e., the magnitude relationship of the S-factor). (d) An example PAT image obtained in an anterolateral thigh. A stem portion of perforator (P) vessels and a bundle representing an artery (A) and a vein (V) were observed. (e) Schematic illustration of the measured tissue in the body of a subject.

Full size image

Figure 4a–c present images of a breast. As mentioned earlier, as the background optical coefficients of subjects were measured, we attempted to quantify the S-factor. In addition, the demarcation of arteries and veins by the colour tone (S-factor image) corresponded to the determination of vessels by the anatomical findings. We calculated the S-factors along three different adjacent vessels (Fig. 4a) and performed three measurements to evaluate the reproducibility. Consequently, the coefficient of variation (CV%) of the S-factor of the arteries and veins was <2% (Fig. S4). Then, we evaluated 23 adjacent arteriovenous S-factors (Fig. 4b). A significant difference was observed (P = 9.93518E-19) between the S-factor values of the arteries and veins (Table S1), and we considered these absolute values to be reasonable results25. Figure 4d,e present examples of the results of quantifying the S-factors of arteries and veins and show the S-factor values for each voxel around A7/V7 and A8/V8, as shown in Fig. 4b,c, respectively. Thus, although the S-factor values evaluated at each point of the blood vessel varied, the range was considered small enough to enable determination of whether they were arteries or veins.

Figure 4

figure 4

Examples of breast PA images obtained from healthy subjects in vivo. (a) Three measured positions for the reproducibility evaluation of the S-factor in the breast. The S-factors of three adjacent vessels (A4–6, V4–6) are shown. The measured ranges of the S-factor are indicated by the yellow dotted lines. (b) Twenty-three measured positions for evaluating the variation in the S-factor values. (c) Examples of the S-factor evaluation of two sets of adjacent vessels (A7/V7 and A8/V8) indicated by a yellow arrowhead in Fig. 3b. The S-factor value in a voxel along a vessel was defined as the average value of the signal at a distance within 2 voxels from the centre of the voxel with a sufficient PA signal intensity to omit background signals. Consequently, the S-factor value calculated at each point of a vessel fluctuated, although the range was very small.

Full size image

Figure 5a–f show examples of the results obtained for a breast cancer lesion. Figure 5a presents a fusion image of the PA obtained at 797 nm and the 3D-US images (red colour). Notably, fine blood vessels surrounding the tumour are clearly visible. Figure 5b presents a US echo B-mode image of a tumour region. We observed a hypoechoic lesion suggesting a breast tumour in the region, as indicated by an arrow. Figure 5c presents a fusion image of the S-factor and 3D-US images (red colour), and Fig. S5 presents the S-factor and the original image of US C-mode around the tumour. Figure 5d presents the S-factor image of the same region with enlargement of the tumour vicinity and without fusion with the US tumour image. Of note, the arterioles and venules show clustering. In addition, yellow- and green-coloured fine vessels have S-factors of approximately 80–90%, suggesting veins. Figure 5e,f show axial images of the S-factor at the edge and the central portion of the tumour, respectively. In this case, a few blood vessels surrounding or inside the tumour exhibited the oxygen saturation of arteries, and several blood vessels with S-factors corresponding to veins surrounded the tumour.

Figure 5

figure 5

Examples of images obtained from a breast cancer patient in vivo. (a) A fusion image of the PA image taken at 797 nm and 3D-US images. Hypoechoic portions of the US were extracted based on the US signal intensity and reconstructed into a 3D volume by stacking the continuous B-mode images. The cubed US data were coloured red. Several clustered fine blood vessels were observed in the vicinity of the tumour. (b) A US B-mode image of the dotted-line portion indicated by 1–1′ in (a). A hypoechoic region is visible and is indicated by the yellow arrowhead. (c) A fusion image of the S-factor and 3D-US images (red colour). The nipple at the centre of this image is coloured light blue because the absorption spectrum of melanin in the nipple and the spectrum of low oxygen saturation are almost the same. (d) An S-factor image of the same region with enlargement of the tumour vicinity but without fusion with the US tumour image. (e) The entire MIP image in the axial direction of the S-factor around the lesion. (f) The slab MIP image at the centre of the lesion corresponding to (e).

Full size image

In addition, we observed some scattered hypoxic points, which most likely represented intratumoural bleeding. The measured diameters of the spots were 700 to 900 µm, as indicated by the yellow arrows in Fig. S6a. This patient underwent tumour resection after PAI observation. The pathology specialist microscopically observed a haematoxylin and eosin (H&E)-stained specimen of the tumour tissue and a CD31-stained specimen of the same area (Fig. S6b–g), which revealed a fibrous core at the centre of the tumour and indicated the decline in endothelial cells of the blood vessels. The central portion of the tumour, where the PA image displayed low and variable SO2, might be avascularised by necrosis and fibrosis.

Discussion

The PAI-04 enabled the visualisation of high-resolution images of arterioles and venules. The improvement in the resolution over that of the conventional prototype could be attributed to the enhanced central frequency. However, it is imperative to carefully consider the use of higher frequency depending on the object measured and the desired application because when the penetration depth decreases because of US attenuation or when the relative bandwidth of the detector is limited, the concern arises that the receivable low-frequency signal is limited and that image quality with good morphological reproducibility cannot be obtained. In addition, with alternating laser irradiation to enhance the S-factor reproducibility, blood vessels anatomically identified as arteries or veins were depicted as relative colours in the S-factor image. Thus, even if the light fluence is unknown, the magnitude relationship of the oxygen saturation can be clinically assessed, demonstrating the validity of the hypothesis that the S-factor reflects the relative SO2 value in the neighbouring region.

In thigh imaging, the PA imaging contributed to visualising not only the branching morphology of perforators but also the arteries and veins in the subcutaneous layer. Precise information on subcutaneous vessels will enable surgeons to plan safer procedures for anterolateral thigh flap surgery. In breast cancer imaging, small tumour-related blood vessels could be visualised. Although in vivo imaging of tumour-related blood vessels has been reported at the xenograft level2, to the best of our knowledge, this study is the first to report this level of detail for tumour-related blood vessels in human cancer tissue. The imaging of fine vessels could be attributed to the fact that the tumour in this study was located near the body surface, with a depth <10 mm from the skin. In the future, we intend to extend the imaging technology for tumour-related blood vessels to deeper positions by enhancing the sensor sensitivity and optimising the laser irradiation method. Because breasts are relatively easy to deform, developing suitable breast positioning to obtain clear PA images using this device is a major goal. We will soon report the results of systemic therapy based on changes in the S-factor in tumour-related blood vessels.

The 3D-US image was markedly effective in identifying the tumour position. By adopting a combination with a US diagnostic system using a linear probe, it was relatively easy to compare the B-mode image usually obtained for diagnosis at the clinical site with the PA image. Conversely, when combining a separate system with a linear probe, a scanning time difference of approximately 2 min occurs, as the US probe is moved after PA imaging. If the subject moves during this period, there is no way to correct for the discrepancy. Thus, physicians should monitor patient movement during the scanning period and carefully assess the correlation between morphology and metabolic function. More reliable data is expected to be obtained from future prototypes by adopting a system configuration that performs the reception of PA signals and US transmission/reception simultaneously7,14,15,21,22. In addition, new clinical knowledge will be generated by improving the device and advancing the analytical methods. As our device can be used to comprehensively analyse the vessel branching structure and oxygen metabolism of the human body, it could contribute to various biological assessments in vivo and provide clinically important findings. For example, it might be feasible to observe the normalisation of blood vessels in response to drug treatment non-invasively. Perhaps the effect of drug treatment could be evaluated at an early stage, which could contribute to individualised medicine. Furthermore, analysis combining PAI and conventional modalities promises to create a new paradigm for disease diagnosis.

In conclusion, we have developed a high-resolution PAT system prototype capable of laser irradiation with pulse-to-pulse wavelength switching that could obtain 3D-PA and 3D-US images. In addition, arterioles and venules could be imaged and distinguished using S-factor images. In the future, clinical studies with a high-resolution PAI-04 system will help to elucidate various aspects of vascular-associated diseases and events.

Methods

Device configuration

We used the PAI-04 system with an HDA (Canon Inc., Tokyo, Japan). We placed the tissue of the subjects intended for measurement in the breast-holding cup on an exam table (Figs 1 and S1). In this study, we referred to the x–y Cartesian plane as the area parallel to the horizontal plane or parallel to the plane of the exam table surface. The x- and y-axes were the short and long axis of the table, respectively, and the z-axis was the vertical axis. We used one laser system per wavelength to emit laser light with a very large pulse. Although a prior study26,27 reported the advantage of spectrum analysis employing multiple wavelengths, we adopted a PAT system with two wavelengths to aim for future practical applications with a straightforward configuration, as in our conventional system7,9. In addition, we verified whether the oxygen saturation of the blood in the living body could be measured using this two-wavelength system. We used two Ti:Sa lasers optically pumped using a Q-switched Nd:YAG laser, which can select wavelengths from 700 to 900 nm. We selected two wavelengths, 756 and 797 nm, to construct an improved haemoglobin saturation distribution image. While 756 nm is the local maximum of the absorption coefficient of deoxidised haemoglobin, 797 nm is the isosbestic point of oxyhaemoglobin and deoxidised haemoglobin.

In this study, the laser lights of these two wavelengths alternated in every irradiation pulse along the z-axis. Notably, the maximum permissible exposure values recommended by the American National Standards Institute for 756 and 797 nm are 12.9 and 15.6 mJ/cm2, respectively. The energy of each laser used in the PAM-04 was 89 and 89 mJ/pulse for 756 and 797 nm, respectively. We set the light energy density for 756 and 797 nm to maximum values of 10.0 and 11.9 mJ/cm2, respectively. Furthermore, an 8-mm-diameter pulsed laser beam was fed through a light fibre that directed the laser beam upward along the vertical axis of the transducer array as the HDA was scanned. A diverging lens placed at the base of the HDA spread the light conically to a diameter of 54 mm at the surface of the holding cup. Each laser light source corresponding to each wavelength was emitted sequentially at 10 Hz. Thus, the pulse sequence of (1) pulse laser irradiation at 797 nm, (2) irradiation at 756 nm after 50 ms, and (3) irradiation at 797 nm after 50 ms was repeated.

We used an internally manufactured broadband capacitive micromachined US transducer (CMUT) as a receive-only transducer element to detect the PA signal. In addition, we constructed the detector array for PA signal detection to correspond to the HDA, and its appearance was almost the same as in our prior study (PAI-03)9,19. Overall, 500 spirally arrayed CMUT elements were present on the HDA; each element measured 1.5 mm in diameter, and the elements were set at approximately 10-mm intervals. The central frequency and fractional bandwidth of the CMUT element were 4 MHz and >100%, respectively. Notably, the average conversion efficiency was 63 mV/kPa at 4 MHz, and its standard deviation was 2.5 mV/kPa over 500 elements. Furthermore, the noise-equivalent pressure (NEP) of this system was 0.5 Pa at the root-mean-square value without signal averaging. We examined the NEP measurement with a low-pass filter with a cut-off frequency of 15 MHz.

The key difference between PAI-04 and PAI-03 was the addition of US echo imaging. PAI-04 used a linear transducer array as a dedicated transducer for US echo imaging because the CMUT device used in the HDA was designed as a receive-only device. In addition, the 256-channel linear US probe was placed adjacent to the HDA. While the US beam was emitted along the z-axis, the electronic scanning direction of the linear probe was along the x-axis. Using probes of the same shape as in routine clinical practice facilitated the interpretation of US images. Water was used between the HDA and the breast-holding cup for acoustic matching and was supplied from a water circulation device set at a predetermined temperature. A drain outlet was provided on the bottom of the HDA. The subjects were inserted into the body-temperature water and scanned. For breasts, the affected breast was scanned first, followed by the contralateral side. Lesions were identified primarily using the other hand-held US machine.

Data acquisition

The HDA was scanned helically in the x–y plane with the same trajectory as described previously19 to obtain a PA image of a large area of the subject. It did not move in the z-axis direction. We acquired a PA signal in a range with a diameter of 140 mm using a spiral scan. The received PA signal data were amplified by a variable gain amplifier (VGA) capable of settings up to 10–42 dB and converted to digital data with a 40 MHz, 12-bit AD converter. We primarily set the amplification rate to 30 dB in this clinical study.

For a 140-mm-diameter area, the measurement time for a PA signal was approximately 2 min. To obtain a US image after PA scanning, we obtained US B-mode images at continuous positions by moving the US linear probe under the subjects and scanning along the y-axis. Then, a 3D-US image was obtained by stacking the B-mode images. The probe was shifted along the x-axis after the one-line scan and the same scan was performed along the y-axis to measure a range greater than the probe width in the x-axis direction. Then, the scanned US data were combined after the measurement and converted into a 3D-US image. The measurement time for a US image was approximately 1 min for a 100-mm-square area. When examining patients with breast cancer, the state of the tumour was confirmed beforehand with a clinical US apparatus, and a breast surgeon selected an appropriate scan area based on the anatomy of the patient’s breast and the extent of the tumour spread. In cases other than breast cancer, we did not obtain 3D-US echo images and evaluated only PA images.

Imaging method

We used universal back-projection (UBP)28 for the PA image reconstruction. The reconstructed voxel size was typically 0.0625 or 0.125 mm and was selected depending on the size of the object and the type of analysis. We created an image for each laser shot and constructed the whole image by combining these images (Supplemental Video 4). Since the water temperature differed between the water above and below the holding cup, we performed image reconstruction by evaluating the sound velocity, which was categorised into two layers based on the difference in the sound velocity because of the temperature difference29.

When imaging a breast, we measured the background optical parameter beforehand with the TRS-20 (Hamamatsu Photonics K.K., Hamamatsu, Japan)30 and obtained an absorption coefficient distribution by roughly evaluating the light fluence distribution. In this study, we hypothesised that the absolute value of SO2 could be approximately evaluated. Conversely, when tissues other than the breast were measured (e.g., the palm), the derivation of the background optical constant was challenging because tissues, such as muscles, tendons and bones, were intertwined with each other. Thus, we performed S-factor analysis based only on the measured initial sound pressure distribution without evaluating the absorption coefficient distribution. Even in that case, the magnitude correlation between the true oxygen saturation and the measured oxygen saturation, which we called the S-factor, was not reversed. Hence, we did not consider the absolute value of oxygen saturation as a subject of discussion and instead considered it as a relative value. The calculation of the S-factor and its validity for relative assessment are described later. In addition, we did not evaluate the light fluence except in the breast evaluation and obtained images using only the measured initial acoustic pressure distribution, with the light fluence ratio assumed to be 1; this calculation or assumption was expressed in the S-factor image. For S-factor images in this study, the S-factor value and the absorption coefficient at 797 nm were indicated by the hue and intensity, respectively.

In this study, as in our previous study, we used two types of body motion correction methods for image reconstruction9. The first was a correction for body motion between different laser pulses (p-BMC), and the second was a correction between two different laser wavelengths (w-BMC) for the S-factor image. We used p-BMC and not w-BMC in this study, except for the comparative experiments (Fig. S3). We selected p-BMC because the movement of the living body at rest during the 50-ms pulse interval was assumed to be negligible, even in the case of the high-resolution PA image, where the resolution was several hundred micrometres. In addition, when scanning the breast, a 3D-US image was acquired using the scanning sequence (Fig. 2c). We constructed 3D-US images by scanning the ultrasonic linear array probe along the elevation direction and creating an image in which successive B-mode images were stacked.

The coordinates of the constructed 3D-US data were set to coincide with the coordinates of the PA image so that two could be overlaid on an image viewer31. We extracted the low-echo area of the US based on the intensity information from the US image to compare and analyse the obtained 3D-US image and PA image. Regarding mammary tumours, physicians familiar with the US diagnosis of breast cancer manually assigned a hypoechoic region inside the tumour. Then, by applying the graph cut method23,24 with the hypoechoic region as the seed point, we extracted the curved boundary surface between the tumour and the nontumour region to derive the 3D structure. After binary processing with the extracted tumour area set to 1 and other areas set to 0, we fused the US image of the tumour volume alone in the red image with the PA image. For the S-factor, colour imaging was adopted in which the value of the S-factor corresponded to the hue and the image intensity value to the absorption coefficient at 797 nm, which provided an image approximating the weight distribution of the total haemoglobin. After obtaining the PA image, we deleted the data near the subcutaneous tissue by image processing as necessary. Next, we applied a cloth simulation method to obtain information on the body surface31,32. Furthermore, the volume data were trimmed from the body surface position to an arbitrary depth to simplify the image interpretation of deeper tissues.

Error between the S-factor and the oxygen saturation

We adopted the following method to evaluate the haemoglobin oxygen saturation using a two-wavelength light source. However, quantitative PA imaging is a major problem that requires a solution. In particular, quantitatively and precisely evaluating the absolute value of the absorption coefficient (μa) is difficult because of the limited view problem7 along with difficulty in estimating the light fluence. If the two μa values of the absorber with wavelengths λ1 and λ2 are obtained correctly, the Hb saturation value (SO2) can be derived from the following equation:

SO2=[HbO2][HbO2]+[Hb]=μλ2a(r)μλ1a(r)⋅ελ1Hb−ελ2Hbελ2ΔHb−μλ2a(r)μλ1a(r)⋅ελ1ΔHb

(1)

where λ1 and λ2 represent wavelengths of 756 and 797 nm, respectively; r is the position to be calculated; εHb is the molar extinction coefficient of deoxy-Hb; and εΔHb is the difference in the molar extinction coefficients between deoxy-Hb and oxy-Hb (HbO2). However, in practice, the absolute value of SO2 is unreliable because the light fluence affects the absolute SO2 value. Hence, we adopted the following method to extract the semiquantitative information in vivo.

To evaluate μa, Eq. 2 was derived as follows:

μλ1a(r)=pλi0(r)Γ⋅φλi(r),i=1,2

(2)

where p0, Γ and ϕ are the initial pressure of PA, the Grüneisen parameter and the light fluence, respectively. In addition, λ1 and λ2 represent wavelengths of 756 and 797 nm, respectively.

Eq. 1 can be transformed as follows:

SO2(r)=pλ20(r)pλ10(r)⋅φλ1(r)φλ2(r)⋅ελ1Hb−ελ2Hbελ2ΔHb−pλ20(r)pλ10(r)⋅φλ1(r)φλ2(r)⋅ελ1ΔHb

(3)

Here, the correlation between the reconstructed initial pressure pmeas(r) and the true initial pressure p0(r) can be expressed as follows:

pλimeas(r)=pλi0(r)⋅Uλi(r)+Aλi(r),i=1,2

(4)

where U is a correction term expressing the variation in the acoustic properties accompanying the propagation and reception of PA waves because of acoustic diffusion, acoustic attenuation, the limited view problem (caused by the small receiving angle), the central frequency and the fractional bandwidth of the probe. In addition, the acoustic correction term U can be considered to have approximately the same value for both wavelengths except for cases where the shape of the initial pressure or the frequency characteristics of the PA signals differ between the two wavelengths because of the extreme wavelength dependence of light absorbance. A is a term expressing noise and artefacts. We used the term S-factor (Sf) instead of oxygen saturation (SO2) for the approximation to avoid misunderstanding. Based on this information, Eq. 3 can be approximated as follows:

Sf(r)=M(r)⋅ελ1Hb−ελ2Hbελ2ΔHb−M(r)⋅ελ1ΔHb

(5)

M(r)=pλ2meas(r)−Aλ2(r)pλ1meas(r)−Aλ1(r)⋅φλ1(r)φλ2(r)

(6)

In Eq. 6, the initial acoustic pressure at each wavelength can be evaluated by the image reconstruction of UBP using the signal obtained by the PAT system. Thus, the S-factor accuracy depends on estimation of the ratio between the light fluence of the two wavelengths in M (r) and the amount of noise and artefacts. When the ratio between the light fluence of the two wavelengths can be regarded as uniform, the following equation is derived:

∂(Sf)∂(SO2)=∂(Sf)∂(φλ1(r)φλ2(r))⋅1∂(SO2)∂(φλ1(r)φλ2(r))

(7)

={φλ1(r)φλ2(r)⋅ελ2ΔHb−pλ20(r)pλ10(r)⋅ελ1ΔHbφλ1(r)φλ2(r)⋅ελ2ΔHb−pλ2meas(r)−Aλ2(r)pλ1meas(r)−Aλ1(r)⋅ελ1ΔHb}2⋅pλ10(r)pλ20(r)⋅pλ2meas(r)−Aλ2(r)pλ1meas(r)−Aλ1(r)

(8)

={φλ1(r)φλ2(r)⋅ελ2ΔHb−pλ20(r)pλ10(r)⋅ελ1ΔHbφλ1(r)φλ2(r)⋅ελ2ΔHb−pλ2meas(r)−Aλ2(r)pλ1meas(r)−Aλ1(r)⋅ελ1ΔHb}2⋅pλ10(r)pλ20(r)⋅Aλ2(r)Aλ1(r)⋅SNRλ2(r)−1SNRλ1(r)−1

(9)

In Eq. 9, the first, second and third terms on the right side are always positive. Thus, ∂(Sf)/∂(SO2) can be considered always positive from term 4, at least when the signal-to-noise ratios (SNRs) are >1 for both wavelengths. In other words, when the ratio between the light fluence of the two wavelengths of subjects is constant (i.e., when the tissue can be assumed to be a homogeneous medium), the magnitude relationship of the oxygen saturation of the obtained blood vessel images can be assessed. In contrast, when the ratio of the light fluence of the subjects is not considered constant, even relative comparisons require attention. For example, comparing S-factors is challenging for far-distant blood vessels because the light fluence is altered by the impact of the tissue on the light paths leading to the respective blood vessels. Conversely, when comparing blood vessels in close proximity to each other, such as adjacent arteries and veins, the light fluence reaching the area is assumed to be constant, and the S-factors can be evaluated. Figure S7 presents the error rate of the light fluence ratio and the numerical value of the S-factor evaluated in Eq. 5. The noise component was ignored in this study. The five lines in the graph represent numerical values when the correct oxygen saturation value is 100%, 90%, 80%, 70% and 60%. This graph illustrates the validity of Eq. 9 because the magnitude relation of the S-factor is upheld if there is an error in light fluence. Furthermore, the gradient of the error as the oxygen saturation increases is gradual.

Data analysis

As mentioned earlier, even if the fluence of light is unknown, the magnitude correlation between the actual SO2 and the S-factor of the blood vessel image does not alter. Thus, we attempted a relative assessment of the S-factor in the limb measurement to distinguish between the vein and the artery. In this clinical study, it was not feasible to invasively confirm either an artery or a vein in humans because the premise was non-invasive measurement. Thus, whether a blood vessel was an artery or a vein was estimated depending on the morphological knowledge and the relationship of the haemoglobin oxygen saturation. Morphologically, it is anatomically known that veins accompany arteries; that is, one or two veins and one artery are often adjacent. Furthermore, the haemoglobin oxygen saturation of venous blood is lower than that of arterial blood25. In this study, a surgeon presumed either an artery or a vein based on the above knowledge. In the clinical setting, we would expect to visually confirm whether the vessel was an artery or a vein at the time of the skin tissue incision.

Regarding clinical studies for breast cancer diagnosis, the background optical coefficient and the light fluence distribution in the breast were evaluated in advance based on the TRS-20 measurement results described earlier. We verified the validity of the quantitative S-factor value as follows. First, the reproducibility was evaluated. Then, we arbitrarily selected arteries and veins adjacent to each other in a range of at least 10 mm at three different positions in a normal breast and evaluated their average S-factor. After confirming that the condition of the subjects was stable, we performed three measurements and obtained the average value and CV% of the S-factor. Next, we selected 23 sets of adjacent vessels (associated with arteries and veins) capable of tracing over 15 pixels in the same breast. Then, the mean value of the S-factor of the blood vessel in that range was obtained and compared between the arteries and veins. In addition, we assessed whether a significant difference was detected between the arteries and veins using Welch’s t-test to detect differences between S-factors. Regarding the validity of the absolute value, we used the average value of a pulse oximeter for the artery. Finally, we validated the vein by considering the average value of the mixed venous blood oxygen saturation measured by a Swan-Ganz catheter33 and the average value of the oxygen saturation at the central vein.

Postoperatively, we sectioned excised specimens at 5 mm intervals perpendicular to the longest axis of a specimen and performed pathological analysis on permanent formalin-fixed paraffin-embedded (FFPE) tissues by conventional H&E staining. In addition, we used the anti-CD31 monoclonal antibody JC70A (Dako, Santa Clara, CA 95051) to visualise microvessels by immunohistochemistry at a dilution of 1:50 overnight after heat-induced antigen retrieval in Dako Target Retrieval solution pH 9 for 30 min, according to the manufacturer’s instructions.

Ethics

This study was approved by the Ethics Committee of the Kyoto University Graduate School of Medicine (UMIN000018893, UMIN000022215 and UMIN000022767), and we obtained written informed consent from all subjects. This study was conducted in accordance with the Declaration of Helsinki. For UMIN000022215, female patients with a breast lesion or lesions who fulfilled all inclusion criteria and who did not meet the exclusion criteria were enrolled. The subjects of the UMIN000018893 and UMIN000022767 studies were healthy volunteers whose palms or thighs were imaged using our PA system. Detailed results of each clinical trial will be reported in a future article.

References

1. Jain, R. K. & Carmeliet, P. F. Vessels of death or life. Sci. Am. 285(6), 38–45 (2001).

Article CAS Google Scholar

2. Vakoc, B. J. et al. Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging. Nat. Med. 15(10), 1219 (2009).

Article CAS Google Scholar

3. Wang, L. V. & Hu, S. Photoacoustic Tomography: In Vivo Imaging from Organelles to Organs. Science. 335(6075), 1458 (2012).

Article ADS CAS Google Scholar

4. Heijblom, M. et al. The state of the art in breast imaging using the Twente Photoacoustic Mammoscope: results from 31 measurements on malignancies. Eur. Radiol. 26, 3874 (2016).

Article Google Scholar

5. Kitai, T. et al. Photoacoustic mammography: initial clinical results. Breast Cancer. 21(2), 146 (2014).

Article Google Scholar

6. Fakhrejahani, E. et al. Clinical Report on the First Prototype of a Photoacoustic Tomography System with Dual Illumination for Breast Cancer Imaging. Plos One 10(10), e0139113 (2015).

Article Google Scholar

7. Asao, Y. et al. Photoacoustic mammography capable of simultaneously acquiring photoacoustic and ultrasound images. J. Biomed. Opt. 21(11), 116009 (2016).

Article ADS Google Scholar

8. Diot, G. et al. Multispectral Optoacoustic Tomography (MSOT) of Human Breast Cancer. Clin. Cancer Res. 23(22), 1 (2017).

Article Google Scholar

9. Toi, M. et al. Visualization of tumor-related blood vessels in human breast by photoacoustic imaging system with a hemispherical detector array. Sci. Rep. 7, 41970 (2017).

Article ADS CAS Google Scholar

10. Matsumoto, Y. et al. Label-free photoacoustic imaging of human palmar vessels: a structural morphological analysis. Sci. Rep. 8, 786 (2018).

Article ADS CAS Google Scholar

11. Tsuge, I. et al. Photoacoustic Tomography Shows the Branching Pattern of Anterolateral Thigh Perforators In Vivo. Plast. Reconstr. Surg. 141, 1288 (2018).

Article CAS Google Scholar

12. Ishida, Y et al. Photoacoustic imaging system visualizes restoration of peripheral oxygenation in psoriatic lesions. J. Eur. Acad. Dermatol. Venereol. 15032 (2018).

13. Shiina, T. et al. Development and clinical translation of photoacoustic mammography. Biomed. Eng. Lett. 8(2), 157 (2018).

Article Google Scholar

14. Horiguchi, A. et al. Pilot Study of Prostate Cancer Angiogenesis Imaging Using a Photoacoustic Imaging System. Urology. 108, 212 (2017).

Article Google Scholar

15. Kim, C. H. et al. Handheld array-based photoacoustic probe for guiding needle biopsy of sentinel lymph nodes. J. Biomed. Opt. 15(4), 046010 (2010).

Article ADS Google Scholar

16. Xu, Y. & Wang, L. V. Reconstructions in limited-view thermoacoustic tomography. Med. Phys. 31, 724 (2004).

Article Google Scholar

17. Deán-Ben, X. L. et al. Adding fifth dimension to optoacoustic imaging: volumetric time-resolved spectrally enriched tomography. Light-Sci. Appl. 3, e137 (2014).

Article Google Scholar

18. Deán-Ben, X. L. et al. Functional optoacoustic human angiography with handheld video rate three dimensional scanner. PACS. 1, 68 (2013).

PubMed Google Scholar

19. Kruger, R. A. et al. Dedicated 3D photoacoustic breast imaging. Med. Phys. 40, 113301 (2013).

Article Google Scholar

20. Deán-Ben, X. L. et al. Functional optoacoustic imaging of moving objects using microsecond-delay acquisition of multispectral three-dimensional tomographic data. Sci. Rep. 4, 5878 (2014).

Article Google Scholar

21. Merčep, E. et al. Hybrid Optoacoustic Tomography and Pulse–Echo Ultrasonography Using Concave Arrays. IEEE trans. ultrason. ferroelectr. freq. control. 62(9), 1651 (2015).

Article Google Scholar

22. Merčep, E. et al. Combined Pulse-Echo Ultrasound and Multispectral Optoacoustic Tomography With a Multi-Segment Detector Array. IEEE trans. med. imag. 36(10), 2129 (2017).

Article Google Scholar

23. Boykov, Y. Y. & Jolly, M. P. Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images. Proc. of ICCV. I, 105 (2001).

24. Boykov, Y. Y. et al. Segmentation of Dynamic N-D Data Sets via Graph Cuts using Markov Models. Proc. of MICCAI. LNCS 2208, 1058 (2001).

25. Walley, K. R. Use of Central Venous Oxygen Saturation to Guide Therapy. Con. Clin. Rev. 184, 514 (2011).

Google Scholar

26. Deán-Ben, X. L. et al. Fast unmixing of multispectral optoacoustic data with vertex component analysis. Opt. Laser Eng. 58, 119 (2014).

Article ADS Google Scholar

27. Tzoumas, S. et al. Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues. Nat. Commun. 7, 12121 (2016).

Article ADS CAS Google Scholar

28. Xu, M. & Wang, L. V. Universal back-projection algorithm for photoacoustic computed tomography. Phys. Rev. E. 71, 016706 (2005).

Article ADS Google Scholar

29. Deán-Ben, X. L. et al. Accounting for speed of sound variations in volumetric hand-held optoacoustic imaging. Front. Optoelectron. 10(3), 280 (2017).

Article Google Scholar

30. Oda, M. et al. Near-infrared time-resolved spectroscopy system for tissue oxygenation monitor. Proc. SPIE 4160, 204 (2000).

Article ADS Google Scholar

31. Sekiguchi, H. et al. Body surface detection method for photoacoustic image data using cloth-simulation technique. Proc. SPIE. 1049459 (2018).

32. Baraff, D. & Witkin, A. Large Steps in Cloth Simulation. Proc. of SIGGRAPH. 43 (1998).

33. Bernard, G. R. et al. Pulmonary Artery Catheterization and Clinical Outcomes. JAMA 283(19), 2568 (2000).

Article CAS Google Scholar

Download references

Acknowledgements

This work was partially supported by the Innovative Techno-Hub for Integrated Medical Bioimaging Project of the Special Coordination Funds for Promoting Science and Technology from the Ministry of Education, Culture, Sports, Science, and Technology, Japan. This work was partially funded by the ImPACT Program of the Council for Science, Technology and Innovation (Cabinet Office, Government of Japan). The authors thank T. Kosaka and K. Kobayashi for coordinating the activities for this clinical research and Canon, Inc., for providing the PAI-04 system.

Author information

1. Author notes

Yoshiaki Matsumoto and Yasufumi Asao contributed equally.

Authors and Affiliations

1. Department of Breast Surgery, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Yoshiaki Matsumoto, Yasufumi Asao, Aya Yoshikawa, Tomoko Ishii, Masahiro Takada & Masakazu Toi

2. Japan Science and Technology Agency, ImPACT Program, Cabinet Office, Kyoto, Japan

Yasufumi Asao & Takayuki Yagi

3. Department of Diagnostic Imaging and Nuclear Medicine, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Hiroyuki Sekiguchi, Masako Kataoka & Kaori Togashi

4. Medical Imaging Development Center, Canon Inc., Kyoto, Japan

Ken-ichi Nagae & Shuichi Kobayashi

5. Department of Plastic and Reconstructive Surgery, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Itaru Tsuge, Susumu Saito & Shigehiko Suzuki

6. Department of Dermatology, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Yoshihiro Ishida & Kenji Kabashima

7. Department of Diagnostic Pathology, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Takaki Sakurai

8. Department of Human Health Science, Graduate School of Medicine, Kyoto University, Kyoto, Japan

Tsuyoshi Shiina

Contributions

Conceived and designed the experiments: Y.M., Y.A., K.N., S.K., I.T., S. Saito, M. Takada, Y.I., M.K., T.Y., K.K., K.T., S. Suzuki and M. Toi. Performed the experiments: Y.M., K.N., S.K., I.T., S. Saito and Y.I. Analysed the data: Y.M., Y.A., H.S., A.Y., T.I., I.T., S. Saito, Y.I., T.S., and T.Y. Wrote the paper: Y.M. and Y.A.

Corresponding author

Correspondence to Masakazu Toi.

Ethics declarations

Competing Interests

Yasufumi Asao, Ken-ichi Nagae, Shuichi Kobayashi and Takayuki Yagi are employees of Canon Inc., Japan. Y.A. and T.Y. temporarily transferred to the Japan Science and Technology Agency. Canon, Inc., designed and invented the photoacoustic imaging system used in this study. The other authors have no conflicts of interest.

Additional information

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Supplementary Video 1

Supplementary Video 2

Supplementary Video 3

Supplementary Video 4

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and Permissions

Non-invasive photoacoustic computed tomography of rat heart anatomy and function

- Li Lin,

- Xin Tong,

- Susana Cavallero,

- Yide Zhang,

- Shuai Na,

- Rui Cao,

- Tzung K. Hsiai &

- Lihong V. Wang

Light: Science & Applications

volume

12, Article number: 12 (2023) Cite this article

- 1388 Accesses

- 1 Citations

- 19 Altmetric

- Metrics

details

Abstract

Complementary to mainstream cardiac imaging modalities for preclinical research, photoacoustic computed tomography (PACT) can provide functional optical contrast with high imaging speed and resolution. However, PACT has not been demonstrated to reveal the dynamics of whole cardiac anatomy or vascular system without surgical procedure (thoracotomy) for tissue penetration. Here, we achieved non-invasive imaging of rat hearts using the recently developed three-dimensional PACT (3D-PACT) platform, demonstrating the regulated illumination and detection schemes to reduce the effects of optical attenuation and acoustic distortion through the chest wall; thereby, enabling unimpeded visualization of the cardiac anatomy and intracardiac hemodynamics following rapidly scanning the heart within 10 s. We further applied 3D-PACT to reveal distinct cardiac structural and functional changes among the healthy, hypertensive, and obese rats, with optical contrast to uncover differences in cardiac chamber size, wall thickness, and hemodynamics. Accordingly, 3D-PACT provides high imaging speed and nonionizing penetration to capture the whole heart for diagnosing the animal models, holding promises for clinical translation to cardiac imaging of human neonates.

Introduction

Cardiovascular disease remains the leading cause of morbidity, responsible for 16% of the world’s total deaths1. Diagnostic imaging and therapeutic targets warrant preclinical investigations of animal models2,3 to understand the cardiac disease mechanisms. Over the last several decades, non-invasive imaging of small animal models provides in vivo insights into the structural and functional phenotypes with physiological and clinical relevance. The current non-invasive imaging modalities allow from a wide range of vertebrate animal models, including light-sheet fluorescent microscopy4, echocardiography5, magnetic resonance imaging (MRI)6, x-ray computed tomography (CT)7, positron emission tomography (PET)8, and/or single-photon emission computed tomography (SPECT)9, and a combination of these imaging techniques10,11,12. Each imaging technique provides the unique degree of tissue penetration, resolution, and contrast to image the specific animal models from zebrafish to mouse to swine models, and the combination of these complementary techniques allows for addressing the spatial and temporal resolution, field of view (FOV), and relative phenotypes in response to the particular imaging needs. For example, echocardiography is a portable tool to assess cardiac function by interrogating the contracting/relaxing heart chambers and open/closure of the valves in real-time. In contrast, CT and MRI are bulky but offer larger FOV and/or finer spatial resolution needed for cardiac anatomy and vascular system7,13.

Complementary to the well-established imaging modalities, photoacoustic (PA) computed tomography (PACT) is an emerging imaging technique that combines functional optical contrast from light illumination and high spatial resolution from acoustic detection. Compared with echocardiography, PACT offers rich optical contrasts that are physiologically relevant14 and suffers a lower level of speckling artifacts15. As compared with MRI, PACT usually provides a higher imaging speed and a more portable imaging platform. Unlike X-ray CT, PACT uses nonionizing illumination to image cardiac vasculature, avoiding injection of contrast agents. Like PET, PACT can reveal molecular information16 with finer spatial resolution obviating the need for ionizing radiation.

Despite these unique advantages, the majority of the published PACT prototypes have not been demonstrated for non-invasive cardiac imaging of the detailed anatomical or functional phenotypes due to suboptimal illumination and detection schemes. This limitation has been recognized by the previous PACT studies, including in vivo imaging of mouse hearts with suboptimal clarity for anatomical/functional details17,18 and ex vivo imaging of excised/perfused heart with no intracardiac flow dynamics19,20. There are 3 challenges for in vivo cardiac PACT. (1) Ribs and lungs surrounding the heart partially block and disturb PA signals. The same limitation exists in echocardiography; however, PACT only suffers from one-way acoustic disturbance since the PA signals are generated in hearts. (2) The high concentration of hemoglobin and myoglobin renders the myocardium highly absorptive to light. This property may produce high PA signal amplitudes; however, it reduces optical penetration into the heart. (3) The periodic heartbeat requires real-time imaging or motion-correction mechanisms to give the imager’s spatial resolution full play for dynamic cardiac structures. For example, the time-gating is widely used in CT21, MRI22, and PET23, dividing every cardiac cycle into multiple phases and collating all the data from a specific phase for imaging reconstruction.

In this context, we seek to overcome these challenges and enhance PACT for pre-clinical and physiological applications. Here, we modified a recently developed three-dimensional PACT (3D-PACT) platform24 with regulated illumination and detection schemes to image the whole rat heart in vivo. This 3D approach allows for a hemispherical detection mesh with an extensive view aperture, detecting PA signals from the heart semi-panoramically to alleviate the acoustic blockage caused by ribs and lungs. In addition, the 1064-nm light suffers less scattering in the biological tissue25 than the wavelengths that are commonly used by PACT within the near-infrared window (680–950 nm). During the scanning, we synchronized the PACT measurement with cardiac cycles by gating with the electrocardiogram (ECG). The 3D-PACT scanned the rat heart in 10 s and reconstructed a series of cardiac images based on a time-gating strategy.

The 3D-PACT along with the ECG-aided synchronization largely resolves the challenges from the previous PACT approaches and reveals the whole rat heart anatomy non-invasively. Specifically, the 3D-PACT significantly enhances the image clarity and captures the dynamic changes in the cardiac structure between obese and control (i.e., lean) rats; namely, chamber size, myocardial wall thickness, and intracardiac flow. Moreover, the 3D-PACT is capable of recording the distinct blood flow in various sizes of the vascular system from healthy, obese, and hypertensive rats, including the aortas, pulmonary arteries, left and right coronary arteries. Notably, the 3D-PACT shows the near-term potential to provide cardiac image quality comparable to that of MRI26, delivering optical contrast for cardiac structure and function, and providing physiological parameters for blood flow.

Results

3D-PACT synchronization with ECG-gating for 4D (3D space + time domain) non-invasive cardiac imaging

In the 3D-PACT platform (Fig. 1a), we applied light illumination with a pulse repetition rate of 50 Hz to the subject through an optical window (Figs. 1a and S1). Four arc-shaped ultrasonic transducer arrays rotated coaxially around the optical window for 90° to form a hemispherical detection matrix (Fig. S2 and Methods section). We covered the imaging aperture with a transparent membrane to separate the imaging system from the subject (rat), which was placed in a prone position over the imaging aperture with chest close to the center (Fig. 1b and Methods section). During the 10-s imaging session, the arrays were scanned for 500 steps (see Methods section), and the heart was periodically beating at about 5 Hz27. Due to the motion artifacts induced by the heart beats, the reconstructed heart image was largely blurred without motion correction (Fig. S3). Therefore, we applied time-gated motion correction guided by the synchronized ECG measurement (see Methods). For example, to reconstruct the heart image when ventricles are at maximum dilation (end of the ventricular diastole), we only used the photoacoustic data acquired at the scanning positions gated to the R waves of the ECG signals (Fig. 1c). Thus, the 3D approach reveals the time-dependent changes in cardiac structure, including the left vs. right ventricle, free wall vs. septum, pulmonary artery vs. aortic arch, and left vs. right carotid arteries.

Fig. 1: 3D-PACT synchronized with ECG for cardiac imaging.

figure 1

a Representative sketch of the 3D-PACT platform with a close-up view of the imaging aperture. b Experimental setup for the rat heart imaging. The rat was mounted in a prone position, and the expanded light beam was directed towards the rat’s chest from the bottom. c Example of the signal synchronization between the 3D-PACT (top) and ECG (bottom) measurements. The 3D-PACT signal was collected from one transducer element on one array at multiple scanning positions (i.e., on the same latitude). In the photoacoustic signal diagram, the horizontal axis represents the scanning steps (i.e., time), and the vertical axis stands for the time-of-flight signal (data acquisition sampling at 20 MHz). The P, R, and T waves in the ECG signals correspond to atrial contraction, ventricular contraction, and ventricular relaxation, respectively. d The rat heart images reconstructed with different ultrasonic apertures, showing the impact of a sizeable ultrasonic detection aperture on image reconstruction. The red arcs on the top right corner of each image represent the cross-sectional view of the ultrasonic detection aperture

Full size image

3D-PACT interrogation of rat cardiac anatomy and vascular system

To address acoustic distortion from the chest ribs and lung air sacs that encompass the heart, we demonstrate a large acoustic detection aperture to collect PA signals from numerous view angles (Fig. 1d). To further improve the PACT reconstruction performance of cardiac anatomy, we synchronized the 3D-PACT with ECG, scanning the heart with ~11 steps (i.e., 11 laser pulses) during each cardiac cycle. Accordingly, we divided every cardiac cycle T into 11 phases and reconstructed a volumetric image for each phase (Fig. 2a and Supplementary Movie 1). Since the image of Phase 1 (i.e., Time 0) was reconstructed using the data acquired at the R waves of the ECG signals, the heart was close to the start of systole. During the 111T−−611T period, the heart was in systole but started to relax from 611T. Such a heartbeat pattern complies with the Wiggers diagram. We further mapped the fluctuations of PA signals across the heart (Fig. S4) to take in the heartbeat as dynamic contrast, showing more spatial movements in the apex of the heart than in the top regions.

Fig. 2: Rat heart anatomy acquired by the 3D-PACT.

figure 2

a Front view of the heart within a cardiac cycle. The heart is identified by a magenta circle at 411T. BA, brachiocephalic artery; ITV, internal thoracic vessels; IV, intercostal vessels. b Cross-sectional images of the heart on the sagittal plane. Each image is a maximum amplitude projection (MAP) of a slice marked by the green dashed lines in (a) (Time 0). c Cross-sectional images of the heart on the coronal plane. Each image is an MAP of a slice marked by the yellow dashed lines in (b) (Time 0). LA, left atrium; LV, left ventricle; PuA, pulmonary artery; RA, right atrium; RV, right ventricle; SVC, superior vena cava. d The same data (at Time 0) shown with color-encoded depths. Shallower structures were peeled away in the lower images to show the posterior anatomy

Full size image

To reveal the internal cardiac structure, we further showed the cross-sectional images on the sagittal (Fig. 2b and Supplementary Movie 2) and coronal planes (Fig. 2c and Supplementary Movie 3). In the cross-sectional images, we can identify the dynamic changes in atria and ventricles during cardiac contraction (systole) and relaxation (diastole). In Fig. 2c, the end ventricular diastolic volume was observed near Time 0, whereas the end atrial diastolic volume was near 611T. We further encoded the imaging depth by colors and peeled away the structures from the chest wall towards the posterior region of the heart, showing cardiac structures at different depths (Fig. 2d). The imaging contrast of the cardiac anatomy acquired by the 3D-PACT is from optical absorption of myoglobin from the myocardium and hemoglobin from the blood28.

3D-PACT interrogation of hypertrophic hearts in the obese rats

It is widely recognized that overweight and obesity are cardiometabolic risk factors for developing acute coronary syndromes and stroke29. Obesity-associated increase in cholesterol levels and vascular inflammation are recognized to cause ventricular remodeling and to reduce blood perfusion to the myocardium30. Here, we demonstrate to capacity of 3D-PACT to interrogate obesity-induced cardiac enlargement or known as hypertrophic abnormalities in cardiac anatomy and function by non-invasively imaging Zucker obese/lean rats of a similar age (Table S1). Compared with the lean rats, whose cardiac cycles typically have 11 phases (11 phases/50 Hz = 0.22 s, i.e., 4.55-Hz cardiac rate), the three obese rats we imaged have 9 phases for each cardiac cycle (i.e., 5.55-Hz cardiac rate), indicating a higher cardiac rate for augmented cardiac output to compensate for the increased metabolic requirements31.

The 3D imaging capability enables the measurement of cardiac anatomy and function. By sectioning the volumetric images on the coronal plane (Fig. 3a, b), we observed concentric hypertrophy32 in the obese rat’s heart with increased free wall thicknesses in both left and right ventricles compared to the lean rat (Fig. 3c). For thickness measurement, we selected the middle slice of the heart on the coronal plane and quantified the full-width-at-half-maximum of the photoacoustic signals across the wall section with a length of about 1.3 mm near the middle (Fig. S5). Ex vivo measurements were performed by splitting the dissected hearts on the coronal plane near the middle after imaging, showing concordance with the in vivo measurements (Fig. S6).

Fig. 3: Differences in cardiac anatomy and function between the Zucker obese and lean rats.

figure 3

a Cardiac anatomy of a Zucker obese rat. Left panel: Color-encoded depth-resolved image of the heart; middle panel: Cross-sectional image of the heart; Right: Longitudinal sectioning of the corresponding hypertrophic heart. LA, left atrium; LV, left ventricle; RA, right atrium; RV, right ventricle. b Cardiac anatomy of a Zucker lean rat (control group). c Free wall thickness measurements of the left ventricle (top: p < 0.01, n = 3) and right ventricle (bottom: p < 0.01, n = 3) between the obese and lean rats. d Relative volume variations of the left ventricles and atriums between the obese and lean rats. The shadow behind curves represents the standard deviation across multiple rats (n = 3)

Full size image

In addition, measurements of the changes in cardiac volumes enabled us to assess both the systolic and diastolic function. We divided each heart image into ten sections with specific thicknesses (1.04 mm) on the coronal plane and segmented ventricular and atrial regions to evaluate the volume variations (Fig. S5). To reduce the individual difference within each group (n = 3), we normalized the volumes of ventricles and atria to their values at Phase 1 (at the end of diastole). The relative changes in the cardiac function between the lean and obese rats (Figs. 3d and S7) were characterized in terms of the end-systolic ventricular volume and end-diastolic atrial volume of the cardiac cycle. Probably due to the reduced strain of ventricles or well-recognized diastolic dysfunction in the hypertrophic hearts32, we observed slightly lower changes in the average of normalized left ventricular volumes in obese rats. Similarly, the left atria in the obese rats developed slightly lower variations in the normalized volume.

3D-PACT interrogation of cardiovascular hemodynamics

Given PACT’s high sensitivity to the hematogenous contrast, the 3D-approach enabled periodic hemodynamic interrogation of the cardiovascular system. In a cardiac image acquired by 3D-PACT (Fig. 4a), the aorta, pulmonary artery, right coronary artery, and left coronary artery were captured. To measure the intravascular hemodynamics, we selected a region of interest (ROI) slightly larger than the vessels’ cross-sections and averaged the photoacoustic signal amplitude within a section of 10-voxel length (1.3 mm) along the vessel. Since oxyhemoglobin dominates the optical absorption at 1064 nm (i.e., μa (oxyhemoglobin) ≈ 10 × μa (deoxyhemoglobin); μa: absorption coefficient)33, the changes in spatially averaged PA signals within the ROI mainly reflect the variations in the amount of oxyhemoglobin, which is determined by the blood volume (or vascular diameter) and oxygen saturation (sO2).

Fig. 4: Cardiac hemodynamics in control (i.e., healthy), hypertensive, and obese hearts.

figure 4

a A heart image acquired by the 3D-PACT shows the aorta, pulmonary artery, right coronary artery (RCA), and left coronary artery (LCA). b PA signal fluctuations in the four cardiac vasculatures among the control, hypertensive, and obese hearts (n = 3 for each group). c Relative changes in ventricular volume of the control (black solid line), hypertensive (blue dashed line), and hypertrophic (red dashed line) hearts during a cardiac cycle (n = 3 for each group)

Full size image

The measurements of hemodynamics in the coronary arteries introduce distinct myocardial perfusion among healthy, hypertensive, and obese rats. The averaged PA signals in the aorta and pulmonary artery, peaked during ventricular systole (0.4–0.5 T), showed slightly higher variations in hypertensive hearts (blue dashed lines in Fig. 4b) than in the control group, yet with insignificant statistical differences. This observation could be supported by measuring relative volume changes of the left ventricles (Fig. 4c). In comparison, we noticed more significant fluctuations in the relative PA signal amplitudes from the left coronary arteries in hypertensive hearts within the cardiac cycle (Fig. S8). We speculate that the abnormal hemodynamics are related to the accumulated fat in the damaged arteries, which may affect blood perfusion34.

Similarly, we measured and plotted the hemodynamics of cardiac vessels in obese rats. Probably due to the augmented cardiac output and obesity-induced arterial hypertension35, we observed higher variations in the averaged PA signals from the four vessels (red dashed lines in Fig. 4b). However, only the right coronary arteries and aortas in the obese rats showed significant statistical differences from the control group (Fig. S8). The combined measurements in cardiac anatomy and hemodynamics depict a potential to provide diagnosis of cardiovascular diseases in individuals with hypertension or obesity.

Discussion

For decades, the photoacoustic field has faced challenges in cardiac imaging, including acoustic disturbance from ribs and lungs, optical attenuation, and motion artifacts from the periodic heartbeat. Nevertheless, preclinical research and clinical practice are in need of innovative imaging methods with complementary advantages to provide high-dimensional physiological information. In this study, we have demonstrated non-invasive 3D-PACT of cardiac anatomy and function in rats, providing substantial improvements in anatomical image clarity and cardiovascular function measurement. To reveal the cardiac anatomy, 3D-PACT employs a hemispherical acoustic detection aperture with large view angles to mitigate the influences induced by ribs and lungs; The 1064-nm light and uniform illumination also facilitate deep penetration of the entire heart (Fig. 2b); although this study has not provided snapshot or real-time imaging, the periodic heartbeat allows the ECG-guided time-gating method to reduce motion-induced artifacts during a 10-second scan. Clear anatomical imaging further enables measurements of the cardiac functions, including the chamber volume variations and cardiovascular hemodynamics. For example, the changes in photoacoustic signals from the cardiac vessels indicate the variations in blood volume and sO2, showing differences in hemodynamics of the hypertensive, hypertrophic, and normal control hearts. In summary, the optical imaging contrast in 3D-PACT depicts the distribution of myoglobin in the myocardium and hemoglobin in the blood, revealing cardiac anatomy and function with physiological relevance36.

While the study was performed based on the modification of a recently published imaging system24, additional technical innovations are critical to producing clear cardiac images: (1) The ECG-guided time-gating strategy synchronized the 3D-PACT system with the ECG measurement by sending a spike signal to one of the ECG’s electrodes when the system starts to run (see Methods). (2) Regulated laser illumination and gated spatial sampling are synchronized to image the heart dynamics during multiple cardiac cycles (see Methods), providing an FOV large enough for whole-heart imaging in 3D space37. (3) In PACT, the light illumination, acoustic detection, and imaging object comprise a comprehensive signal excitation and detection system. Therefore, animal positioning that can reliably match the illumination and detection is usually overlooked but critical for high-quality imaging (see Methods section).

Based on the improvements in this study, PACT of the heart can be further developed towards the following directions: (1) Using ultrasonic arrays with a higher center frequency (e.g., 5–10 MHz) could help reveal detailed cardiac anatomy in small animals with a higher resolution. However, the unwanted effects from acoustic attenuation and rib disturbance will be more significant. (2) To maintain the lateral FOV with an improved imaging speed, rearrangement of the ultrasonic transducers (e.g., eight arc arrays, each with 128 elements) would be implemented to provide an equally dense azimuthal sampling with less scanning steps. (3) The data amplification and acquisition circuits can be modified to integrate ultrasonography, offering well-interpreted imaging contrast (i.e., information) in coregistered images for reference.

3D-PACT is an open imaging platform for preclinical research and clinical translation by simultaneously providing cardiac anatomy and function. Compared with established imaging modalities, PACT of the heart offers additional physiology-related information with high imaging speed and resolution, yet without the need for ionizing radiation or invasiveness, making PACT a promising tool for cardiac imaging of human neonates.

Materials and methods

Optical illumination and acoustic detection

For optical illumination, PACT requires a relatively uniform energy deposition within the field of view. While for acoustic detection, high-quality imaging requires a large view angle (i.e., synthetic aperture) and dense sampling in time and space. In 3D-PACT, we directed a 1,064 nm laser beam (DLS9050, Continuum; 50 Hz pulse repetition rate; 5–9 ns pulse width) to the optical window mounted on the scanning axis. The optical window is composed of an engineered diffuser (EDC-15, RPC Photonics Inc.) and a condenser lens (ACL25416U-B, Thorlabs Inc.), expanding the beam to a diameter of ~4.8 cm (Fig. S2). The illumination area is slightly larger than the imaging object (i.e., heart) to lower the requirement of precise animal positioning. Because the heart is filled with blood (i.e., high optical absorption), we used 1064-nm light in this study for less optical attenuation in biological tissues to achieve deep penetration. The optical fluence on the tissue surface (~20 mJ/cm2 at 50 Hz) was limited by the American National Standards Institutes’ safety standards38.

The acoustic detection module comprises four ultrasonic arc arrays with a separation of 90 degrees. Each array has 256 transducer elements with a central frequency of 2.25 MHz, and each element has a dimension of 0.6 × 0.7 mm2, generating a divergence angle (i.e., far-field acoustic diffraction) around 60 degrees covering the imaging object semi-panoramically. Here, we applied a hemispherical detection matrix to enlarge the view angle so the system can acquire acoustic signals propagating in most independent directions. In addition, the 2.25-MHz center frequency of the array also allows the detection of sufficiently low frequency signals that suffer less distortion from the ribs. Notably, this acoustic detection scheme would not avoid distortions from bones and lungs, but would reduce their effects on the reconstructed images. Connecting to the ultrasonic arrays, 1024-channel pre-amplification and analog-to-digital conversion circuits (PhotoSound, Inc.) provide one-to-one mapped data acquisition24.

ECG-guided time gating and spatial sampling

The data acquisition modules and mechanical scanner of the 3D-PACT system were triggered by the laser’s “Lamp Sync Out” trigger output. At the same time, the computer controlled the trigger signal transmission via an electronic relay in the system’s control box. We turned on the laser and ECG systems shortly before closing the relay, which coordinated the data acquisition, laser pulses, and mechanical scanning. We further directed a coaxial cable from the system’s control box to one of the ECG probes (left or right arm) and output a spike signal much higher than the ECG spikes when the relay was closed (Fig. S9). Accordingly, we pinpointed the start time of the photoacoustic data acquisition in the ECG signal, achieving the data synchronization between ECG and 3D-PACT.

We then grouped the PACT data into multiple heartbeat phases based on the synchronized ECG cardiac cycle and reconstructed a volumetric heart image for each phase (Fig. 1c). Since the heart rate was around 5 Hz and the pulsed laser illuminated at 50 Hz, each cardiac cycle can be divided into 9–11 phases, enabling the reconstruction of 9–11 volumetric images to depict the beating heart. The variations in cardiac cycle duration were caused by differences in weight (e.g., obese versus regular rats) and anesthesia (i.e., light versus deep anesthesia). However, the ECG signal shows that the heart rate was relatively stable during the 10-second scan.

For spatial sampling, each of the four arrays was scanned for 500 steps (50 Hz × 10 s) during the 10-second imaging, providing 2000 azimuthal sampling positions. Accordingly, each heart image is reconstructed by ~200 (2000/10 phases) azimuthal positions times 256 polar positions (i.e., each arc array has 256 elements), which could provide an effective FOV with a lateral diameter of 15 mm37 and elevational radius of 38 mm.

Animal preparation

Four strains of rats were imaged in this study (Table S1): (1) Hsd:Sprague Dawley SD rats (120–150 g body weight), (2) Zucker obese rats (200–250 g), (3) Zucker lean rats (120–150 g), (4) Spontaneously hypertensive inbred SHR/NHsd rats (120–150 g). We first imaged the SD rats to optimize the experimental setup and demonstrate the imaging capability (Fig. 2). The Zucker lean rats were imaged as the control group to compare the cardiac anatomy and function with the obese rats (Fig. 3). We compared the cardiovascular hemodynamics in the hypertensive, obese, and SD rats (Fig. 4).

The imaging aperture was covered by a transparent membrane (disposable plastic wrap) that separated the rat from the imaging system and served as a support for the rat. We then covered a feedback-controlled heating pad with a zip bag for waterproofing and placed it above the membrane to maintain the animal’s body temperature (Fig. 1b). Before imaging, we removed the hair on the chest and placed the rat above the heating pad. A 3D-printed nose cone and tooth bar were positioned to provide air with 2% vaporized isoflurane. During imaging, we placed two cross lines above the imaging aperture to mark the center position of the FOV. We further bent the rat’s front limbs towards the back of the body so that the heart would be prominent towards the chest wall. After the imaging session was completed, we euthanized the rat, collected the heart, and obtained a photo of the dissected heart for reference. All the animal experiments followed the protocol approved by the Institutional Animal Care and Use Committee (IACUC) of the California Institute of Technology.

Image processing and visualization

We used the dual-speed-of-sound back-projection algorithm39 implemented in C + + to reconstruct all images with mitigated artifacts induced by the acoustic inhomogeneity between the biological tissue and the coupling fluid (D2O). Each volumetric image was reconstructed with a voxel size of 0.13 × 0.13 × 0.13 mm3. We then processed each image to improve the contrast through the following steps: (1) A depth compensation (e0.81×depth (cm)) method was applied to enhance the PA amplitude in deep tissue. (2) We then used a high-pass filter to suppress the low-frequency background. (3) The high-passed image was further denoised using sparse 4D transform-domain collaborative filtering40. (4) To enhance the contrast of blood vessels, a Hessian-based Frangi vesselness filter41 was applied to the denoised image. (5) Finally, we added the vesselness-filtered images (self-normalized) with a weighting factor of 0.5 back to the high-passed images (with a weighting factor of 0.5) and obtained the presented images. The color-encoded images in Figs. 2 and 3 were acquired using a specialized 3D visualization software package, namely 3D PA Visualization Studio42.

Data availability

All data are available within the Article and Supplementary Files, or available from the authors upon request.

References

1. Word Health Organization. The Top 10 Causes Of Death. https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death (2020).

2. Neeland, I. J. et al. Effects of liraglutide on visceral and ectopic fat in adults with overweight and obesity at high cardiovascular risk: a randomised, double-blind, placebo-controlled, clinical trial. Lancet Diabetes Endocrinol. 9, 595–605 (2021).

Article Google Scholar

3. Lobelo, F. et al. Routine assessment and promotion of physical activity in healthcare settings: a scientific statement from the American Heart Association. Circulation 137, e495–e522 (2018).

Article Google Scholar

4. Fei, P. et al. Cardiac light-sheet fluorescent microscopy for multi-scale and rapid imaging of architecture and function. Sci. Rep. 6, 22489 (2016).

Article ADS Google Scholar

5. Chetboul, V. Advanced techniques in echocardiography in small animals. Vet. Clin. North Am.: Small Anim. Pract. 40, 529–543 (2010).

Article Google Scholar

6. Nahrendorf, M. et al. Cardiac magnetic resonance imaging in small animal models of human heart failure. Med. Image Anal. 7, 369–375 (2003).

Article Google Scholar

7. Bertolini, G. & Angeloni, L. Vascular and cardiac CT in small animals. in Computed Tomography-Advanced Applications (ed Halefoğlu, A. M.) (Rijeka: IntechOpen, 2017).

8. Wang, Q. et al. Application of animal and human PET in cardiac research. Am. J. Cardiovasc. Dis. 8, 24–30 (2018).

Google Scholar

9. Golestani, R. et al. Small-animal SPECT and SPECT/CT: application in cardiovascular research. Eur. J. Nucl. Med. Mol. Imaging 37, 1766–1777 (2010).

Article Google Scholar

10. Koba, W., Jelicks, L. A. & Fine, E. J. MicroPET/SPECT/CT imaging of small animal models of disease. Am. J. Pathol. 182, 319–324 (2013).

Article Google Scholar

11. Nanni, C. & Torigian, D. A. Applications of small animal imaging with PET, PET/CT, and PET/MR imaging. PET Clin. 3, 243–250 (2008).

Article Google Scholar

12. Izaguirre, E. W. et al. Dual modality micro-SPECT and micro-CT for small animal imaging: technical advances and challenges. Proceedings of the SPIE 5923, Penetrating Radiation Systems and Applications VII. (SPIE, 2005).

13. Botnar, R. M. & Makowski, M. R. Cardiovascular magnetic resonance imaging in small animals. Prog. Mol. Biol. Transl. Sci. 105, 227–261 (2012).

Article Google Scholar

14. Yao, J. J., Xia, J. & Wang, L. V. Multiscale functional and molecular photoacoustic tomography. Ultrasonic Imaging 38, 44–62 (2016).

Article Google Scholar

15. Guo, Z. J., Li, L. & Wang, L. V. On the speckle-free nature of photoacoustic tomography. Med. Phys. 36, 4084–4088 (2009).

Article Google Scholar

16. Weber, J., Beard, P. C. & Bohndiek, S. E. Contrast agents for molecular photoacoustic imaging. Nat. Methods 13, 639–650 (2016).

Article Google Scholar

17. Deán-Ben, X. L. et al. Spiral volumetric optoacoustic tomography visualizes multi-scale dynamics in mice. Light Sci. Appl. 6, e16247 (2017).

Article Google Scholar

18. Ivankovic, I. et al. Volumetric optoacoustic tomography enables non-invasive in vivo characterization of impaired heart function in hypoxic conditions. Sci. Rep. 9, 8369 (2019).

Article ADS Google Scholar

19. Lin, H. C. A. et al. Ultrafast volumetric optoacoustic imaging of whole isolated beating mouse heart. Sci. Rep. 8, 14132 (2018).

Article ADS Google Scholar

20. Özsoy, Ç. et al. Ultrafast four-dimensional imaging of cardiac mechanical wave propagation with sparse optoacoustic sensing. Proc. Natl Acad. Sci. USA 118, e2103979118 (2021).

Article Google Scholar

21. Desjardins, B. & Kazerooni, E. A. ECG-gated cardiac CT. Am. J. Roentgenol. 182, 993–1010 (2004).

Article Google Scholar

22. Nacif, M. S. et al. Cardiac magnetic resonance imaging and its electrocardiographs (ECG): tips and tricks. Int. J. Cardiovasc. Imaging 28, 1465–1475 (2012).

Article Google Scholar

23. Lassen, M. L., Kwiecinski, J. & Slomka, P. J. Gating approaches in cardiac PET imaging. Pet. Clin. 14, 271–279 (2019).

Article Google Scholar

24. Lin, L. et al. High-speed three-dimensional photoacoustic computed tomography for preclinical research and clinical translation. Nat. Commun. 12, 882 (2021).

Article ADS Google Scholar

25. Jacques, S. L. Optical properties of biological tissues: a review. Phys. Med. Biol. 58, R37–R61 (2013).

Article ADS Google Scholar

26. Price, A. N. et al. Cardiovascular magnetic resonance imaging in experimental models. Open Cardiovasc. Med. J. 4, 278–292 (2010).

Article MathSciNet Google Scholar

27. Azar, T., Sharp, J. & Lawson, D. Heart rates of male and female sprague-dawley and spontaneously hypertensive rats housed singly or in groups. J. Am. Assoc. Lab. Anim. Sci. 50, 175–184 (2011).

Google Scholar

28. Lin, L. et al. In vivo photoacoustic tomography of myoglobin oxygen saturation. J. Biomed. Opt. 21, 061002 (2015).

Article Google Scholar

29. Piché, M. E., Tchernof, A. & Després, J. P. Obesity phenotypes, diabetes, and cardiovascular diseases. Circ. Res. 126, 1477–1500 (2020).

Article Google Scholar

30. Clark, L. T. Cholesterol and heart disease: current concepts in pathogenesis and treatment. J. Natl Med. Assoc. 78, 743–751 (1986).

Google Scholar

31. Alpert, M. A. & Hashimi, M. W. Obesity and the heart. Am. J. Med. Sci. 306, 117–123 (1993).

Article Google Scholar

32. Aurigemma, G. P., De Simone, G. & Fitzgibbons, T. P. Cardiac remodeling in obesity. Circ.: Cardiovasc. Imag. 6, 142–152 (2013).

Google Scholar

33. Pantaleone, C. Cancer Immunotherapy and Biological Cancer Treatments (ed. Arnouk, H.) (IntechOpen, 2019).

34. Mayet, J. & Hughes, A. Cardiac and vascular pathophysiology in hypertension. Heart 89, 1104–1109 (2003).

Article Google Scholar

35. Artham, S. M. et al. Obesity and hypertension, heart failure, and coronary heart disease-risk factor, paradox, and recommendations for weight loss. Ochsner J. 9, 124–132 (2009).

Google Scholar

36. Hendgen-Cotta, U. B., Kelm, M. & Rassaf, T. Myoglobin functions in the heart. Free Radic. Biol. Med. 73, 252–259 (2014).

Article Google Scholar

37. Hu, P. et al. Spatiotemporal antialiasing in photoacoustic computed tomography. IEEE Trans. Med. Imaging 39, 3535–3547 (2020).

Article Google Scholar

38. American National Standards Institute. American National Standard for the Safe Use of Lasers ANSI, z136. 131-2014 (Laser Institute of America, 2014).

39. Li, L. et al. Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution. Nat. Biomed. Eng. 1, 0071 (2017).

Article Google Scholar

40. Maggioni, M. et al. Nonlocal transform-domain filter for volumetric data denoising and reconstruction. IEEE Trans. Image Process. 22, 119–133 (2013).

Article ADS MathSciNet MATH Google Scholar

41. Frangi, A. F. et al. Multiscale vessel enhancement filtering. Proceedings of the First International Conference on Medical Image Computing and Computer-Assisted Intervention. (pp. 130–137. (Springer, Cambridge, 1998).

Google Scholar

42. Cho, S. et al. 3D PHOVIS: 3D photoacoustic visualization studio. Photoacoustics 18, 100168 (2020).

Article Google Scholar

Download references

Acknowledgements

This work was sponsored by the United States National Institutes of Health (NIH) grants R35 CA220436 (Outstanding Investigator Award) and U01 NS099717 (BRAIN Initiative).

Author information

1. Author notes

Li Lin

Present address: College of Biomedical Engineering and Instrument Science, Zhejiang University, Hangzhou, China

2. Li Lin

Present address: The First Affiliated Hospital, Zhejiang University School of Medicine, Hangzhou, China

3. These authors contributed equally: Li Lin, Xin Tong

Authors and Affiliations

1. Caltech Optical Imaging Laboratory, Andrew and Peggy Cherng Department of Medical Engineering, Department of Electrical Engineering, California Institute of Technology, Pasadena, CA, USA

Li Lin, Xin Tong, Yide Zhang, Shuai Na, Rui Cao & Lihong V. Wang

2. Department of Bioengineering, UCLA, Los Angeles, CA, USA

Susana Cavallero & Tzung K. Hsiai

3. Division of Cardiology, Department of Medicine, UCLA, Los Angeles, CA, USA

Tzung K. Hsiai

Contributions

L.V.W. and L.L. conceived and designed the study. L.L. constructed the hardware system. S.N. and L.L. developed the control program. Y.Z. and R.C. modified the control program for ECG synchronization. L.L. and X.T. performed the experiments. L.L., X.T., and Y.Z. analyzed the data. L.V.W., T.K.H., and S.C. supervised the study. All authors wrote the manuscript.

Corresponding authors

Correspondence to Tzung K. Hsiai or Lihong V. Wang.

Ethics declarations

Conflict of interest

All the authors declare no competing interests. L.V.W. has a financial interest in Microphotoacoustics, Inc., CalPACT, LLC, and Union Photoacoustic Technologies, Ltd., which, however, did not support this work.