Introduction

Systemic imaging stands as a pivotal tool in advancing life sciences, offering various modes for comprehensive development1,2. Non-invasive three-dimensional (3D) visualization of intrinsic physiology and dynamics allows a macroscopic understanding of organ-level biological mechanisms. However, in-vivo systemic imaging grapples with challenging compromises involving spatiotemporal resolution, limited field of view (FOV) or penetration depth, and acceptable image quality amidst the presence of abundant contrast sources. Presently, the widespread use of magnetic resonance imaging (MRI) allows researchers to obtain quality 3D whole-body images, albeit with construction costs and strong magnetic fields over minutes due to physical limitations3,4. X-ray computed tomography (X-CT) lacks soft tissue contrast and is unsuitable for dynamic imaging due to low temporal resolution and prolonged radiation exposure5. Positron emission tomography (PET) specializes in invasive diagnoses of cardiac-cerebral vascular diseases and cancer screening but suffers from poor spatial resolution6. Ultrasonic imaging (USI) non-invasively provides structural images but lacks endogenous molecular contrasts, impairing functional imaging capabilities7,8,9. These deficiencies hinder deeper insights into biological studies. Although optical imaging (OI) offers multifarious functional visualizations, conventional OI faces challenges caused by strong optical scattering in biological tissue, limiting systemic imaging quality beyond the optical scattering limit of 1–2 mm in depth10. Overcoming these challenges is crucial for advancing systemic imaging in life sciences.

Photoacoustic (PA) Computed Tomography (PACT), as an emerging imaging modality, stands out for its distinctive ability to visualize optical absorption in deep tissues11. Utilizing light pulse illumination and acoustic detection, PA signals carrying information on light absorption propagate from tissues several centimeters deep, creating images with acoustic resolution and optical contrast12. This approach holds great promise for whole-body imaging of small vertebrates, with various PACT systems demonstrating their capacity to generate high-fidelity structural images of cross-sections of the trunk or brain13,14,15. The unique combination of high-velocity and multiwavelength imaging abilities in two-dimensional (2D) PACT previously extended its performance to the tracking of metabolism dynamics16, monitoring physiological parameters, and recording exogenous probes for specialized applications17,18,19. For example, existing studies have already been able to achieve dynamic assessment of metabolic functions in specific organs using 2D-PACT20. Despite these advancements, PACT systems designed for 3D imaging face challenges, including low spatiotemporal resolution, confined FOV, and image quality without optimization. In one approach, a PACT system equipped with 128 ultrasonic transducers integrated into a hemispherical surface was introduced21, however, artifacts persist outside the FOV due to angular sparse sampling. To address these challenges, some researchers scanned a dome-type array with 256 elements in a helical line, successfully obtaining 3D images of the trunk with clear vessels and recording the pharmacokinetics in the relatively shallow tissue21. As a further optimization, other researchers adopted the 512-element dome-type array with higher center frequency for raster and angular scanning to achieve whole-trunk imaging22. However, the application in functional imaging of this system was limited by slow imaging speed and shallow penetration depth. In addition, from the perspective of in-vivo imaging quality, the increase of the center frequency of the transducers did not appear to achieve higher spatial resolution as expected. Another significant development involved a 3D PACT system with superior performance achieved by rotating an arc-shaped transducer array at high speed, allowing for a large FOV and a relatively high frame rate22. While this design was successfully demonstrated in rat brain and human breast imaging, the system’s molecular imaging capability remains unproven, necessitating further exploration and validation through multispectral unmixing techniques23,24. These advancements underscore the evolving landscape of PACT, showcasing its potential for transformative contributions to biomedical imaging. Overall, the inherent trade-off between the FOV and spatiotemporal resolution limits the application of PACT in biomedical research, particularly in studying the dynamics of multiple organs and cross-regional systems throughout the body.

