In vivo microvascular flow quantification in the mouse brain using row–column ultrasound localization microscopy and directed graph analysis

Abstract

Brain perfusion relies on a complex vascular network of arteries, veins, and capillaries to meet its constant demand for oxygen and nutrients. Disruption of this microvascular system is a hallmark of many neurological disorders, including small vessel disease, stroke, and brain tumors. As such, high-resolution in vivo imaging of cerebral microvascular flow and structure remains critical to understanding these pathologies. Among them, ultrasound localization microscopy (ULM) allows noninvasive imaging of microvascular network down to small arterioles and venules at subwavelength resolution using injected microbubbles, but the approach remains mainly limited to 2D imaging with few volumetric implementations. In this study, we explore in vivo transcranial 3D ULM of the mouse brain using row–column arrays (RCA) and introduce an analysis framework to build a flow-directed vascular graph from the ULM microbubble tracking data, allowing to differentiate between subgraphs of artery-like and vein-like vascular segments. Using this framework, we are able to quantify flow and radius relationships for each subgraph in different anatomical regions. This high-sensitivity framework enables in vivo microvascular imaging and quantification in mice and provides a scalable platform for preclinical neurovascular studies in health and disease.

Introduction

Preclinical neuroimaging has become an indispensable tool in neuroscience research for imaging the brain’s vascular network, a highly intricate system responsible for regulating cerebral blood flow, oxygen and nutrient delivery, and metabolic waste clearance. This network is fundamental to maintaining brain health, yet many aspects of its structure and function remain poorly characterized. Emerging evidence suggests that alterations in microvascular structure and function may serve as early biomarkers for a range of neurological disorders. ¹ For instance, small-vessel diseases and neurodegenerative disorders—such as Alzheimer’s disease—are often preceded by microvascular alterations, ² and such alterations may contribute to disease progression before overt symptoms appear. In vascular stroke, phenomena such as vasospasms, aneurysms, and delayed ischemia disrupt blood flow across multiple spatial and temporal scales and are critical areas of study for developing and testing effective recovery and intervention strategies. ³ In brain tumors, the formation of aberrant and complex vascular networks is essential for supporting rapid tumor growth and survival. Imaging these networks is not only vital for understanding tumor biology but also for developing targeted therapeutic approaches to disrupt their blood supply and limit progression.

Several imaging modalities have been developed for neurovascular mapping, including magnetic resonance angiography (MRA)^4,5 and computed tomography angiography (CT) using contrast agents but also optical imaging techniques. ⁶ However, MRA and CT imaging techniques are limited in their ability to provide both high-resolution imaging and blood flow dynamics, while optical imaging cannot image deeper brain regions. While advanced ex vivo techniques, such as optical clearing methods (e.g. iDISCO, CLARITY) ⁷ and synchrotron-based phase-contrast tomography,^8–10 excel in providing unparalleled spatial resolution for microvascular imaging, they come with significant limitations. These techniques are highly effective for detailed structural visualization of cerebral vasculature, often achieving resolutions at the sub-micron to micron scale but lack the ability to capture dynamic flow information and are not suitable for longitudinal studies.

These limitations highlight the need for novel imaging modalities that can provide both high-resolution structural imaging and dynamic functional assessment of the brain’s microvasculature in vivo. Ultrasound localization microscopy (ULM) is an emerging imaging technique that can map microvascular flow at high spatial and temporal resolution. ULM relies on the imaging, localization, and tracking of FDA-approved gas microbubbles circulating in the blood flow ¹¹ using ultrafast ultrasound imaging. The technique has been applied to map cerebral microvascular flow in rodents and in humans using 2D imaging approaches allowing subwavelength resolution in depth and providing quantitative velocity estimates. ¹² Recently, functional ULM (fULM) has been introduced to map functional hyperemia in the rodent brain at high spatial resolution and in time. ¹³

While ULM has significant potential to increase our understanding of cerebral microvasculature, its generic two-dimensional (2D) implementation presents important field of view limitations but also introduces systematic biases in velocity measurements, as only two components of the velocity vector can be reliably captured. These constraints have motivated the development of three-dimensional (3D) ULM approaches using dense ultrasonic matrix arrays, composed of thousands of piezoelectric elements arranged in 2D grid for rodent brain imaging including stroke and tumor models.^14–18 Despite early proofs of concepts, matrix arrays face substantial technical challenges. Manufacturing remains costly and technically demanding, particularly at the high frequencies and small element pitch. ¹⁹ This makes scaling to larger fields of view challenging without ever increasing the number of elements, escalating costs, and system complexity. While electronic multiplexers can mitigate some of these requirements by dynamically switching element,^17,20 this approach reduces the effective volume rate due to the sequential addressing of array elements, compromising microbubble tracking efficiency ²¹ and overall image quality. ²²

