Introduction

In recent years, benefiting from high-frame-rate imaging techniques and speckle tracking algorithms, ultrasound particle image velocimetry (uPIV, also known as echo-particle image velocimetry) has emerged as a representative tool for measuring blood flow fields based on contrast-enhanced ultrasound (CEUS) images1. High frame rate imaging (up to ten thousand frames per second) enables the reliable tracking of blood flow speckles2, and pulse-reversal imaging sequences enhances the signal-to-noise ratio of microbubble (MB) signals1,3,4,5,6,7. Compared to traditional Doppler imaging, uPIV measures flow velocity and direction with higher sensitivity and accuracy, unconstrained by imaging angles1,8. However, uPIV is unable to depict flow fields in microvessels due to the the diffraction limit.

Inspired by optical super-resolution imaging techniques9,10, super-resolution ultrasound (SRUS) imaging technology, including ultrasound localization microscopy (ULM) and non-localization-based SRUS strategies, employs microbubble to achieve non-invasive structural and functional imaging of in vivo microvasculature beyond the diffraction limit11,12,13,14. In contrast to CEUS and Doppler imaging techniques, SRUS enables measurements of flow velocity at the microvascular level13,14,15, providing sensitive microvascular biomarkers for early disease progression16,17,18. Therefore, SRUS holds high potential across a broad range of preclinical and clinical applications.

However, to achieve robust blood flow velocity measurements, precise localization and tracking of sparse microbubbles are required for ULM, inevitably leading to a fundamental trade-off between imaging speed and localization accuracy12,19. Meanwhile, long data acquisition time are necessary for the reconstruction of comprehensive microvasculature12,20, introducing challenge for ULM in clinical applications12,17,21,22,23.

High-concentration microbubble perfusion can enhance a larger proportion of the vascular system within the same duration12,21,24,25, but the spatial overlap of microbubble signals may increase false localizations and reduce accuracy. To enhance the speed of SRUS imaging, a series of nonlocalization-based SRUS strategies have been proposed, aiming to reconstruct microvasculature under high microbubble concentration.

For example, by introducing temporal auto-correlation to ultrasound data, such as super-resolution optical fluctuation imaging (SOFI), the spatial resolution of CEUS can be elevated by a factor of 1.4 within 0.3 second24,26. To further improve spatial resolution, sparsity-based ultrasound super-resolution hemodynamic imaging (SUSHI) further exploits the structural sparsity of microvessel, achieving a frame rate of 25 Hz25,27,28. Inspired by the optical super-resolution radial fluctuations (SRRF), ultrasound diffraction attenuation microscopy (UDAM) and entropy-based super-resolution (ERSR) leverage spatial gradient distribution characteristics and temporal fluctuations29,30. This combination results in notable improvement in spatial resolution within a 10-second acquisition time. Unfortunately, they do not provide measurements of blood flow dynamics, which limits the utility of these strategies in diagnosing microvascular-related diseases. Yu et al. utilized deconvolution to separate microbubbles from densely clustered microbubbles, achieving high temporal and spatial resolution imaging of the vasa vasorum of the rabbit femoral artery31. Considering that the point spread function (PSF) of the ultrasound imaging system is influenced by varying imaging, Qian et al. further employed pre-estimated PSFs at different imaging depths to reconstruct the peripapillary micrometer-scale vascular network32. However, these deconvolution methods are often sensitive to noise, and may suppress true signals under relatively low signal-to-noise ratio conditions, such as in high-speed plane-wave imaging32.

Meanwhile, deep learning-based approaches have made substantial progress in improving the localization accuracy of high-concentration microbubbles21,33,34,35,36. However, convolutional networks trained on simulated data have fundamental limitations in generalization performance across different organs23,35,36.

Beyond microbubble contrast agents, Zhang et al. proposed the Acoustic Wave Sparsely Activated Localization Microscopy (AWSALM) technique based on acoustically controllable nanodroplet contrast agents, achieving super-resolution reconstruction with sub-second temporal resolution37,38. However, the high mechanical index associated with this method may pose potential risks in clinical applications32.

Therefore, achieving rapid dynamic imaging of microvessels is a critical challenge.

In this study, we attempt to combine SRUS with uPIV, aiming to enhance the spatial resolution of uPIV while addressing the limitation of ULM in capturing dynamic blood flow information. To achieve this, we realize that by the construction of high-resolution ‘speckle fields’ reflecting motion characteristics of moving microbubbles, it may be possible to achieve dynamic evaluation of microflow speeds and directions. Consequently, we decompose the problem of microvascular blood flow imaging into two distinct steps: ‘reconstruction of high-resolution MB trajectory images’ and ‘speckle tracking/uPIV’, introducing a measurement tool for microvascular blood flow dynamics under high microbubble concentrations, denoted as enhanced ultrasound particle image velocimetry (E-uPIV).

To validate the feasibility of E-uPIV, we conducted experiments using simulated CEUS images, in vitro tube model with various flow rates, and in vivo data from mouse brain. These experiments demonstrate the feasibility and effectiveness of E-uPIV in enhancing spatial resolution and dynamically evaluating blood flow characteristics in high-density microbubble distributions.

Results

Simulation

The diagram of E-uPIV is illustrated in Fig. 1. Based on ultrasound simulation data with varying scatterer concentrations, we try to validate the advantage of E-uPIV in clearly imaging adjacent line structures. Simulated data cover a variety of MB concentrations (Supplementary Video 1).

Fig. 1: Workflow of enhanced ultrasound particle image velocimetry (E-uPIV) and conventional ultrasound particle image velocimetry (uPIV).
figure 1

A comparison of vector flow fields and velocity maps obtained using E-uPIV and uPIV is presented. In E-uPIV, contrast-enhanced ultrasound images are processed using a sliding time window of W frames across T consecutive frames. The in-phase and quadrature (IQ) data, generated after beamforming, undergo singular value decomposition (SVD) to produce contrast-enhanced images highlighting flowing microbubbles. MIP maximum intensity projection.

Full size image

Figure 2 displays images of maximum intensity projection (MIP) of simulated images, and SRUS results of second-order SOFI, ULM, UDAM, multiple signal classification algorithm (MUSICAL), mean-shift super resolution (MSSR), and E-uPIV under high MB concentration. ULM exhibits numerous false localizations due to MB overlap interference, therefore incapable of distinguishing neighboring streamlines. Supplementary Fig. 2 illustrates an increase in false localizations in ULM images with rising concentrations. Although MSSR, SOFI, MUSICAL and UDAM show improved discrimination between streamlines, their spatial resolution remains insufficient (Supplementary Fig. 3).