In this study, we present a state-of-the-art 3D panoramic photoacoustic computed tomography (3D-PanoPACT), representing a cutting-edge advancement in PACT technology. Leveraging a high-density transducer array arranged according to the Fibonacci grid and rigorous engineering design and calibration, 3D-PanoPACT achieves high-fidelity, high-spatiotemporal-resolution imaging within scalable FOVs, allowing analysis of physiological dynamics from whole-organ to whole-body level. We used 3D-PanoPACT to clearly image the whole liver at a frame rate of 25 Hz, extracted the heartbeat and respiratory signals, and mapped the 3D arterial network and the phase gradient of the liver lobe related to the pulse wave, employing a single wavelength. We further obtained the authentic whole-brain vascular anatomy, with particular emphasis on the vivid depiction of pivotal basicranial arterial structures, the Circle of Willis. Based on this, whole-brain functional dynamics, from the Circle of Willis to the cortex, can be recorded and analyzed at a frame rate of up to 10 Hz. To underscore the 3D-PanoPACT’s imaging capacity at whole-body scale, we introduced 10 Hz dynamic reconstruction with high quality of the whole trunk, from the thoracic cavity to the reproductive system, visually presenting the 3D anatomy with elaborate vascular networks as well as rhythmic organ motions in nature, with a penetration depth beyond 20 mm. We further tracked the metabolic pathways of a small-molecule probe across multiple organs with high spatiotemporal resolution using 3D-PanoPACT, thereby demonstrating its significant practical value in cross-regional dynamic studies at the whole-body scale. Therefore, we illustrated the comprehensive imaging performance of 3D-PanoPACT by extending from single-organ to whole-body FOV and from structural to functional imaging. Given that high-speed 3D imaging of the whole brain and whole body are significant research areas, presenting the whole-brain and whole-body imaging results separately helps to highlight the broad applicability of 3D-PanoPACT in addressing diverse imaging needs and its potential for solving key biological questions. For a clear illustration, we further provided a detailed description of the implementation and effectiveness of the spatiotemporal-integration (STINT) method, which is extensively used in this study. Consequently, 3D-PanoPACT not only explores but fully harnesses the superiority of the PA modality, attracting increasing attention from biomedical practitioners.

Results

Development and calibration of 3D-PanoPACT

The constructed 3D-PanoPACT system is illustrated in Fig. 1a. The foundational framework of the imaging system comprises four modules: (1) 1024 piezoelectric ultrasonic transducers are integrated into a hemispherical housing with a diameter of 20 cm (Fig. 1a). Each transducer element is non-focusing with a central frequency of 3.16 MHz and a one-way bandwidth of 60% (Fig. S1, see Supplementary Information for details). The transducer elements are arranged according to the Fibonacci grid25 to form a high-density 3D array with uniform acoustic detection field. Considering the mechanical geometrical deviations, the actual position of each element was calibrated so as to affect the image reconstruction quality (see “Methods”). (2) Four sets of encapsulated 256-channel parallel analog-digital conversion modules are attached to four connectors of the ultrasonic array (Fig. 1a). Each module contains 256-channel pre-amplification circuits and 256-channel data acquisition (DAQ) circuits interlinked with an adapter plate for facilitating interface correspondence. This direct plug-in architecture is conducive to minimizing coupling noise. (3) The illumination light from the lasers is delivered through spatial light paths. One of the Nd:YAG lasers provides fundamental 1064 nm output and the other one equipped with an optical parametric oscillator (OPO) generates 670-900 nm tunable output. The 3D-PanoPACT offers two illumination modes. 1064 nm output reaches a maximum laser pulse repetition rate of 25 Hz when single-wavelength imaging mode is applied (Fig. 1c). The tunable and 1064 nm outputs work in pairs to provide a 10 Hz dual-wavelength imaging mode, whose pulse repetition rate is limited by the damage threshold of the crystal inside OPO (Fig. 1c). The laser pulse interval for each group was set to 20 µs, preventing aliasing of the back-and-forth arriving PA wavefronts. To achieve a 20 μs delay in the output of laser 2 (1064 nm output) relative to laser 1 (tunable output), we delayed the flash signal of laser 1 by 50 μs relative to that of laser 2, resulting in a laser output interval of 160 + 50 − 230 = −20 μs, corresponding to a 20 μs delay. Two outputs are beam-combined and pass through an engineered diffuser sealed on the bottom of the hemispherical housing to obtain a compatible illumination inside the imaging region (Fig. 1a). (4) The customized animal holders accommodate the systemic imaging tasks. For frontal and backside trunk imaging, the anesthetized animals are secured to the holder via limb restraints, and the holder is mounted onto the triaxial translation stage, allowing precise movement of the entire body in space. The animal’s trunk is submerged in deionized water, with its head above the surface, connected to a gas anesthesia apparatus to maintain respiration (Fig. 1d). For brain imaging, the animal is placed with the brain region adjusted to the center of the array and, the interface between the mask and nose is sealed by parafilm (Fig. 1e).