To address these limitations, row–column array (RCA) probes have emerged as a promising alternative for 3D high volume rate ultrasound imaging. Conceptually, RCA probes consist of two overlapped orthogonal one-dimensional arrays which offers a simplified, yet efficient approach compared to matrix arrays. RCA enables larger field-of-view imaging with reduced complexity while maintaining improved sensitivity and can be manufactured with finer pitches even at high frequencies. They further enable a significant reduction of channels compared to matrix arrays (N + N rather than N², N being the number of elements along one side of the array), and remains compatible with standard 256-channel preclinical ultrasound scanners. The versatility of RCA probes in preclinical neuroimaging has been demonstrated across various modalities, including Doppler imaging, ²³ functional ultrasound imaging (fUS) ²⁴ and non-linear contrast imaging of gas vesicles. ²⁵ Their potential has been explored in proof-of-concept studies of 3D ULM^26,27 including recently 3D imaging of the mouse brain by Wu et al. ²⁸ and of the rat brain by Sun et al. ²⁹ While RC-ULM shows clear promise for 3D ULM imaging in rodent brains, further advancements in acquisition protocols, reconstruction methods, and analytical frameworks are still needed toward achieving a high-throughput, benchside microvascular flow quantification technique suitable for in vivo preclinical studies.

Our study explores transcranial 3D ULM in mouse models using a dedicated RCA probe through the combination of optimized acquisition sequences, advanced reconstruction algorithms, and automated flow-directed graph-based analysis of vessels. We first demonstrate the capability of this approach to reconstruct the 3D cerebral microvasculature through the skull with high sensitivity and spatial resolution, down to small arterioles and venules, including the local estimation of blood velocities, and registration to the Allen Brain Atlas. We further show that these data can enable the construction of a flow-directed vascular graph which locally encodes flow direction and can then be used for labeling artery-like and vein-like segments based on local branching and flow gradients. We further extend the analysis to the automated estimations of dynamic features (flowrate, microbubble velocity) and vessel structural characteristics (segment length, tortuosity) across different scales and anatomical regions of the brain.

Materials and Methods

Ethics