Fig. 2: Comparative analysis of simulation data results under high-density microbubbles.
figure 2

The intensity distributions at red lines of all enlarged rectangles are compared with the ground truth. CEUS MIP represents the maximum intensity projection (MIP) of contrast-enhanced ultrasound (CEUS) images. ULM ultrasound localization microscopy, UDAM ultrasound diffraction attenuation microscopy, SOFI super-resolution optical fluctuation imaging, MUSICAL multiple signal classification algorithm, MSSR mean-shift super resolution, E-uPIV enhanced ultrasound particle image velocimetry.

Full size image

In comparison, E-uPIV effectively avoids the introduction of false structures within streamlines, successfully separating closely situated structures at a distance of 64.2 μm (half wavelength: 98.7 μm). Notably, E-uPIV maintains consistent performance across both high and low MB densities (Supplementary Figs. 2 and 3), demonstrating robustness to varying MB concentrations.

For velocity measurement, as depicted in Supplementary Fig. 4, E-uPIV shows more accurate flow velocity measurements and better distinguishes local neighboring bifurcations compared to ULM under higher MB concentrations. The similarity coefficients between the E-uPIV/ULM flow velocity maps and the actual flow velocity map were 0.89 and 0.55 (p < 0.05), respectively. This highlights the feasibility of E-uPIV in reliably measuring flow velocity and direction, particularly under high MB concentration.

In vitro

We validated the accuracy of E-uPIV in measuring flow velocity using an in vitro flow channel model. The true velocity was calculated based on the volumetric flow rate provided by the infusion pump21. To assess the reliability of E-uPIV in constructing low and high-velocity flow fields, we established ten different flow rate gradients, ranging from 4 ml/h to 10 ml/h in 1 ml/h increments to simulate blood flow velocities ranging from 7.2 mm/s to 18 mm/s.

As depicted in Fig. 3, the MIP image of acquired CEUS images (a) indicates the complete perfusion of MBs throughout the entire channel. Under high MB concentrations, the localization and tracking results in ULM causes relatively sparse MB tracks, showing incomplete reconstruction within the tube (b, marked by white arrows). This is primarily attributed to the elimination of overlapping MBs and incorrect MB localizations.

Fig. 3: Comparative results of in vitro flow velocity measurements.
figure 3

a Maximum intensity projection (MIP) of acquired contrast-enhanced ultrasound images, depicting the flow channel phantom’s structure and position. Velocity measurements obtained under multiple flow rate conditions are shown for b ultrasound localization microscopy (ULM) and d enhanced ultrasound particle image velocimetry (E-uPIV). c The measured velocities (mean ± SD, N = 3 trials) are plotted as a function of channel flow rate for E-uPIV (green), ULM (yellow), and ground truth (blue), with error bars representing standard deviation.

Full size image

In contrast, E-uPIV images demonstrate better alignment of velocity measurements with the gold standard. We sampled multiple locations at a flow rate of 10 ml/h, obtaining diameter measurements of 411 ± 16 μm using E-uPIV and 504 ± 25 μm using ULM. Considering that the tube diameter is approximately 450 μm, these results indicate the reliability of the E-uPIV.

In Fig. 3c, velocity measurements obtained by E-uPIV and ULM are compared with the set channel flow rate. We selected a region with relatively stable flow distribution (indicated by the white box in Fig. 3d) to measure the average flow velocity. A clear linear relationship exists between E-uPIV results and channel flow rate (R2 > 0.99), while ULM fails to accurately reflect the true flow velocity. The overestimation of flow velocity in ULM may be attributed to the interference from wrong microbubble localization points and mismatches in tracking. This underscores the reliability of E-uPIV, especially in scenarios characterized by high MB concentrations.

In vivo

The comparison of E-uPIV and uPIV

E-uPIV has demonstrated satisfactory performance in the flow phantom experiment. Then, to validate the feasibility of E-uPIV in imaging in vivo microvasculature and monitoring microflow dynamics, we applied E-uPIV to the highly complex vascular structures in the mouse brain. The densely packed microvasculature poses a critical challenge for vascular flow imaging.

Figure 4 presents comparative analysis of the flow field obtained through conventional uPIV and E-uPIV. uPIV results were generated through performing particle image velocimetry algorithm on CEUS images. The flow dynamics measured by E-uPIV is shown in Supplementary Video 2.

Fig. 4: Local blood flow field reconstructed with E-uPIV and traditional ultrasound particle image velocimetry (uPIV) with 4-second acquisition time.
figure 4

The flow velocity of mouse brain vessels, measured by E-uPIV and uPIV, is shown in c and a, respectively. The red/yellow color in a and c indicate upward flow directions, while the blue/cyan color indicate downward flow directions. The blood flow fields measured by E-uPIV (c) and uPIV (a) in regions of interest, marked as 1 and 2 in a and c, are enlarged for detailed comparison in b and d. The arrows in the vector fields indicate the flow direction and the color reflects the flow velocity. The flow velocities measured by uPIV and E-uPIV along line L1 [marked as white lines in a and c] are plotted for comparison in e. The normalized intensity distribution of uPIV and E-uPIV along white lines L2 [marked in b and (d)] are compared in f.

Full size image

Due to the diffraction limit, uPIV is unable to clearly distinguish complex local arterial and venous blood flows (d). However, E-uPIV demonstrates substantial improvements in spatial resolution, revealing local complex microflow with varying flow directions (Fig. 4d). Further resolution estimates in Supplementary Table 2 shows that E-uPIV achieves a 5-fold enhancement in spatial resolution.

Additionally, the measured flow speeds of E-uPIV and uPIV are comparable in larger vessels (Fig. 4e), validating the feasibility of E-uPIV in tracking high-speed blood flows.

Overall, E-uPIV can measure local microblood flow fields with higher resolution while achieving reliable measurements of blood flow velocities in large-diameter vessels.

Spatial resolution improvement

Previous studies have demonstrated that accurate reconstruction of cerebral microvasculature using ULM typically requires very low microbubble (MB) concentrations and prolonged data acquisition times12,14,39. As illustrated in Fig. 5, under sparse MB conditions, ULM reconstructs only a limited fraction of vessels within the first 2 seconds, due to inefficient MB localization. Extending the acquisition time to 10 seconds significantly improves the completeness of the vascular network reconstruction.