Fig. 1: System design and imaging examples of 3D-PanoPACT.

a The schematic diagram of the 3D-PanoPACT system. HUTA hemispherical ultrasonic transducer array, AMP amplification circuits, DAQ data acquisition module, RS rotation stage, BC beam combiner, FSM front-silvered mirror, ED engineered diffuser, OPO Optical Parametric Oscillator, TS translation stage. b Real-time 3D imaging of the whole liver anatomy produced by the 3D-PanoPACT with separate displays of the different liver lobes attached. The image was acquired with a single laser pulse. MLL middle liver lobe, LMLL left middle liver lobe, LLLL left lateral liver lobe. c The sequence diagram of 25 Hz single-wavelength mode and 10 Hz dual-wavelength mode. d The schematic diagrams of in-vivo trunk imaging in 3D-PanoPACT. e The schematic diagrams of in-vivo brain imaging in 3D-PanoPACT.

Considering that acoustic detection is crucial for imaging quality, the determination of the array element size has undergone special design and validation. We first used the k-Wave toolbox to investigate the relationship between detection sensitivity and transducer element size. We placed an ideal PA point source at the center of FOV with the transducer set to the center frequency and bandwidth used and observed the relationship between the detected PA amplitude and the radius of the transducer element (Fig. S2a). The simulation results indicate that the PA amplitude exhibits a quadratic relationship with the radius (i.e., is linearly proportional to the area) when the radius is small. However, when the radius is large, the simulated values fall below the ideal values, and the quadratic relationship is no longer maintained (Fig. S2b). The critical point where these two relationships intersect is around Radius = 2.5 mm. In other words, designing the transducer element radius to be 2.5 mm can maximize the detection sensitivity at the highest cost-effectiveness ratio. We further used Field II to simulate the receiving aperture angles for different element radii. The simulations show that while a transducer with a radius of 0.5 mm has a large aperture angle (Fig. S2c), its sensitivity is 25 times lower than that of a transducer with a 2.5 mm radius (Fig. S2b), which is not conducive to deep tissue imaging. For comparison, the transducer with a radius of 1.5 mm has an aperture angle comparable to that of the transducer with a 2.5 mm radius, but its sensitivity is 2.78 times weaker. After evaluating both detection sensitivity and receiving aperture angle, we determine that a transducer element radius of 2.5 mm (i.e., 5 mm diameter) offers the optimal balance for 3D-PanoPACT.

In 3D-PanoPACT, the deviation of the detection coordinates is influenced by both the element position deviations of the hemispherical array itself and the deviations of the rotation axis, which can affect the quality of image reconstruction26,27. To calibrate the spatial position of each element, we used a point absorber (50 μm diameter black polystyrene microspheres) embedded in an agarose phantom. The array housing, machined with high precision for an exact 20 cm diameter and ideal Fibonacci grid distribution, still introduced radial deviations due to individual element insertion, necessitating calibration. The calibration procedure involved three steps (see “Methods” for details). First, we centered the FOV at the origin and adjusted the point absorber’s position while repeatedly imaging it at 690 nm until it was centered in the field with PA raw data acquired (Fig. S3a). Second, we measured the water temperature with a high-precision thermocouple and calculated the speed of sound in the medium using the classic temperature-velocity formula28. Third, we positioned the signals from the transducer surface and the point absorber in the raw data, identified the sample points where the PA signals peaked (Fig. S3b), and calculated the radial deviation for each element using the travel time and acoustic velocity29 (see “Methods” for details). The corrected positions of the elements were then determined in Cartesian coordinate with the radial deviations kept within ±0.5 mm (Fig. S3c). We then calibrated the system rotation axis due to assembly errors. To correct this, we positioned a point absorber 15 mm from the origin and imaged it while rotating the array by 70°. By co-reconstructing 15 evenly-spaced frames, we formed a 3D image where the center of the reconstructed points clearly defined a circle (Fig. S3d). The normal through the center was identified as the calibrated rotation axis of the system, allowing us to correct the misalignment and improve image quality. The calibration revealed that the actual rotation axis was parallel to the ideal axis but shifted by 0.07 mm along the X-axis and 0.29 mm along the Y-axis (Fig. S3d). Applying these calibrations, we observed improvements in resolution and contrast, particularly in regions away from the center of the FOV (Fig. S3e, f).