The institutional and regional committees approved the study for animal care (Comité d’éthique pour l’expérimentation animale/#59 – “Paris Centre et Sud,” Protocole/#2017-23) and adhered to the ARRIVE recommendations. All animals received humane care in accordance with the European Union Directive of 2010 (2010/63/EU).

Mice were naive, 7 weeks old at the time of the experiment, and randomly placed in the cage; none were excluded. The individual animal is the experimental unit in this investigation.

Animal preparation

Five male C57Bl/6 mice (Janvier Labs, Le Genest St Isle, France) were utilized in this study. The mice were housed under controlled conditions (22 °C ± 1 °C, 60% relative humidity, 12/12 h light/dark cycle, ad libitum food and drink) for a week prior to the start of the experiment. The mice were anesthetized with an initial administration of ketamine (Imalgene, 100 mg/kg) and Xylazine (Rompun) intramuscularly to render them unconscious, and anesthesia was maintained with a 1.5% isoflurane supply and the animal physiology was monitored. The animals were then placed in a stereotaxic frame. The skin and the skull were kept intact. The hair was shaved from the mouse’s head using a trimmer (ISIS GT 421 Aesculap, PhyMep, France) and depilatory cream (Vichy, France). The eyes of the mouse were protected using an ointment (Ocry-gel, TVM, UK). To ensure the stability and reproducibility of the animal state during imaging, body temperature, heart rate, and breathing rate were monitored and maintained within specified ranges throughout the imaging sessions. More specifically, body temperature was controlled with a rectal probe connected to a heating pad (TCAT-2DF; Physitemp, Clifton, USA) set at 37 °C and maintained between 36.6 °C and 37.2 °C. Respiratory and heart rate were monitored using a PowerLab data acquisition system with the LabChart software (AD Instruments, USA) and maintained in range between 105–160 breaths/min (mean 130) and 200–320 heart beats/min (mean 250), respectively. Additional IP ketamine/xylazine doses (25 and 2.5 mg/kg, respectively) were infused intermittently (every 90–120 min), as deemed necessary based on changes in physiological parameters. Microbubbles (SonoVue; Bracco Imaging, Switzerland) were injected at a concentration of 8 µl/ml in the tail vein using a 27 G × 1/2″ winged infusion set (DREXCO Medical, France). Three bolus of 80 µL were successively injected.

RC-ULM acquisition

Row–column array, sequences, and scanner

A piezoelectric 15 MHz RCA (Icoprime 4D RCA; Iconeus, Paris, France) consisting of 80 rows and 80 columns with a 0.110-mm pitch was used to perform real-time 3D power Doppler imaging and subsequent 3D ULM acquisitions with microbubbles. The array was connected to a 256-channel functional ultrasound scanner (Iconeus One; Iconeus, Paris, France) driving the probe at 12.5 MHz.

For ultrafast power Doppler, the imaging sequence was defined by transmitting 40 plane waves interleaved on both rows and columns at a PRF of 20 kHz (angles listed in Supplementary Table 1). A four-half-cycle pulse with a duty cycle of 67% and a driving voltage of 25 V was used with a volume rate of 500 Hz. Each block of 200 images (0.4 s) was acquired at 1 Hz to limit the probe heating.

The probe was positioned perpendicular to the animal head’s surface, and a four-axis motor module (Iconeus, Paris, France) was used to ensure proper positioning of the probe using a real-time 3-view power Doppler imaging mode.

The ultrafast RC-ULM sequence was defined by transmitting 28 plane waves interleaved on both rows and columns at a PRF of 28 kHz (angles listed in Supplementary Table 1) yielding a high-volume rate of 1000 Hz. A two-half-cycle pulse was sent at a frequency of 12.5 MHz and a driving voltage of 20 V. The ultrafast sequence comprised 500 images acquired in 500 ms, repeated 400 times. The acquisition time was 10 min, with 1.5-s intervals between each. The three bolus were injected before (at T = −60 s) and after the start of the acquisition at T = 180 and 400 s.

Beamforming and XDoppler reconstruction

RCA beamforming was implemented on two GPUs (A6000 boards, Nvidia, USA) on the scanner and optimized to reach real-time 3D beamforming of the ultrafast volumes including GPU-based clutter filtering. This allowed to perform both the live power Doppler volumetric imaging with a three-view rendering for the positioning but also perform the online reconstruction of the microbubble volumes using the XDoppler scheme ³⁰ down to the prelocalization of microbubbles.

More specifically, the volume RC (transmission along rows and reception along columns) was beamformed using a conventional delay-and-sum beamformer optimized for the dual GPUs (Figure 1). The beamformer grid was set to (λ × λ × λ∕2) and the different angles transmitted along the row were coherently compounded. A SVD clutter filter (cutoff = 15) was then applied on this stack with an ensemble length of 500 frames. The same processing was applied to the CR volume (transmission along columns and reception along rows). The filtered RC and CR stacks were then cross correlated with a window size of five frames for ensemble length using the XDoppler scheme. This correlation approach allows to reduce the signal from the lobes, compared to a synthetic summation approach like orthogonal plane wave compounding (OPW). ³⁰

Figure 1.

3D ULM acquisition and processing workflow using the RCA probe: (i) acquisition setup with the Iconeus One system driving a 160-element RCA probe positioned over the anesthetized mouse brain, with the skull and skin intact, (ii) image formation using the RCA probe, where columns (blue) and rows (red) are used for plane-wave imaging. Beamforming is applied to both RC and CR transmissions, generating compounded volumes. The XDoppler strategy combines these volumes after clutter filtering, ensuring an isotropic PSF, and (iii) 3D ULM processing pipeline: the ultrafast B-mode movie of MB undergoes peak detection, interpolation, Gaussian filtering, and parabolic fitting for tracking. Tracks are created, smoothed, realigned, and rasterized, followed by vesselness filtering and skeletonization to refine the final output.

Within each ultrafast block, MB were prelocalized in real time using a local maximum algorithm and a neighborhood of (5 × 5 × 5) voxels around each maximum was extracted and stored on disk for further processing and refined localization. All those steps were implemented on the GPU in CUDA to run in real time.

3D ULM processing

The localizations were then refined offline using a 3D interpolation with a gaussian kernel of each neighborhood volumes. The sub voxel localization was performed using a 3D paraboloid fit, the resulting position table was stored and a tracking algorithm was used to pair microbubbles from one frame to another within given physical constraints on pairing distance and velocity ranges using a Hungarian linker. ³¹ The tracks were subsequently rasterized to generate microbubble (MB) count volumes and MB velocity component volumes (10 μm isotropic) (Figure 2).

Figure 2.

3D ULM visualization: (i) sagittal MIP of the 3D MB density volume (i, a) and the 3D signed MB velocity volume (i, b), rendered with Amira software. ULM tracks are visualized in 2D slices using IcoStudio software, (ii) axial slices of the 3D MB density (ii, a) and signed MB velocity (ii, b) volumes (slice thickness: 2 mm). Scale bars represent 1 mm, (iii) coronal slices of the 3D MB density (iii, a) and signed MB velocity (iii, b) volumes (slice thickness: 1.5 mm), and (iv) sagittal slices of the 3D MB density (iv, a) and signed MB velocity (iv, b) volumes (slice thickness: 1.5 mm).

Comparison with multiplexed fully populated matrix arrays (MUX–FPM)

We compared the RC-ULM scheme with a MUX–FPM (32 × 32 elements, 0.3 mm pitch, 15 MHz; Vermon, Tours, France) on the same animal. After the RC-ULM acquisition, a pause of 30 min allowed the mouse to recover and to clear the remaining MB from the circulation. During this time, the RCA probe was replaced with the MUX–FPM array. A second imaging session of the same acquisition duration (400 blocks of 0.5 s each acquired in 600 s) was then performed using the same ultrasound frequency and the same voltage. The imaging sequence was adapted from the work of Chavignon et al. ²⁰ where 10 transmit/receive combinations were repeated per tilted plane wave. The PRF was fixed at 12 kHz and six tilted plane waves were transmitted using the following angular sequence: (−5°, 0°); (0°, 0°); (0°, 5°); (0°, −5°); (0°, 0°); and (0°, 5°). The resulting framerate of 200 Hz allowed the accumulation of 100 compounded frames/0.5 s block. The SVD cutoff for the clutter filter was set to 10/100 after optimization. After localization and tracking, we generated volumetric density and velocity maps (resolution 10 µm). The resulting ULM volume was realigned to the RCA volume automatically using the transformation found by intensity-based registration between the two corresponding power Doppler volumes. The MUX–FPM ULM density map was resampled in the RCA space for side-by-side visual comparisons (see Figure 3).

Figure 3.

Comparison between the RCA and the MUX–FPM in the same mouse, following identical injection protocols: (i) registered MB density volumes: The RC-ULM map (blue, right) reveals a greater number of vessels and higher MB density compared to the MUX–FPM map (red, left), (ii) cumulative tracked distance over time: the RCA acquisition demonstrates a faster increase in cumulative tracked distance, indicating more efficient tracking compared to MUX–FPM, and (iii) particle velocity distributions: The RCA acquisition shows a broader and higher count of MB detections across the entire velocity range, particularly for higher velocities, compared to MUX–FPM.

Spatial resolution assessment

To evaluate the spatial resolution of the proposed vascular imaging modality, we investigated the capacity to separate the MB velocity profiles from small and close cortical vessels in the mouse brain (Figure 4). The velocity profile selection was performed on two penetrating arterioles taken at three different depths to intersect different branches of the same penetrating arteriole. Across each profile, the MB velocity was extracted using the Line Probe function of Amira Avizo, using a sampling step of 5 microns.

Figure 4.

Spatial resolution assessment: characterization of MB velocity at various branching levels of cortical penetrating arterioles: (i) zoomed view of the mouse brain cortex (corresponding to the dashed green box in Figure 2(i)), (ii) two arterioles extracted from a region of interest (4λ × 4λ × 10λ) in the ACA territory (white box in (i)), (iii) transverse velocity profiles extracted along dashed colored lines for three depths: at z = 325 µm (green), arteriole α and arteriole β profiles are 137 µm apart. At z = 635 µm (blue), arteriole α splits into branch γ and branch δ, separated by 53 µm. At z = 895 µm (red), branch δ further branches into branch ε and branch φ, separated by 53 µm, and (iv) FRC indicates a frequency-based spatial resolution of 20.26 µm.

Another conventional assessment of the resolution was performed by measuring the Fourier shell correlation (FSC) between two MB count volumes, that were obtained by rasterizing two independent sub-groups of MB tracks split from the same acquisition.

Microvascular graph reconstruction and analysis

Graph construction from ULM maps

As a first step, the MB density volumes were filtered using a hessian-like 3D spatial filter to enhance vessels (3D Jerman enhancement filter ³² ). The resulting volume was then thresholded and skeletonized to extract the vessel centerline coordinates. Radius estimation was performed on these centerline coordinates using the VesselVio ³³ toolbox. The binary vesselness mask and the corresponding skeleton are represented in Supplementary Figures 3 and 5. VesselVio allowed to build a graph structure representing the vascular tree, described as a list of edges (describing a vascular segment) connecting vertices with a connectivity of two between branching nodes (vertices with a connectivity of at least 3). Vascular segment morphological features such as radius, length, surface, volume and tortuosity features such as MB density and 3D MB velocity vector were extracted and stored in the graph data. The radius and centerline MB velocity estimations were used to further estimate the flowrate in each segment of the vascular graph, as proposed in our previous work ¹⁴ assuming a Poiseuille flow applying the following formula:

$$\begin{array}{l}Q=\frac{V_{max}}{2}\pi R^2\\ (Q:flowrate,V_{max}:maximumvelocity,R:radius)\end{array}$$

Graph orientation from local blood flow direction

To orient the graph of vascular segments, we took advantage of the locally measured flow direction and the graph geometry (Figure 5(i)). The orientation of each edge followed the direction of the flow, meaning that the sign of the scalar product between the average MB velocity vector field and the orientation vector of each segment was always positive. By running the Breadth First Search (BFS) algorithm (igraph, python), perfused or drained territories could be easily revealed from any seed segments. A graph traversal algorithm was also implemented in order to reveal the tree structure of any cluster of connected components and store the branching level of each segment in this tree.

Figure 5.

Mapping vascular territories and classifying artery-like or vein-like segments using graph-based analysis: (i) graph creation and orientation: edges represent vascular segments, vertices represent branching points, and orientation follows the 3D velocity vector field. Vessels are classified based on branching behavior, with artery-like segments diverging and vein-like segments merging. Morphometric and flowrate parameters are extracted at the segment level, (ii) MCA-perfused territory: a seed upstream of the MCA bifurcation is used, and downstream segments are extracted. Colors indicate branching levels (left) and flowrate (middle). The right plot shows a negative flowrate slope versus branching levels and positive branching differences, reflecting diverging artery-like topology. Measurements were stopped after the 24th branching level, where radius estimation becomes less accurate. The error bars represent the standard deviation of the flowrate at each branching level, and (iii) GV-drained territory: a seed near the GV confluence node is used, and upstream segments are traced. Segments merge at successive levels, with increasing flowrate and negative branching differences, highlighting converging vein-like patterns.

Artery-like and vein-like labeling from graph analysis

Vascular segments were categorized as artery-like or vein-like using a two-step approach. First, seed segments with a flow rate above 3 µL/min were identified as upstream initiators (arteries) or downstream attractors (veins). The flow path and gradient were traced: segments with a negative gradient were labeled artery-like, and those with a positive gradient were labeled vein-like. Second, to confirm and complement the first step, the branching pattern at each node was analyzed. Nodes with more descendants than ancestors (diverging) were classified as artery-like, while nodes with more ancestors than descendants (merging) were classified as vein-like.

Alignment to Allen Brain Atlas and region labeling

The geometric transformation aligning the brain vascular volume to the Allen Mouse Brain Atlas ³⁴ was automatically estimated using the Iconeus Brain Positioning System, as described by Nouhoum et al., ³⁵ and implemented in the IcoStudio software (Iconeus, Paris, France). This process involved aligning a power Doppler volume extracted from the first 20 s of the ULM acquisition to a mouse brain power Doppler template already pre-aligned with the Allen Brain Atlas. ³⁴

Vascular segments were then assigned to different anatomical regions of interest based on the registration. We chose 10 regions in the study (see Figure 6(ii)). For each chosen region, we tested if the coordinates of the graph centerline points were inside or outside the closed mesh of said region using intriangulation. ³⁶

Figure 6.

Automatic extraction of flowrate and morphometric metrics in brain anatomical regions: (i) the top section illustrates the validation of Murray’s law in vein-like (blue) and artery-like (red) subgraphs, extracted based on branching behavior and color-coded by signed flowrate. A 3D visualization of the brain vasculature highlights vessel radius, while scatter plots of vessel radius versus signed flowrate for veins and arteries confirm adherence to Murray’s law, with coefficients of 2.48 (R²: 0.805) and 2.69 (R²: 0.828), respectively and (ii) atlas-based segmentation of the brain, with 2D histograms showing the density and distribution of flowrate data, illustrating the relationship between vessel size and flow characteristics in each brain region. Symmetrical distributions with balanced densities of arterioles and venules are observed in the isocortex, cortical subplate, thalamus, striatum, and midbrain. Asymmetrical distributions, with a predominance of either arterioles or venules, are found in the olfactory areas, hippocampal formation, pallidum, and hypothalamus. High densities of small venules are prominent in the hippocampal formation, olfactory areas, and hypothalamus, while small arterioles are more frequent in the olfactory areas, pallidum, and striatum. Larger arteries are more prominent in the olfactory areas, hippocampal formation, pallidum, and midbrain, whereas larger veins are particularly noticeable in the hippocampal formation and striatum.

Vascular metrics extraction

Vascular segments were classified according to vascular type (artery-like or vein-like) and anatomical regions. We then extracted in vivo metrics encompassing both morphology and flow-related parameters. The metrics ³⁷ included radius, length, tortuosity, microbubble (MB) velocity and flow rate. We further estimated the flowrate-radius relationship in log scale and fitted the values for the whole brain to compare to Murray’s law ³⁸ and in the different regions of the brain.

Measurements dynamic range

Finally, to assess the dynamic range of the measurements, we plotted the histograms of the velocity and radius distribution of the whole graph, for each independent vascular segment of the graph (Supplementary Figure 1) and traced the minimal and maximal reported values.

Data visualization

The 3D renderings of the MB density and signed velocity volumes (Figures 2 and 4) were performed using the Amira software (Amira v.6.0.1 software; Visualization Sciences Group, USA). The 2D slices represented in the Figure 2 were reconstructed with a prototype of IcoStudio software featuring 3D ULM tracks visualization (Iconeus, Paris, France). The 3D vascular graph renderings were rendered using PyVista visualization toolkit (Figures 5, 6(i), and 7(i)) or using MATLAB (Figures 6(ii) and 7(ii)). Flowrate for artery-like subgraph was rendered in red whereas it was rendered in blue for vein-like subgraphs (Figures 6 and 7).

Figure 7.

Validation of vessel labeling in different brain regions using axial velocity sign and branching classification: (i) three regions are presented: cortical penetrating arterioles and venules (left), cortical vessels around the MCA (middle), and subcortical vessels including the GV and the AZc (right): (top row) vessels are colored based on the sign of their axial velocity (red for downward flow and blue for upward flow). In the cortical penetrating arterioles and venules region, this method separates arterioles and venules due to the anatomy of these cortical vessels. In the other regions, vessel orientation is more complex, and axial flow direction cannot be used to separate arterioles and venules, (bottom row) vessels are labeled based on branching behavior, with artery-like vessels in red and vein-like vessels in blue. The branching classification correctly identifies arterioles and venules, matching the axial velocity sign method with 85% accuracy in the first region and successfully recovering the correct classification of the MCA, GV, and AZc in the other regions. The white arrow in (a, middle) shows that while the sign is flipped when an arteriole starts to penetrate the cortex, it is consistently labeled as an artery at the pial level in the second row. The white arrows in (b, right) shows a similar behavior. For white arrows (c and d, right), we can see that ascending splitting vessels appear in blue in the first row, and in red in the second row, and that descending merging vessels appear in red in the first row, and in blue in the second row and (ii) density of penetrating arterioles and venules in the somatosensory cortex (axial slice). A density of 21 arterioles/mm² and 18.6 venules/mm² was found using the axial velocity and a density of 20.9 arterioles/mm² and 19.1 venules/mm² using the branching-based labeling.

Results

3D RC-ULM imaging

The proposed RC-ULM approach allows for a precise localization and tracking of microbubbles, enabling the generation of detailed microbubble density and velocity volumes that are displayed in 3D and in 2D slices in Figure 2.

We compared those results to the matrix-array MUX-based acquisition on the same mouse with the same injection protocol. Results shown in Figure 3 illustrate the gain in sensitivity and the approximately three times faster cumulative tracked distance with the RC-ULM acquisition. Additionally, RC-ULM was able to capture and track faster microbubbles up to 120 versus 75 mm/s for the MUX-based acquisition.

Using the RCA, important arteries and veins were identified visually and are indicated in Figure 2. Thanks to the high temporal resolution, we successfully reconstructed the arteries of the circle of Willis (CoW) and high-velocity arteries such as the internal cerebral artery (ICA) or the posterior cerebral artery (PCA), where microbubble (MB) velocities reached up to 80 and 56 mm/s, respectively (Figure 2(ii, a) and (iv, b)). For the ICA, our findings are in great agreement with previous phase contrast MRI studies. ³⁹ In larger arterial branches, MB velocities ranged between 36–60 mm/s for the anterior cerebral artery (ACA) and 37–75 mm/s for the middle cerebral artery (MCA; see Figure 2(ii, b), (iii, a), and (iii, b)). Venous structures, including the Galeno vein (GV) and superior sagittal sinus (SSS), exhibited velocities of 27 mm/s and 35 mm/s, respectively. Meanwhile, lower flow velocities (1–10 mm/s) were measured in small penetrating arterioles originating from the MCA and ACA.

Spatial resolution and velocity profiles assessment

To illustrate the ability to measure the flowrate in small cortical vessels, we also selected one region of interest (ROI) around one arteriole and extracted the velocity profile at three different depths, chosen at three different levels of branching. The region of interest and the corresponding velocity profiles are shown in Figure 4(iii). These velocity profiles demonstrated the ability to estimate the velocity in arterial subbranches, separated from less than half the acoustic wavelength. The FWHM and velocity amplitude peaks for each vessel crossed by each profile are reported in the legend. The minimum FWHM among the three examined profiles was 22.9 µm.

The spatial resolution of the RC-ULM volumes was also estimated at 20.26 µm, using the Fourier shell correlation and a 0.5 resolution cutoff (Figure 4(iv)).

Flow-directed graph reconstruction

To illustrate and investigate our flow-directed graph analysis and labeling method, we first focused on two vascular territories of interest (Figure 5(ii) and (iii)). By visually comparing our data with other microvascular volumes from the ex-vivo mouse brain vascular imaging literature,^40–42 we selected a seed segment just before the bifurcation node of the medial cerebral artery (MCA) into three sub-branches, and another seed segment near to the confluence node of the two Galeno veins (GV). These two seed segments were used to initiate the BFS algorithm, that revealed 353 descendants downstream of the MCA seed segment, delineating its perfused vascular territory, and 53 ancestor segments upstream to the GV seed segment, defining its drained vascular territory. Each resulting subgraph is presented in two different representations in the Figure 5(ii) and (iii), the first one displaying the branching level of each vascular segment with respect to the seed position, and the second one showing the corresponding segment flowrate. While the flowrate decreases in the case of the MCA sub-graph, it increases in the case of the GV sub-graph, and the same tendency is also verified in the plots representing the evolution of the average flowrate by branching level, with an average flowrate gradient of −1.35 µL/min/mm in the MCA sub-graph case, and an average flowrate of +0.22 µL/min/mm in the MCA sub-graph case.

Artery-like and vein-like labeling

Branching behavior and labeling accuracy

The average branching difference was consistently positive in the arterial subgraph (e.g. the middle cerebral artery (MCA)) and negative in the venous subgraph (e.g. the Galenic vein (GV)). This reflects the diverging behavior of arteries and the merging behavior of veins along the vascular path. By combining these metrics across the complete graph, 90% of the vascular segments were successfully labeled. Of these, 5% of the labels were automatically corrected to prevent biologically implausible vein-to-artery transitions, while the remaining 10% were discarded due to ambiguity. The resulting artery-like and vein-like subgraphs are visualized in red and blue, respectively, in Figure 6.

Validation using axial velocity

To validate our vessel labeling method, we focused on a cortical region of interest where penetrating arterioles and venules can be distinguished based on the sign of their axial velocity. We compared the vessel labels derived from our branching categorization algorithm with those inferred from axial velocity alone. This comparison yielded an 85% matching score between the two approaches in the cortex, confirming the robustness of our method (Figure 7(i)).

In more complex vascular territories such as the MCA, GV, and azygos cerebral vein (AZc) only the branching categorization algorithm successfully recovered the correct vessel label in regions where flow direction alone is insufficient for accurate classification (Figure 7(i)).

Density and spatial distribution

We estimated the density of penetrating arterioles and venules within the first millimeter of the somotasensory cortex in an axial slice and found values of 20.86 arterioles/mm² and 19.06 venules/mm² (Figure 7(ii)).

Automatic features extraction

Validation of Murray’s law

Log–log plots of signed flowrate versus vessel radius (Figure 6(i), bottom right) demonstrate adherence to Murray’s law for both vein-like (coefficient = 2.48, R² = 0.805) and artery-like (coefficient = 2.69, R² = 0.828) vessels. These results confirm the physiological consistency of the segmentation and flow quantification methods.

Region-specific dynamics

Regional log-log plots (Figure 6(ii)) reveal the signature of flowrate as a function of vessel radius across distinct brain regions, including the isocortex, olfactory areas, hippocampus, cortical subplate, striatum, pallidum, thalamus, hypothalamus, and midbrain. The results highlight regional variations in vascular dynamics, providing a comprehensive characterization of cerebrovascular architecture and function.

More specifically, symmetrical distributions of arterioles and venules were observed in the isocortex, cortical subplate, thalamus, striatum, and midbrain. Asymmetrical distributions, with a predominance of either arterioles or venules, were found in the olfactory areas, hippocampal formation, pallidum, and hypothalamus. High densities of small venules were prominent in the hippocampal formation, olfactory areas, and hypothalamus, while small arterioles were more frequent in the olfactory areas, pallidum, and striatum. Larger arteries were more prominent in the olfactory areas, hippocampal formation, pallidum, and midbrain, whereas larger veins were particularly noticeable in the hippocampal formation and striatum.

Average segment tortuosity values across all brain regions were tightly clustered between 1.1 and 1.2 (Supplementary Table 2). The hypothalamus exhibits the highest tortuosity (1.2), while the olfactory areas and pallidum show the lowest (1.1).

Dynamic range of measurements

To characterize the dynamic range of our measurements, we estimated the distribution of each segment velocities and radii within the complete reconstructed vascular graph (Supplementary Figure 1). Our analysis revealed a minimum detectable velocity of 0.4 mm/s and a minimum resolvable vessel radius of 5 µm, which aligns with the ULM rasterization grid. These values, averaged at the segment level, define the current lower limits of our technique’s sensitivity.

Discussion

The cerebral microvasculature is central to both neurophysiology and pathology, yet in vivo characterization remains a formidable technical challenge even in preclinical research. ULM has recently emerged as a promising microvascular imaging technique, offering the potential to non-invasively map microvascular networks with high-resolution including local blood velocity measurements. Recent advancements have demonstrated the feasibility of 3D ULM in preclinical models, leveraging matrix array technologies^14,16,19 and more recently multiplexed matrix arrays^20,43 and row–column (RC) arrays^28,29 to reduce channel count requirements. While those works proved 3D ULM in vivo feasibility, including in rodent brains, 3D ULM remains challenging in practice due to the massive raw data-rate throughput from the arrays, intensive reconstruction processing and the complex analysis needed to make sense of the 3D ULM data. Building on these recent advances, we thus developed an optimized RCA-based ULM framework, compatible with 256-channel scanners, that enable high quality imaging, full microvascular graph reconstruction, flow-direction labeling, and region-specific quantification. We demonstrated high-resolution, high-sensitivity imaging of cerebral vessels down to 10 µm and the reconstruction of a flow-directed microvascular graph from RCA-based 3D ULM to extract quantitative morphological and hemodynamic signatures across anatomical regions in the mouse brain.

Our approach is well-adapted to reconstruct a diverse range of vessel diameters (10–200 µm and minimum 10 µm) corresponding to small arterioles and venules and flow velocities (about 1–55 mm/s and minimum 0.4 mm/s) as estimated on vessel segments (Supplementary Figure 1). The smallest vessels in this dataset have diameters above the typical capillary range (~5–8 μm) and thus likely correspond to precapillary arterioles and postcapillary venules rather than capillaries. As such it is not possible to establish a direct link with perfusion parameters such as cerebral blood flow or cerebral blood volume. RC-ULM remains complementary to established optical and anatomical imaging techniques. For instance, optical coherence tomography (OCT) and two-photon microscopy offer excellent spatial resolution and flow measurement capabilities but are limited by their penetration depth. ⁴⁴ Conversely, micro-CT ⁹ provides a high field of view and spatial resolution but lacks sensitivity to flow. Other widely used techniques, such as blood oxygen level-dependent (BOLD) fMRI ⁴⁵ and power Doppler ultrasound, ⁴⁶ are unable to resolve most individual arterioles or venules due to their limited spatial resolution but provide indirect sensitivity to cerebral blood flow (CBF) and cerebral blood volume (CBV).

We further registered the reconstructed vessel volumes to the Allen Brain Atlas using a preregistered mouse brain Doppler template. The high contrast of the 3D RC-ULM data then allowed us to automatically reconstruct the vascular graph using standard optical imaging tools. ³³ To categorize vessels as “artery-like” or “vein-like,” we assessed the branching index along the vascular path, rather than relying solely on upward and downward flow direction, which can be ambiguous. By incorporating blood velocity direction into the vascular graph, we oriented each segment and separated artery-like and vein-like subgraphs based on their branching or merging behavior along the flow. This classification was validated in the cortex, where the direction of vertical venules (upward) and arterioles (downward) is well established, allowing us to cross-check our labeling with more than 85% match. When counting the density of penetrating arterioles and venules in the somatosensory cortex (Figure 7(ii)), we found values of 20.86 arterioles/mm² and 19.06 venules/mm², values close to the one found in the literature^47,48 (~16.1 arterioles/mm² and ~23.8 venules/mm²), although authors generally report a larger ratio of venules compared to arterioles in the mouse cortex ⁴⁵ which might point to a lower sensitivity toward venules due to lower velocities.

We could quantify essential microvascular parameters such as flow velocity, vessel diameter, vessel tortuosity, and flow rate for the different vascular scales and regions. The analysis of segment tortuosity across brain regions in naive mice (Supplementary Table 2) revealed a narrow range of values, spanning 1.1–1.2. consistent with literature values. ³³

We further hypothesize that adding such dynamic flow information into directed microvascular graphs will further enable researchers to introduce robust and scalable tools for brain microvascular studies beyond raw ULM data volumes. In the present study, measurements were integrated over the entire duration of the acquisition without accounting for pulsatility. However, the methodology could be further expanded to recover pulsatile flow dynamics by either triggering acquisitions using an ECG or by applying post-processing techniques to extract pulsatility from the existing dataset.^49,50 Such analyses would provide deeper insights into the temporal characteristics of cerebrovascular flow as possible with other techniques ⁵¹ and are an important direction for future work.

The preclinical applications of RC-ULM extend from neurovascular research to drug discovery. It could enable high-throughput microvascular phenotyping in genetic mouse models, quantification of drug-induced vascular changes, and integration of local perfusion data with molecular profiles. RC-ULM could also provide localized microvascular assessment in acute neurovascular disorders (stroke, aneurysms), neurodegenerative diseases (Alzheimer’s, Parkinson’s), tumor angiogenesis, and cerebrovascular small vessel disease.

The translational potential of this technology for clinical investigations of microvasculature is also significant, especially in perioperative settings where a skull piece is removed during the surgery. Those applications include assessing tumoral vascular network in the brain to help surgeons better define the tumor margins and guide the resection or tools in deep and small brain regions without hitting larger vessels.

A key benefit of the RCA technology is its capacity to deliver high-quality 3D ULM using fewer channels (<256) as opposed to traditional matrix arrays that require specialized research scanner with more than a thousand electronic channels. This simplification not only reduces the complexity but also increases the technology’s accessibility and adaptability for broader applications. Compared to a multiplexed array, also the RCA does not require any additional electronics and can be plugged directly to an ultrafast scanner with a smaller cable. Furthermore, using a multiplexor requires to transmit and to receive the channels in different steps which leads to a strong diminution of the volume rate by a factor up to 16 when using a 1:4 multiplexor. In our own tests, RC-ULM consistently outperformed the multiplexed matrix array of similar frequency in tracking efficiency based on cumulative tracking distance and in imaging quality, especially in high flow regions like the CoW (Figure 3). The RC-ULM approach allowed to track higher velocity microbubbles due to its higher volume rate (Figure 3).

There remain limitations to the RC-ULM approach. Firstly, the current spatial resolution is insufficient to assess capillaries, limiting mapping to small arterioles and venules as seen in the velocity and radius histograms (Supplementary Figure 1). This reflects the limited sensitivity of clutter filters to separate very slow or stationary microbubbles from tissue. Improving the clutter filtering and especially SVD clutter filters through multi-thresholds approaches could lower the detection threshold for slower flows. ⁵² Other research groups have proposed alternative methods, including exploiting the non-linear response of microbubbles to detect them even when stationary, ²⁵ or retracing the slowest, undetected parts of their paths by reconnecting precapillary and postcapillary trajectories. ⁵³

Secondly, high microbubble concentrations can cause overlap, compromising detection and tracking, particularly in larger vessels. To mitigate these effects, more robust separation algorithms have been proposed, enabling more efficient microbubble tracking through advanced Kalman filter ⁵⁴ or statistical features. ⁵⁵

Thirdly, the proposed estimation for the flowrate assumes laminar Newtonian flow. Because blood is non-Newtonian, this assumption induces a systematic overestimation of the flow rate and accuracy could be improved by estimating the mean velocity directly over a local vessel cross-section.

While many steps from beamforming to prelocalization could be implemented in real time on the scanner GPUs in this study, performing 3D ULM including tracking and motion correction remains computationally intensive. To address these challenges, we are exploring the development of deep learning algorithms and online GPU-based tracking algorithms.

Another limitation is the restricted “en face” field of view of the RCA due to intrinsic inability to steer the row and column beams on the same side. To overcome this limitation especially for clinical applications such as cardiac imaging, researchers have proposed the use of diverging lens^56–58 or curved RCA arrays. ⁵⁹ Finally, another limitation of the RC-ULM approach is its reduced flexibility in correcting skull-induced aberrations. While individual row and column corrections is possible along each dimension, full 2D adaptive correction as with dense matrix arrays remain impossible. In our case, we observed reduced imaging quality below the main skull sutures, it remains to be understood if this is due to wavefront distortions from said structures and if they could be further corrected to improve imaging in those areas.

Overall, RC-ULM with directed graph reconstruction demonstrates robust capability for high-resolution imaging and whole-brain microvascular quantification in mice compatible with preclinical functional ultrasound systems. This framework could enable large-scale acquisition and analysis of 3D ULM, facilitating comprehensive microvascular flow quantification down to small arterioles and venules for preclinical neurovascular research and complementing existing imaging techniques.

Supplemental Material

sj-docx-1-jcb-10.1177_0271678X261438569 – Supplemental material for In vivo microvascular flow quantification in the mouse brain using row–column ultrasound localization microscopy and directed graph analysis

Supplemental material, sj-docx-1-jcb-10.1177_0271678X261438569 for In vivo microvascular flow quantification in the mouse brain using row–column ultrasound localization microscopy and directed graph analysis by Adrien Bertolo, Jeremy Ferrier, Oscar Demeulenaere, Alexandre Dizeux, Tanguy Delaporte, Bruno Osmanski, Mickael Tanter, Mathieu Pernot and Thomas Deffieux in Journal of Cerebral Blood Flow & Metabolism