Fig. 5: Vascular maps of the mouse brain reconstructed using E-uPIV and ultrasound localization microscopy (ULM) at varying acquisition times and microbubble (MB) concentrations.
figure 5

The vascular maps generated by E-uPIV are reconstructed using high-density MB perfusion. ULM results are shown for both high-density (HD) and low-density (LD) MBs. Acquisition times are indicated in each corresponding image. s=second. White scale bar= 1 mm.

Full size image

In contrast, E-uPIV enables rapid vascular reconstruction under high MB concentrations, with shorter acquisition times. While ULM enhances vascular details with longer acquisitions—ultimately achieving superior spatial resolution—E-uPIV and ULM yield highly similar reconstructions when processing the same low-MB-concentration dataset. Notably, E-uPIV provides clearer and more continuous vascular structures in short acquisition windows, whereas ULM excels in vascular continuity and capillary detail with prolonged imaging, surpassing E-uPIV in spatial resolution.

These findings indicate that under low-concentration microbubble conditions, ULM can obtain the highest-quality super-resolved microvascular images by extending the acquisition time. Meanwhile, in short acquisition durations, E-uPIV demonstrates notable advantages.

When ULM is applied under high microbubble concentration, overlapping MBs may lead to incorrect localizations, generating false structures and distortions that severely compromise the spatial resolution of traditional ULM21,23. E-uPIV circumvents this issue: even with a short 4-second acquisition, it preserves large-vessel integrity while enhancing microvascular branch continuity through spatiotemporal MB trajectory analysis.

For in vivo SRUS imaging, ULM results reconstructed from long-term sparse microbubble acquisition data serve as a reference due to the absence of a gold standard. A similarity comparison was conducted between E-uPIV images acquired within a short time frame and the reference ULM image (Fig. 6).

Fig. 6: Comparative analysis between E-uPIV images of high-concentration microbubbles acquired in a short time and ultrasound localization microscopy (ULM) images acquired under the low-density microbubbles over a longer duration.
figure 6

The white boxes indicate regions of interest, which are magnified in the middle panel. One-dimensional distributions at three different cross sections (marked with green annotations as I, II, and III) are compared. CEUS MIP represents the maximum intensity projection (MIP) of contrast-enhanced ultrasound (CEUS) images.

Full size image

Overall, the vascular structures in both ULM and E-uPIV images appear similar, while E-uPIV image displaying more continuous vascular structures. Detailed information extracted from locally sampled one-dimensional distributions shows that the vascular structures recovered by E-uPIV closely resemble those in the reference ULM image. The 2D correlation coefficient between the reference ULM (low concentration, 12 second) and E-uPIV (4 second) is 0.79 (p < 0.05), compared to 0.65 between ULM at low and high concentrations (12 seconds). This demonstrates E-uPIV’s advantage in achieving efficient microvascular imaging.

Feasibility in measuring microflow dynamics

In Fig. 7, velocity maps reconstructed by ULM and E-uPIV within a short acquisition time (4 seconds) are compared. To ensure the accuracy of velocity quantification, ULM was acquired under sparse MB distribution. E-uPIV was acquired under dense MBs.

Fig. 7: Ultrasound localization microscopy (ULM) and E-uPIV velocity maps reconstructed within different acquisition periods.
figure 7

Conventional color doppler image is shown. The red/yellow color in (a) and (c) indicate upward flow directions, while the blue/cyan color indicate downward flow directions. s=second.

Full size image

Within the short acquisition time, the flow direction and speed estimated by ULM appears highly chaotic due to limited MB localizations and unreliable velocity estimations19,21.

As the acquisition time extended, the velocity estimates for individual vessels gradually converge, demonstrating more homogeneous velocity profiles within the same blood vessels21,23.

In comparison, it can be observed that the dynamic range of blood flow velocity distribution reconstructed by E-uPIV does not exhibit evident changes with extended acquisition time, and the blood flow distribution within a single vessel remains stable. Whether in large vascular regions or microvascular areas, E-uPIV proves capable of reconstructing high-quality velocity maps with less than 4 seconds of acquisition data.

To assess E-uPIV’s capability in monitoring blood flow dynamics, as illustrated in Fig. 8, we conducted blood flow velocity measurements at several randomly sampled cross sections. In ROI 1, we observed a gradual decrease in peak velocity along the flow direction (Fig. 8b). In ROI 2, a reduction in peak velocities within the vein compared to the adjacent artery was evident (Fig. 8c). These observations are consistent with established physiological principles. These findings demonstrate E-uPIV had the advantage of providing a higher temporal resolution by enabling efficient high-resolution reconstruction within a short period of time.

Fig. 8: Local blood flow velocity fluctuations measured using E-uPIV.
figure 8

a Two regions of interest (ROIs) are selected. The color coding of flow direction is indicated with arrows. b Enlarged velocity maps and corresponding velocity profiles at selected cross sections within the ROIs. The sampled cross sections are marked as black lines. The light-shaded red region indicates the range of flow velocity fluctuations across the vessel’s cross-sectional area, while the solid line represents the mean velocity. c Speed measurements of a typical artery and vein at sampled cross sections (position indicated with white triangles) in ROI 2. The light-shaded red/blue region indicates the range of flow velocity fluctuations across the vessel’s cross-sectional area, while the solid line represents the mean velocity.

Full size image

In contrast, observing blood flow dynamics using ULM within short acquisition times proves challenging19,21. This is mainly due to the limited microbubble localization and the cumulative tracking required for ULM to reliably capture blood flow velocities (Fig. 7).

The results demonstrate that E-uPIV enables the reconstruction of microvascular dynamic changes with high feasibility and short imaging time.

Discussions

To enhanced the spatial resolution of uPIV and facilitate visualization of microvasculature, we introduced the idea of SRUS imaging into traditional uPIV and proposed a strategy termed E-uPIV. This strategy includes two steps: ‘reconstruction of high-resolution MB trajectory images’ and ‘speckle tracking/uPIV’. Initially, noise suppression and formation of MB trajectories are achieved through feature analysis of spatiotemporal MB data. Subsequently, we employ a spatial refinement scheme based on the mean shift theory to further shrink the initially enhanced MB trajectories, resulting in the high-resolution MB trajectory images. Through sliding window operations, the resolved MB trajectories exhibit dynamic behavior along the direction of blood flow, providing crucial features for subsequent tracking process. It is evident that the reconstruction of high-resolution MB trajectories with high signal-to-noise ratio is essential for achieving dynamic microflow field reconstruction in E-uPIV.