To enhance the quantitative accuracy of the PA images, we addressed the issue of limited bandwidth in the transducers and subsequent circuits, which affects the fidelity of the acquired signals (see “Methods”). We employed the edge-emission method to measure the system’s electrical impulse response (EIR)30,31, using a 10 μm thick polyimide film cut into a 2 × 15 mm rectangle as the PA signal source (Fig. S4a, b). This method produced a pair of negative and positive unipolar PA pulses at the edges of the film, which were used to determine the system EIR (Fig. S4c). To demonstrate the effectiveness of EIR deconvolution, we imaged three square polyimide films with side lengths of 3 mm, 6 mm, and 9 mm. The deconvolution process significantly improved the reconstruction of low-frequency components (Fig. S4d, e), which are crucial for accurate unmixing of oxyhemoglobin saturation (sO2) and visualizing neurovascular coupling responses in the cortex.

In PA imaging, accurate optical fluence compensation is crucial for enhancing image fidelity and ensuring reliable quantitative analysis32, especially for multi-wavelength unmixing. We employed empirical testing to estimate optical attenuation in freshly excised tissue for subsequent compensation. A simple setup was first constructed using two quartz plates, a collimated laser beam, a 1:9 beam splitter, and a pair of optical energy meters. One optical energy meter was used to detect the light passing through the quartz plates, while the other was used to monitor the energy fluctuations of the split beam as a normalization reference (Fig. S5). By comparing these measurements with and without a sample, we calculated the optical attenuation coefficient (({\mu }_{{\mbox{eff}}})) for brain and liver tissues at different wavelengths based on the Lambert-Beer law33, which relates optical attenuation to the properties of the material (see “Methods” for details). For brain tissue, ({\mu }_{{\mbox{eff}}}) was 0.065 mm−1 at 1064 nm and 0.077 mm−1 at 800 nm, which were utilized for imaging whole-brain anatomies and dynamic functions. For liver tissue, ({\mu }_{{\mbox{eff}}}) was 0.075 mm−1 at 1064 nm and 0.133 mm−1 at 690 nm, which were applied in dynamic liver imaging and whole-body imaging. These coefficients were then used in image processing to apply exponential depth compensation, enhancing the PA amplitude from deep tissues.

The spatial resolution of the 3D-PanoPACT system was tested with a 10-µm-diameter tungsten wire, which was imaged at various radial positions along the Z-axis (0 mm, 5 mm, 10 mm, and 15 mm) to quantify the variation in spatial resolution across different regions in the FOV (Fig. S6a–d). Gaussian fitting of the profiles of the X-Z and X-Y sections was performed, and the full width at half maximum (FWHM) values were used to characterize lateral and axial resolutions (Fig. S6e, f). It was observed that, as moving radially from the center, the lateral resolution degraded from 178.6 to 312.3 μm, while the axial resolution remained stable at approximately 250 μm (Fig. S6g). This phenomenon was attributed to the angular sampling in the transverse plane becoming sparser with increasing radial distance. On the other hand, within a radial distance of 15 mm (corresponding to an FOV of 30 mm diameter), the resolution could be regarded as near-isotropic, thereby benefitting the imaging quality across a broad FOV. In addition, the high-density acoustic detection enables an ample volumetric imaging field captured by each single laser excitation, empowering the 3D-PanoPACT with the ability to perform real-time volumetric imaging for living animals. By utilizing the weak scattering properties of 1064 nm laser output in tissues, 3D-PanoPACT generated high-fidelity whole-liver vascular network at 25 Hz within the effective FOV size of ~6 cm3 and an imaging depth of ~15 mm (Fig. 1b). For ease of observation, individual lobes of the liver are separately displayed, showcasing rich vascular morphology and details (Fig. 1b). The FOV and the imaging depth of the 3D-PanoPACT are further extended through the spatiotemporal-integration (STINT) method, enabling dynamic visualization of deep tissue ranging from organ to whole-body scales with high fidelity.

Real-time visualization of liver dynamics in 3D-PanoPACT