It is noteworthy that traditional ULM exhibits high sensitivity to the spatial distribution of MBs, generating chaotic flow estimates in regions characterized by MB overlap, and incomplete vascular reconstruction in areas with inadequate MB perfusion. Therefore, prolonged accumulation of sparse localizations is crucial to attain accurate estimations of vascular structure and flow velocity. In contrast, our findings show that E-uPIV is more robust under different MB concentrations. With the infusion of high-concentration MBs and the preservation of MB trajectories, E-uPIV can yield high-quality vascular structure reconstruction and precise flow velocity quantification with only 4 seconds data acquisition. This capacity offers technical support for the rapid and accurate analysis of vascular networks. In regions with multiple flow directions or high MB concentrations (e.g. large-diameter vessels), conventional ULM often suffers from evident structural artifacts and incomplete reconstruction40,41. Similar effects can be observed in the results of this study (Figs. 3 and 6).

We believe that the proposed E-uPIV emphasize the importance of preserving MB trajectory information for rapid ultrasound microvascular imaging. The previously-proposed nonlocalization-based SRUS studies have revealed that the exploration of the spatiotemporal characteristics of MBs can effectively preserves MB trajectories in the presence of high MB concentrations, while simultaneously reducing the size of PSFs29,30. These strategies prove to be instrumental in enhancing the efficiency of microvascular structure reconstruction. By applying traditional uPIV to high-resolution MB trajectory images, E-uPIV bypasses MB localization process, expediting the measurements of microflow dynamics through dense MB perfusion. In current non-parallel CPU configurations, E-uPIV requires extended computation time, often ranging from 40 minutes to an hour. When implemented with GPU acceleration, E-uPIV can achieve higher computational efficiency, and real-time blood flow imaging can be achieved

Given the challenges in clinical applications, such as short acquisition time, motion interference, and high MB concentration, ULM faces critical challenges in clinical implementation. Although many nonlocalization-based methods proposed by previous researchers have notably expedited the speed of SRUS imaging, they are incapable of measuring blood flow velocity21,24,25,29,30,41, which is a disadvantage compared to ULM. The proposed E-uPIV method, which could enable real-time tracking of microvascular hemodynamics, offering a promising strategy. Furthermore, E-uPIV also presents a practical solution for functional ultrasound imaging, particularly those involving the real-time monitoring of microvascular blood flow dynamics. While technologies like dynamic ULM are currently employed for tracking dynamic changes in cerebral blood flow42, repetitive task-induced stimuli are required to achieve reliable measurements of microvascular blood flow dynamics, thus unsuitable for dynamic velocity monitoring.

For PIV algorithms, high-intensity regions may lead to non-uniform flow estimation results43. Therefore, it is necessary to apply intensity upper limit filtering to high-resolution MB trajectory images, replacing all pixels exceeding the specified grayscale threshold with the threshold value. The size of the maximum inquiry window is determined through the careful sampling of MB displacements within high-velocity vessels.

Despite the advantages of E-uPIV, there are also limitations. The theoretical spatial resolution limit of E-uPIV is not as high as that of ULM under ideal conditions. In fact, for non-localization-based SRUS strategies, increasing the imaging speed requires sacrificing spatial resolution and vascular reconstruction integrity. Nonetheless, E-uPIV still offers evident improvement in spatial resolution compared to uPIV. Meanwhile, there is currently no well-established method that serves as the gold standard for the measurement of in vivo microvascular flow velocity. Therefore, we utilize the flow maps reconstructed with ULM, as a reference gold standard. To enhance the accuracy and robustness of the ULM velocity estimation, we employed sparser MB concentration and longer acquisition times. However, when comparing ULM results with E-uPIV, it is difficult to fully match them. On the one hand, the spatial resolution of E-uPIV is lower than ULM. On the other hand, E-uPIV encounters inherent physical limitations in CEUS imaging, whereby certain vessels with very low flow velocities may not encompass MBs within the acquisition time25. Despite the use of high-concentration MBs to maximize blood perfusion within a given time, comprehensive imaging of the microvascular system remains unattainable12,19,20,21.

Regarding imaging artifacts, E-uPIV, like conventional medical imaging methods, also have the probability of introducing artifacts. For example, there is a difference between the vessel geometry by the proposed approach and the ULM in the flow phantom results (Fig. 3). To minimize the generation of artifacts, incorporating new constraint factors to filter reliable microbubble trajectories, similar to the continuous microbubble tracking strategy used in ULM, can help exclude unreliable velocity estimates. Moreover, it is essential to balance parameter trade-offs and adjust them based on different application scenarios to minimize artifact interference. Meanwhile, it is essential to adjust the size of the time sliding window (w and T) and \({\sigma }_{0}\), selecting the optimal parameter combination based on visual assessments of the recovered image quality.

Following previous research44, this study selected a relatively simple hydrodynamic simulation model to validate the feasibility of the proposed method, which is sufficient to demonstrate E-uPIV’s ability to enhance spatial resolution and accurately distinguish neighboring structures. However, more advanced simulation studies are needed to consider additional factors, such as the dynamic velocity distribution across the vessel cross-section and the pulsatility of blood flow.

When imaging other moving organs, the performance of SVD may be compromised. Firstly, SVD operates under the assumption that tissue, blood flow, and noise signals are uncorrelated and lie in orthogonal subspaces. However, this assumption may break down in the presence of tissue motion or when imaging slow microvascular blood flow45. Secondly, the performance of SVD is sensitive to the selection of singular value thresholds8,46. These factors result in imperfect separation of blood flow and noise8, thereby reducing the accuracy of E-uPIV and ULM. To address this issue, multiple researchers has proposed a series of new strategies to suppress clutter and improve the signal-to-noise ratio47,48,49. These advancements contribute to enhancing the reliability of SRUS imaging.

In recent years, several studies have explored the value of high-frame-rate ultrasound imaging in revealing vascular structure and functions, introducing a series of high signal-to-noise ratio and high-resolution Doppler blood flow imaging techniques, which enable contrast-free observation of hemodynamics in small vessels50,51. Although these technologies struggle to surpass the diffraction limit, their unique advantages of contrast-free and rapid imaging make them more feasible for effective clinical application52. We believe that E-uPIV holds promise in achieving contrast-free microflow imaging. Furthermore, E-uPIV is expected to offer new perspectives for research in fluid dynamics or aerodynamics.

Methods

E-uPIV

E-uPIV offers a methodology for assessing micro-blood flow dynamics with exceptional resolution. This method bypasses MB localization, characterized by the reconstruction of high-resolution microbubble trajectory images, and the estimation of blood flow field by measuring speckle displacement fields between adjacent high-resolution microbubble trajectory images (Fig. 1).

Step 1. Reconstruction of high-resolution microbubble trajectory images

The reconstruction of high-resolution microbubble trajectory images is a crucial step. This study achieves this through the enhancement of spatio-temporal MB trajectories and the refinement of spatial resolution.

Enhancement of spatio-temporal MB trajectories

To enhance MB trajectories while reducing interference, E-uPIV employs sub-pixel confidence assessment based on statistical analysis of microbubble spatiotemporal data (CEUS images, Fig.1). Our approach involves the assumption that each microbubble PSF corresponds to a confidence region, with the highest confidence located at the center. Pixels outside of this region are considered to introduce additional noise. Inspired by multiple signal classification algorithm (MUSICAL)53, the confidence assessment is achieved through singular value decomposition (SVD) and indicator function computation53.

We consider a stack of CEUS images, denoted as \({I}_{{{{\rm{c}}}}}\,(c=1,\,2,\ldots ,T)\), where \(T\) represents the total number of frames. To reliably preserve the motion trajectories of microbubbles, we consider only small spatial window around a pixel and slide this window across all the pixels, instead of considering the whole image stack at once. The size of the spatial window is the approximate size of MB PSF, which is modeled as a symmetric Gaussian function with a standard deviation of \(\sigma\). To set the value of \(\sigma\), Gaussian fitting is applied to representative MBs randomly selected from the image stack. Then, the size of spatial window is set as (\(\omega {\times} \omega\)), where \(\omega =2\sigma\) and \(\omega \times \omega =N\). The center pixel of the sliding spatial window (\(\omega \times \omega\)) is \({{{{\boldsymbol{r}}}}}_{p}\). Meanwhile, temporal sliding window (\(W\) frames, where \(W\, < \,T\)) was employed for confidence analysis as well.

By conducting confidence assessment on these local spatiotemporal stacks (\(\omega \times \omega \times W\)), we generate local confidence patches across the entire image, which are then integrated to construct a final confidence map, exhibiting enhanced MB trajectories within \(W\) frames. Consequently, we obtain \(T-W\) confidence maps \({C}_{n}\,(n=1,2,\ldots ,\,T-W)\).

The implementation of this procedure is detailed in the following steps.

Step 1.1. The generation of confidence patches

We denote the centers of the regularly arranged pixels in the image plane as \({{{{\boldsymbol{r}}}}}_{i};i=1\) to \(N\), where \(N\) is the number of pixels and equal to the square of \(\omega\). It is assumed that there are fluctuations in the pixel values across frames within the temporal window. Microbubbles flowing in different streamlines move independently while microbubbles flowing in the same streamline move in a correlated fashion. Then, the local spatiotemporal stack (\(\omega \times \omega \times W\)) centered at the pixel \({{{{\boldsymbol{r}}}}}_{p}\) is denoted as \({{{{\boldsymbol{I}}}}}_{p}\). In this context, \({{{{\boldsymbol{I}}}}}_{p}\) is reorganized into a \(N\times W\) two-dimensional matrix, in which each row corresponds to a pixel and each column corresponds to a frame. The complete local spatiotemporal stack which contains images can be written as a matrix:

$${{{{\boldsymbol{I}}}}}_{p}=[{{{{{\boldsymbol{I}}}}}_{p}}^{1}{{{{{\boldsymbol{I}}}}}_{p}}^{2}\ldots {{{{{\boldsymbol{I}}}}}_{p}}^{W}]$$
(1)

where the intensity measurement in the k th frame at all pixels can be written as:

$${{{{{\boldsymbol{I}}}}}_{p}}^{k}={[{{I}_{p}}^{k}\left({{{{\boldsymbol{r}}}}}_{1}\right){{I}_{p}}^{k}\left({{{{\boldsymbol{r}}}}}_{2}\right)\ldots {{I}_{p}}^{k}\left({{{{\boldsymbol{r}}}}}_{N}\right)]}^{T}$$
(2)

Due to the use of sliding window, the MB trajectories is most reliable in the center pixel. Then, we applied soft window function to transform \({{{{\boldsymbol{I}}}}}_{p}\), scaling the weights of the pixels according to their proximity to the center pixel. The spatial distribution of soft window function is given as:

$${{{\boldsymbol{h}}}}=[h\left({{{{\boldsymbol{r}}}}}_{1}\right)h\left({{{{\boldsymbol{r}}}}}_{2}\right)\ldots h\left({{{{\boldsymbol{r}}}}}_{N}\right)]$$
(3)

where \(h\left({{{\boldsymbol{r}}}}\right)\) represents the weights of symmetric Gaussian function centered at the pixel \({{{{\boldsymbol{r}}}}}_{p}\):

$$h\left({{{\boldsymbol{r}}}}\right)=exp \frac{{||r}-{r}_{p}{{||}}^{2}}{2{\sigma }^{2}}$$
(4)

We denote the diagonal soft window transformation matrix as \({{{\boldsymbol{H}}}}={diag}({{{\boldsymbol{h}}}})\), and the soft-windowed local spatiotemporal stack is obtained as:

$${{{{\boldsymbol{J}}}}}_{{{{\boldsymbol{p}}}}}={{{\boldsymbol{H}}}}{{{{\boldsymbol{I}}}}}_{p}$$
(5)

Then, SVD is used to compute the eigen-images and corresponding eigenvalues of the soft-windowed local spatiotemporal stack \({{{{\boldsymbol{J}}}}}_{{{{\boldsymbol{p}}}}}\), where each feature image \({{{{\boldsymbol{\mu }}}}}_{i}\) correspond to different structural details of the MB signals in the local spatiotemporal stack \({{{{\boldsymbol{I}}}}}_{p}\), and the corresponding eigenvalue \({\sigma }_{i}\) represent the energy. It can be inferred that the eigen-images with high singular values represent the local MB trajectories, and statistically less likely patterns (include noises), are characterized by small singular values. Therefore, we can choose first few eigen-images with the large eigenvalues as the representative of the structure. Meanwhile, we assume that the presence of noise may corrupt the details in the eigen-images with smaller eigenvalues.