The 3D-PanoPACT possesses considerable single-pulse 3D imaging capabilities within a certain FOV. Its imaging frame rate depends on the laser repetition rate when single-pulse imaging is applied, enabling real-time 3D imaging. Here, we demonstrate that 3D-PanoPACT was capable of achieving real-time 3D imaging of the liver non-invasively, enabling the study of rapid hemodynamic processes at the organ level, without exogenous contrast agents (Movie S1, see Supplementary Information for details). At an imaging frame rate of 25 Hz (the 1064 nm laser repetition rate was set to 25 Hz), respiratory motions and heartbeats were fully captured under 1064 nm excitation at well above the Nyquist sampling rate. By placing the left lateral lobe, which is the largest lobe of the liver, at the center of the imaging FOV of the hemispherical array, obvious displacement can be observed during the expansion and contraction of the thoracic cavity (Fig. 2a). By recording the signal changes inside the selected voxel of the vessel (the blue arrow in Fig. 2a), in addition to the low-frequency signal changes caused by respiratory motion, the high-frequency heartbeat signal was also captured by the 3D-PanoPACT (Fig. 2b). Fourier analysis showed that two main frequency components existed in the recorded signal of the selected voxel, corresponding to the heartbeat frequency (around 3.9 Hz) and respiratory frequency (around 0.8 Hz), respectively. By applying a finite-impulse-response high-pass filter on the recorded signal, the heartbeat component can be extracted while keeping the phase constant, which was shown as a quasi-sine wave (dash line in Fig. 2b).

Fig. 2: Real-time visualization of liver dynamics in 3D-PanoPACT.

a The 1064 nm real-time 3D imaging of the liver at 25 Hz. The images are taken at different points in a respiratory cycle, which was labeled in (b) with gray lines. b The raw signal of PA amplitude change (dark blue solid line) of the blood vessel (blue arrow in a) and the extracted heartbeat signal (light blue dashed line). c The Fourier transform of the raw data showing the respiratory frequency and heartbeat frequency, respectively. d The comparison of heartbeat signals extracted from blood vessels at different locations (the blue and orange arrows in a) with the phase delay clearly visible. The corresponding time windows are indicated by shaded areas in (b). e Heartbeat-encoded arterial network mapping overlaid on the anatomy background in gray scale. f Mapping of pulse wave phase delay overlaid on the gray-scale anatomy background. SV superficial blood vessels, HPV hepatic portal vein, HA hepatic artery, MA mesenteric artery.

Additionally, we chose another vessel voxel that was far away from the previous location (the orange arrow in Fig. 2a) and performed a similar analysis on the recorded signal. By comparing the heartbeat component waves obtained from the two locations, a stable phase delay was revealed in the zoomed-in panel (Fig. 2d). In one heartbeat cycle, the contraction of the ventricles leads to the rapid influx of blood into the aorta, causing the vessel walls to contract, subsequently generating pressure waves that propagate throughout the arterial network, which may be an explanation for the phase delay between two locations. With the imaging sequence, voxel-wise extraction of the amplitude at heartbeat frequency was performed (see “Methods”) and the arterial network can be mapped on the 3D liver anatomy (Fig. 2e). If equipped with a higher-repetition-rate pulsed laser, the 3D-PanoPACT could potentially capture cardiac substructure signals and extract information using more complex filters, as these signals often require a higher sampling frame rate for accurate analysis34. It is intuitively seen that, although the superficial vessels have stronger amplitude, they cannot be identified as arteries by this method, which is consistent with the fact that the superficial vessels are almost entirely composed of venous branches. On the contrary, the hepatic artery and the mesenteric artery are highlighted by heartbeat encoding. Further, by extracting the phase of each voxel’s time sequences at the cardiac frequency component, we are able to map the 3D phase gradient of hepatic vasculature (Fig. 2f) (see “Methods”), which may be related to the spatial gradient of blood pressure35. This unique mapping relies on the high-quality video-rate 3D imaging enabled by 3D-PanoPACT, which may serve as a reference for assessing overall hepatic hemodynamics or as an indicator for certain liver diseases, such as acute liver injury36,37. Benefiting from the imaging capability and high spatiotemporal resolution of 3D-PanoPACT, we employed single-wavelength imaging to non-invasively map the 3D arterial networks and quantify the relative phase difference of pulse waves between arteries within an FOV covering the entire organ. This example demonstrates the utility of the 3D-PanoPACT in the label-free imaging and diagnosis of chronic cardiovascular diseases and liver lesions.

Cerebrovascular network anatomy and real-time functional dynamics in 3D-PanoPACT