Quantitively, similar to MUSICAL and MUSIC algorithms, we set a threshold \({\sigma }_{0}\) such that the corresponding eigen-images can be designated as null subspace N: {\({{{{\boldsymbol{\mu }}}}}_{{\sigma }_{i} < {\sigma }_{0}}\)} and signal subspace R: {\({{{{\boldsymbol{\mu }}}}}_{{\sigma }_{i}\ge {\sigma }_{0}}\)}. The choice of the value of \({\sigma }_{0}\) is illustrated in the Supplementary Note 1.

Based on this information, we can evaluate the confidence level of MB trajectories presented in certain point based on separated spaces N and R. Specifically, we computed indicator function value \(f\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)\) at every test point \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\). Initially, we assumed the PSF image at \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\) as \({G}^{{\prime} }({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} })\), represented as:

$${G}^{{\prime} }\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)={[G\left({r}_{1},{r}_{{test}}^{{\prime} }\right)G\left({r}_{2},{r}_{{test}}^{{\prime} }\right)\ldots G\left({r}_{N},{r}_{{test}}^{{\prime} }\right)]}^{T}$$
(6)

Here, \(G\left(r,{r}_{{test}}^{{\prime} }\right)\) denotes the value of the pixel \({{{\boldsymbol{r}}}}\) if there is a single microbubble at point \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\), derived based on the previously established PSF model. Here, \({{{{\boldsymbol{r}}}}}_{i};i=1\) to \(N\) represents the centers of the discrete pixels in the image space. The test point \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\) are set at subpixel level, for MB locations are not confined to discrete grid points. Then, we calculate the projections of \({G}^{{\prime} }({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} })\) onto the signal subspace R and null subspace N, denoted as \({d}_{V}({r}_{{test}}^{{\prime} })\) and \({d}_{N}({r}_{{test}}^{{\prime} })\) respectively:

$${d}_{V}\left({{{\boldsymbol{r}}}}_{{test}}^{{\prime} }\right)=\sqrt{{\sum }_{{\sigma }_{i}{\ge \sigma }_{0}}{||G}^{\prime} ({{r}^{{\prime} }}_{{test}})\cdot {{{\boldsymbol{\mu }}}}_{i}{{||}}^{2}}$$
(7)
$${d}_{N}\left({{{\boldsymbol{r}}}}_{{test}}^{{\prime} }\right)=\sqrt{{\sum }_{{\sigma }_{i} < {\sigma }_{0}}{||G}^{\prime} ({{r}^{{\prime} }}_{{test}})\cdot {{{\boldsymbol{\mu }}}}_{i}{{||}}^{2}}$$
(8)

These projections indicate whether the PSF at the test point \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\) is associated with the patterns represented by the range. Subsequently, an indicator function \(f\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)\) is formulated as:

$$f\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)={\left(\frac{{d}_{V}\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)}{{d}_{N}\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)}\right)}^{2}$$
(9)

The indicator function’s value estimates the confidence level at each sub-pixel, which is high at pixels located within MB trajectories and significantly low where microbubbles are absent. Theoretically, MB locations \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\) can be any points in the sliding window. However, for further computation and stitching convenience, we decimate every pixel into a sub-pixel grid with even number of grid points along each direction and take the centers of the sub-pixels as the test point \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\).

Repeatedly applying the same procedure at every sub-pixel yields a confidence patch \({{{{\boldsymbol{M}}}}}_{p}\) centered at pixel \({{{{\boldsymbol{r}}}}}_{p}\), revealing enhanced MB trajectories hidden within local spatiotemporal stack \({{{{\boldsymbol{I}}}}}_{p}\) with a size of \(\omega \times \omega \times W\).

Step 1.2. Stitching all confidence patches

Through the use of \({d}_{V}\left({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\right)\), confidence patch \({{{{\boldsymbol{M}}}}}_{p}\) of each sliding window is automatically scaled with reference to the local spatiotemporal stack \({{{{\boldsymbol{I}}}}}_{p}\) and the confidence patches positioned in the relevant pixel can simply be added for stitching all patches.

We assume \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\) is located in the q th grid pixel: \({{{{\boldsymbol{r}}}}}_{q}\). To calculate the final confidence level at \({{{{\boldsymbol{r}}}}}_{{test}}^{{\prime} }\), we gathered all values of indicator functions at positions where pixel \({{{{\boldsymbol{r}}}}}_{q}\) is located within corresponding sliding windows. Then, the final confidence map \({C}_{n}\) representing microbubble information within \(W\) frames, is obtained.

As temporal-window analysis is performed using a defined sliding window (size: \(W\) frames, \(T\) frames in total, \(W < T\), step: 1 frame), the above procedure results in a sequence of two-dimensional confidence maps \({C}_{n}(n={1,2},\ldots ,{T}-W)\), displaying dynamic MB trajectories. Detailed algorithm was provided in53. As depicted in Fig.1, the MIP of MB trajectories demonstrate more refined MB trajectories.

Refinement of spatial resolution

To further enhance the spatial resolution of confidence maps \({C}_{n}\,(n={1,2},\ldots ,{T}-W)\), we employed a spatial refinement approach based on the mean shift theory to further shrink the initially-enhanced MB trajectories. Specifically, inspired by mean-shift super resolution (MSSR)54, we calculate the mean shift (MS) vector for each pixel locally by a kernel window that slides throughout the entire image using a spatial-range neighborhood.

Assuming one pixel as \({p}_{0}=[{x}_{0},y\left({x}_{0}\right)]\),where \({x}_{0}\) represents the coordinate and \(y\left({x}_{0}\right)\) represents the grayscale of \({C}_{n}\) at position \({x}_{0}\). We define the spatial-range bandwidth as \(h=({h}_{s},{h}_{r})\), in which \({h}_{s}\) denotes the spatial bandwidth and \({h}_{r}\) denotes an intensity range. Then the neighborhood around pixel \({p}_{0}\) is defined as \({B}_{h}({p}_{0})\), which excludes only the host pixel:

$${B}_{h}\left({p}_{0}\right)=\{p=\left[x,y\right]\!\!:0\, < \,{||}{x}_{0}-{x||}\, < \,{h}_{s},0\, < \,{||y}\left({x}_{0}\right)-{x||}\, < \,{h}_{r}\}$$
(10)

We define \({p}_{{{{\rm{n}}}}}\in {B}_{h}\left({p}_{0}\right)\). Then the distance between \({p}_{{{{\rm{n}}}}}\) and \({p}_{0}\) is measured using a common profile \(k\) and denoted as \({k}_{h}({p}_{0},{p}_{0})\):

$${k}_{h}\left({p}_{0},{p}_{{{{\rm{n}}}}}\right)={k}_{s}({t}_{s})\cdot {k}_{r}({t}_{r})=k\left(\left|\left|\frac{{x}_{0}-{x}_{n}}{{h}_{s}}\right|\right|^{2}\right) \cdot k\left(\left|\left|\frac{{y}({x}_{0})-y({x}_{n})}{{h}_{r}}\right|\right|^{2}\right)$$
(11)

\(k\) is defined as gaussian profiles. For a given neighborhood, \({p}_{0}\) are \(h\) are fixed. As the position moves away from the central pixel \({p}_{0}\), the distance \({t}_{s}\) increases, and the value of \({k}_{s}\) decreases. The same phenomenon happens with \({k}_{r}\).

The partial derivatives of \({k}_{h}\left({p}_{0},{p}_{{{{\rm{n}}}}}\right)\) is denoted as:

$${g}_{h}\left({p}_{0},{p}_{{{{\rm{n}}}}}\right)={{k}_{s}}^{{\prime} }\cdot {{k}_{r}}^{{\prime} }=g\left(\left|\left|\frac{{x}_{0}-{x}_{n}}{{h}_{s}}\right|\right|^{2}\right)\cdot g\left(\left|\left|\frac{{y}({x}_{0})-y({x}_{n})}{{h}_{r}}\right|\right|^{2}\right)$$
(12)

The MS vector at \({p}_{0}\) can be derived through calculating partial derivatives of \({k}_{h}\) respect to \({t}_{s}\) and \({t}_{r}\) over a spatial-range space, formulated as:

$${{MS}}_{h}\left({p}_{0}\right)=\frac{{\sum }_{n}{g}_{h}({p}_{0},{p}_{{{{\rm{n}}}}}){y}({x}_{{{{\rm{n}}}}})}{{\sum }_{n}{g}_{h}({p}_{0},{p}_{{{{\rm{n}}}}})}-{y}({x}_{0})$$
(13)
$${{SM}}_{h}\left({p}_{0}\right)=\frac{{\sum }_{n}{g}_{h}({p}_{0},{p}_{{{{\rm{n}}}}}){y}({x}_{{{{\rm{n}}}}})}{{\sum }_{n}{g}_{h}({p}_{0},{p}_{{{{\rm{n}}}}})}$$
(14)

The term \({{SM}}_{h}\left({p}_{0}\right)\) represents the weighted local mean in the spatial-range space at pixel \({p}_{0}\). Therefore, the MS vector refers to the difference between the weighted mean or sample mean (SM) of all elements of the analyzed region \({B}_{h}\left({p}_{0}\right)\). The magnitude of MS depends on the difference between the central pixel of the neighborhood and the surrounding pixels, and always points towards the direction of the intensity gradient. It reflects that the MS theory was designed to move away from minimum local density points and converges towards higher density regions54.

For digital image, we can conclude that points of minimum local density are assigned to large values of MS. For an isolated microbubble PSF, the local maximum density point corresponds to the peak of the Gaussian distribution. Therefore, the minimum MS value calculated from the Gaussian distribution matches the point of maximum intensity of the initial distribution, while positive values of MS indicate the absence of MB PSF54. Previous studies demonstrate that, if the negative MS values are inverted and those less than 0 are discarded, it helps reduce the size of the PSF, resulting in an image with higher spatial resolution54.

Therefore, by calculating the value of MS vector at every pixel in the confidence maps \({C}_{n}(n={1,2},\ldots ,{T}-W)\), we would get high-resolution MB trajectory images \({T}_{n}(n={1,2},\ldots ,{T}-W)\), denoted as:

$${T}_{n}=-{S}_{n}$$
(15)

As illustrated in Fig. 1, compared to confidence maps \({C}_{n}\), artifacts in \({T}_{n}\) are suppressed and spatial resolution is improved, generating images with refined microbubble trajectories.

Step 2. Particle image velocimetry algorithm

After constructing the high-resolution MB trajectory images \({T}_{n}(n={1,2},\ldots ,{T}-W)\), we further measured blood flow field based on typical PIV algorithms55,56, which was implemented using MATLAB (R2015a, The MathWorks Inc.). With typical PIV algorithms integrated, the toolbox provided in55,56 has been extensively applied in scenarios such as cell motion particle flows57,58,59, as well as in ultrasound images60.

In this study, we achieve flow field reconstruction through image cross-correlation analysis. Cross-correlation is a pattern matching technique that calculates the discrete cross-correlation function of sub-images (interrogation regions) to determine the most likely particle displacements within the local area61. To enhance processing speed, this study employs traditional discrete Fourier transform (DFT) to perform cross-correlation in the frequency ___domain55.

Simultaneously, to obtain high spatial resolution and precise displacement vectors, we adopted a multi-pass window deformation approach58. In each iteration, the interrogation region is offset and refined based on the previous iteration’s displacement estimate55,62.Previous studies have found that this multi-pass window deformation approach is particularly beneficial to minimize information loss caused by particle displacement, aid convergence to true auto-correlations, and yield a high spatial resolution in the final vector map55,61,62.

After each iteration, a bicubic interpolation is applied to deform the previous inquiry region55. Between the passes, the velocity information is smoothed and validated and missing information is interpolated55. To mitigate noise interference, a 3\(\times\)3 Gaussian kernel is used to smooth the correlation matrix. After obtaining the correlation matrix, the displacement of two interrogation regions can be directly obtained by the peak positions55,63. In this study, by fitting standard Gaussian function to the integer intensity distribution, the ___location can be refined with sub-pixel precision55,64.

ULM

Ultrasound data is saved as radiofrequency (RF) data and subsequently processed offline using MATLAB. After beamforming, In-phase and quadrature (IQ) data were generated.

To extract microbubble signals from each ultrasound IQ dataset, a clutter filter based on spatiotemporal SVD is applied41,65,66. The adaptive determination of the low-order singular value threshold, representing tissue signals, is achieved67, generating CEUS images.

To model the PSF of the system, we manually identified isolated MB signals and fitted them to a multivariate Gaussian function as empirical PSF. To localize MB centers, we evaluated the 2D normalized cross-correlation between the empirical PSF and each CEUS frame. A threshold is applied to the cross-correlation map to remove pixels with correlation coefficient values below 0.6, thereby suppressing noise and retaining high-confidence single MB images. Subsequently, the intensity-weighted centroid of each individual MB image is calculated to locate the MB center. Then, we employed the Hungarian algorithm for frame-by-frame matching of moving MBs to calculate flow velocity based on MB displacement and to generate MB trajectories. Accumulating microbubble trajectories for vascular structure reconstruction and reconstructing blood flow velocity maps through time-averaging of MB flow velocities.