Whole-brain imaging has become an indispensable approach for conducting in-depth research in neuroscience38, and 3D-PanoPACT proves to be competent in providing high-quality whole-brain anatomy and dynamic functional images non-invasively. To minimize acoustic loss and wavefront distortion, we removed a portion of the parietal bone with a diameter of ~8 mm and sutured the scalp without adding replacement. During the modeling process, to avoid overpressure, we carefully removed only the skull without causing any damage to the dura mater and employed low-dose mannitol to further alleviate brain swelling that may be induced by static overpressure (see “Methods”)39,40. The reliability of this modeling method for short-term brain imaging studies has been confirmed through structural imaging and pathological sections of the control groups41. A 10-Hz dual-wavelength illumination of 800 nm and 1064 nm was adopted to reveal more abundant cerebrovascular structures and obtain a larger penetration depth. To attain high spatiotemporal resolution within whole-brain FOV, we employed a method termed STINT, which involves the continuous and rapid rotation of the transducer array at small angular increments to synthesize an equivalent high-density detection array (refer to the final part of the Results section for details). With this, a comprehensive 3D representation of the whole brain anatomy was reconstructed, showcasing distinct main vessels and an abundance of branches (Movie S4). For visual clarity, the entire cerebrovascular network was categorized into cortical and deep-brain regions with multiple perspectives provided (Fig. 3a). Subsequently, a cut-open view of the deep-brain-region vessel network was depth-encoded along the Y-axis (Fig. 3b), and the image of each region was depth-encoded along the Z-axis (Fig. 3c, d), revealing intricate details with the names of vascular structures marked39,40. More strikingly, the 3D-PanoPACT provided clear observation of the Circle of Willis, which was representative basicranial vessels and acted as a crucial hub for ensuring blood supply and balancing blood flow throughout the entire brain. The arterial vascular composition (Fig. 3e) and the 3D location (Fig. 3f) of the Circle of Willis were distinctly observable at the base of the brain42,43, an application unattainable with other optical imaging modalities in living animals. The imaging depth of the whole brain in the real sense can be illustrated more intuitively by displaying different coronal sections from imaging results (Fig. 3g). These imaging results underscore the capability of 3D-PanoPACT in imaging the deep-brain anatomy, providing a high-quality structural reference for functional imaging.

Fig. 3: Label-free whole-brain anatomy in 3D-PanoPACT.

a Schematic diagrams illustrating regional segmentation of whole-brain imaging results from various observation perspectives. b The cut-open view of the deep brain with the cutting position and view orientation labeled in a (orange arrow). The image was depth-encoded along the Y-axis with the cutting position set to zero for better visualization. c The image of the cortical region depth-encoded along the Z-axis (yellow region in a). d The image of the deep-brain region depth-encoded along the Z-axis (blue region in a). e The 3D view of the Circle of Willis (displayed in red for emphasis) along with its constituent vascular structures. f The in situ location of Circle of Willis from the perspective of the base of the brain. g Six coronal sections from the image (labeled in a) showing the imaging depth in 3D-PanoPACT. Each cross-section was obtained by averaging a 1 mm thick 3D image in the thickness direction. OA olfactory artery, AACA azygos of the anterior cerebral artery, LHV longitudinal hippocampal vein, ISS inferior sagittal sinus, DTV dorsal thalamic vein, ACA anterior cerebral artery, IC internal carotid, RRV rostral rhinal vein, SV supraorbital vein, SSS superior sagittal sinus, SS sigmoid sinus, TS transverse sinus, CRV caudal rhinal vein, MCA middle cerebral artery, RE right eye, THA transverse hippocampal arteries, PCA posterior cerebral artery, ACoA anterior communicating artery, PCoA posterior communicating artery.