Statistics and Reproducibility

The 2D correlation coefficient between the reference ULM (low concentration, 12-second acquisition) and E-uPIV was quantified using the 2D Pearson correlation coefficient (R). The Pearson R values were calculated directly from pixel-wise intensity comparisons of the two images, with values ranging from −1 (perfect inverse correlation) to +1 (perfect positive correlation). A value of R = 0 indicates no linear relationship. The analysis was performed in MATLAB, using the built-in function (corr2 in MATLAB).

Materials

Simulation study

To compare the accuracy of different SRUS strategies at various concentrations, we simulated moving scatterers (MBs) within a Y-shaped flow ___domain and generated simulated images using the Verasonics simulator (Verasonics Inc., Kirkland, Washington). Detailed parameters are provided in Supplementary Table 2. The Y-shaped flow simulates the cross-section of bifurcating vessels. MB velocities in the vessels on the right and left sides were set at 8.70 cm/s.

In our simulated data, scatterers move downward at a constant speed from above a vertical line. At bifurcations, they flow into either the left or right branches with equal probability and at the same velocity. We control the concentration of scatterers within the imaging area by adjusting the absolute entry time (delay).

In the high concentration scenario defined in our study, 509 microbubbles entered within 1000 frames, with a spacing range of 0 to 2λ (λ: wavelength). For medium concentration, 276 microbubbles entered with a spacing range of 0 to 5λ, and for low concentration, 159 microbubbles entered with a spacing range of 0 to 10λ. Supplementary Video 1 illustrates these scenarios, highlighting the distinct differences in scatterer density. Notably, in high concentration scenarios, microbubbles overlap at bifurcation zones, with overlap of adjacent microbubbles within individual streamlines.

In vitro study

In order to validate the accuracy of E-uPIV in flow velocity measurements, we fabricated a phantom model. We drilled a pair of orifices on the side walls of a transparent plastic container, and a rubber tube with a diameter of 500 μm was passed through these orifices. Agarose gel was poured into the container until it was approximately 5 mm above the tube’s surface. The phantom was then placed in a 4 °C refrigerator until the agarose gel had completely solidified, at which point the tube diameter had contracted to approximately 450 μm. To create different flow velocities, we connected the inlet of the flow channel to a microfluidic pump to provide controlled and constant volume of injected liquid within a fixed time period. After dissolving the MBs in 5 mL of 0.9% sterile sodium chloride, it was further diluted 1000-fold and perfused into the channel through a syringe fixed on a microfluidic pump device.

In this experiment, we established injection rates ranging from 4 ml/h to 10 ml/h in 1 ml/h increments to simulate blood flow velocities ranging from 7.2 mm/s to 18 mm/s. The mean velocity was calculated using the flow rate of the connected syringe pump and the known inner diameter of the tube. Ultrasound data were acquired using the Verasonics Vantage 128 platform and the L35-16vX linear array transducer, with a center frequency of 28 MHz and a pulse repetition frequency (PRF) of 25000. Specifically, we placed the L35-16vX transducer at the top of the model with a transducer holder to enable imaging of the longitudinal cross-section of the pipeline. For each flow rate, data acquisition commenced after approximately 10 seconds to guarantee a stable flow, and a total of 2000 frames with a 500 Hz frame rate were collected.

Ultrasound data is saved as RF data and subsequently processed offline using MATLAB. After beamforming, IQ data were generated, and singular value decomposition was utilized to enhance microbubble signals. The spatio-temporal MB trajectories were reconstructed with a subpixel size of 0.0086 mm and a temporal window size of 50 frames. Then, two-dimensional uPIV algorithm is employed, with four interrogation window size, including 60 pixels (0.52 mm), 30 pixels, 10 pixels and 5 pixels, with a 50% overlap between iterations.

In vivo study

All experimental procedures involving laboratory animals were reviewed and approved by the local animal care committee of Peking University (AAIS-ZhangJ-8). In this study, female KM mice (8 weeks old; Charles River Laboratories) were used for cerebral blood flow imaging. Prior to imaging, mice were anesthetized using a 3% isoflurane gas induction chamber. Subsequently, the mice were positioned within a stereotaxic frame equipped with a nose cone, allowing for head fixation. Anesthesia was maintained at 2% isoflurane. Once the mice achieved a stable condition, the scalp was removed, and a cranial window was created. Ultrasonic imaging was performed using the Verasonics Vantage 256 system. The L35-16vX transducer was fixed to maintain the stability of data acquisition (center frequency: 28 MHz). Ultrasound gel was applied directly to the surface of the mouse brain. 59 mg of Sonovue was diluted in 5 mL of saline and administrated via the tail vein at a dosage of 0.1 mL/kg and an injection volume of 200 μL. Imaging is performed using a 7-angle compound plane wave sequence at a 1-degree increment, resulting in a compounded frame rate of 500 Hz. 60,000 frames of data are acquired within 120 seconds acquisition time. After microbubble injection, we utilized the data acquired initially as a typical sequence of high-concentration CEUS images (0-20 s). As the acquisition time progressed and microbubbles became sparse, we selected CEUS data from the interval when microbubbles could be distinctly differentiated (after 60 s) as a typical sequence of low-concentration CEUS images (Supplementary Fig. 1).

Ultrasound data is saved as RF data and subsequently processed offline using MATLAB. After beamforming, IQ data were generated and SVD filter was applied to enhance microbubble signals. Then, the spatio-temporal MB trajectories were reconstructed with a subpixel size of 0.002 mm and a temporal window size of 50 frames.

To track the local displacement of moving MBs in high-resolution MB trajectory images, two-dimensional PIV algorithm is employed. This algorithm utilizes cross-correlation analysis of iterative deformation window with four iterations, progressively reducing the interrogation window size from 50 pixels (0.1 mm) to 5 pixels (0.01 mm) with a 50% overlap between iterations. Window deformation and Gaussian sub-pixel estimator are employed to enhance the accuracy of displacement estimation.

Referring to previous study37,68,69, we measured the spatial resolution at the corresponding imaging frequency. Specifically, we sampled 10 independent MB PSFs and calculated the axial and lateral full width at half maximum (FWHM), resulting in values of 143 ± 6 μm and 173 ± 27 μm, respectively.

Reporting summary

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