Furthermore, to demonstrate the superiority of imaging functional dynamics, the 3D-PanoPACT was applied in monitoring real-time hemodynamic response introduced by sodium nitroprusside (SNP) administration in the whole brain (Movie S2). As a commonly-used clinical antihypertensive hypotensor, SNP dilates the systemic vascular wall to alleviate the cardiac load44,45. In this experiment, a 10 Hz laser output (800 nm and 1064 nm) was adopted throughout the imaging process for unmixing operation (see “Methods”). To improve the signal-to-noise ratio (SNR) of deep-brain signals for high-precision unmixing results, the STINT method was adopted and the dual-wavelength image sequence at 0.5 Hz was produced, which fully met the requirements for capturing non-ultrafast changes in real time (refer to the final part of the Results section for details). Given that low-frequency information often accompanies functional responses such as cerebral blood oxygenation changes, we adopted the precisely measured EIR and applied deconvolution to the signals, enabling more accurate reconstruction of PA image amplitude (see “Methods”). A fixed dose of SNP was delivered from the caudal vein with a retention syringe over 180 s, and the whole brain was constantly monitored for 20 min. The estimation of sO2, as a pivotal physiological indicator, were demonstrated from the cortex to the Circle of Willis and displayed in parts acquired at 1 s, 290 s, and 1000 s after SNP injection (corresponding to Fig. 4a, b, and c respectively), to facilitate observation. To highlight the advantages of imaging deep-brain functional dynamics, two parameters, the sO2 and the total hemoglobin (HbT), were calculated in parallel for the Circle of Willis. In the whole-brain view, the average sO2 change for each major vessel is respectively plotted in Fig. 4d. Veins exhibited varying degrees of sO2 decrease, reaching a nadir around 300 s after SNP injection. The caudal rhinal vein (CRV) showed the maximum percentage decrease, approximately 40%, with sO2 values recovering to normal conditions in ~200 s. The reduction in cerebral sO2 resulted from an increase in oxygen extraction fraction caused by a simultaneous reduction in capillary perfusion46,47. In contrast, sO2 was maintained at almost 100% in the anterior cerebral artery (ACA), azygos of the anterior cerebral artery (AACA), and middle cerebral artery (MCA) due to their direct connection to the carotid artery. For the monitoring of the hemodynamics in the Circle of Willis, the ACA consistently maintained a high sO2 value, accompanied by a more pronounced change in HbT value (Fig. 4e). This could be attributed to the greater elasticity of the arterial wall. In contrast, the posterior communicating artery (PCoA) exhibited a relatively steep decline in sO2 values, while the HbT displayed a more gradual alteration (Fig. 4e). This observation suggested a potential role of the Circle of Willis in balancing and coordinating blood supply between the anterior and posterior parts of the whole brain. Therefore, the capability of imaging whole-brain functional dynamics in 3D-PanoPACT was underscored. Additionally, upon discontinuation of the anesthetic input, we monitored changes in sO2 in the major cerebral vessels during the recovery process to normal conditions (Fig. S7a–c), utilizing similar system setup and imaging methods. As illustrated, venous sO2 values gradually rebounded from hypoxia due to the progressively reawakening neural activities throughout the body, while the AACA consistently maintained a hyperoxic state, supporting fundamental cerebral activities (Fig. S7d). The inverse correlation between the average gradient of each sO2 curve and the physiological significance of the corresponding brain region may provide insights into the functional importance of the vessel’s location. In general, the 3D-PanoPACT demonstrates competence in noninvasively elucidating vessel-specific hemodynamic mechanisms throughout the entire brain, owing to its deep penetration and high spatial resolution.

Fig. 4: Real-time monitoring of whole-brain functional dynamics in 3D-PanoPACT.

a–e The 0.5-Hz-frame-rate imaging of cerebral hemodynamic response to sodium nitroprusside (SNP) administration. a–c The first column: the oxyhemoglobin saturation (sO2) map in the whole brain in response to SNP administration at 1 s, 290 s, and 1000 s (labeled in d) after SNP injection superimposed on the anatomy background; the second and third columns: the sO2 and total hemoglobin (HbT) maps of the Circle of Willis at the corresponding time point. d The sO2 change in whole-brain vascular (labeled in a) over 20 min following SNP administration (the administration time window is indicated by the shaded gray area). e The sO2 and HbT change in the Circle of Willis (labeled in a) over 20 min following SNP administration. TS transverse sinus, CRV caudal rhinal vein, ISS inferior sagittal sinus, LHV longitudinal hippocampal vein, DTV dorsal thalamic vein, SSS superior sagittal sinus, AACA azygos of the anterior cerebral artery, SS sigmoid sinus, ACA anterior cerebral artery, MCA middle cerebral artery, PCoA posterior communicating artery. f–k The 10-Hz-frame-rate imaging of whole-brain neurovascular coupling response under external electrical stimulation. f,** g** The map of sO2 change (%baseline) introduced by right and left forelimb stimulation, respectively. The top view and rear view of the responsive brain regions are displayed. MC motor cortex, SC sensory cortex. h The schematic curve of the electrical stimuli, the hemodynamic response function (HRF), and the convoluted hemodynamic response function (CHRF) used for extracting the responsive regions. i–k The averaged