Functional ultrasound speckle decorrelation-based velocimetry of the brain

We present a high-speed, contrast free, quantitative ultrasound velocimetry (vUS) for blood flow velocity imaging of the whole coronal section of rodent brain, complementing the high-speed, non-quantitative functional ultrasound imaging (fUS) and the low-speed, quantitative velocimetry by ultrasound localization microscopy (vULM). We developed the theory for analyzing the normalized first order temporal autocorrelation function of ultrasound speckle dynamics to quantify ‘angle-independent’ blood flow velocity. Further, by utilizing an ULM spatial constraint on the bulk motion rejected data, vUS provides high resolution blood flow velocity images at a frame rate of 1 frame/s, compared to ∼ 2 mins per frame with vULM. vUS was validated with numerical simulation, phantom experiments, and in vivo measurements against vULM. We demonstrated the functional imaging ability of vUS by monitoring blood flow velocity changes during whisker stimulation in awake mice.


Introduction
Functional quantitative in vivo imaging of the entire brain with high spatial and temporal resolution remains an open quest in biomedical imaging. Current available methods are limited either by shallow penetration such as two photon microscopy, optical coherence tomography and photoacoustic microscopy that only allow imaging of superficial cortical layers, or by low spatiotemporal resolution of functional magnetic resonance imaging or positron emission tomography. The emerging technology of ultrafast ultrasound 1 paves the way for ultrasound to be applied for functional brain imaging, promising to fulfill the unmet demands of imaging the cerebral hemodynamics of the entire rodent brain with 10-100 resolution, and even holds the promise of measuring neuronal activity with the advent of acoustic reporter genes 2 .
Since the introduction of Power Doppler-based functional ultrasound imaging (fUS) 3 , an increasing number of studies are exploiting the truly impressive capabilities of fUS for functional brain imaging studies [4][5][6][7][8] . However, the exact relationship between the fUS signal and the underlying physiological parameters is not quantitatively established. Although the fUS signal is dominated by blood volume 9 , it's also a mixture of blood flow velocity and the reflectivity of moving particles and it does not provide an unambiguous signal of any of these hemodynamic parameters. On the other hand, ultrasound Color Doppler is able to measure the axial blood flow velocity but suffers from unstable estimations of mean speed due to the presence of noise 3,10,11 and from incorrect estimation if opposite flows exist within the measurement voxel. The microbubble tracking-based ultrasound localization microscopy (ULM 12 ) is able to map with superresolution the whole mouse whole brain vasculature (coronal plane) and quantify the inplane blood flow velocity (vULM 12,13 ) with ~10 resolution. However, it suffers from a fundamental limitation of low temporal resolution as it requires extended data acquisition periods (~150 seconds for 75,000 images 12 ) to accumulate sufficient microbubble events to form a single vascular image and corresponding velocity map, hindering its potential for functional brain imaging studies.
Here, we report a novel ultrasound speckle decorrelation-based velocimetry (vUS) method as a quantitative alterative to fUS and a faster velocimetry method than vULM. In contrast to speckle correlation-based ultrasound imaging velocimetry (UIV) 14,15 and Bmode intensity signal-based speckle decorrelation method (SDC) 16,17 , vUS measures the 'angle-independent' directional blood flow velocity by quantifying the decorrelation of the normalized first order temporal autocorrelation function ( 1 ( )) of the complex ultrasound quadrature signal (sIQ). We developed a comprehensive ultrasound-based 1 ( ) analysis theory and data processing methodology for cerebral blood flow velocity measurement. We validated vUS with numerical simulations, phantom experiments, and in vivo measurements compared with ULM velocimetry (vULM 12,13 ). We demonstrated the functional imaging ability of vUS by quantifying blood flow velocity changes during whisker stimulation in awake head restrained mice.

Theory
The complex ultrasound quadrature signal from a measurement voxel Using a Gaussian shaped pulse wave, the complex point spread function of an ultrasound imaging system is approximated by, where, ( 0 , 0 , 0 ) is the central position of the measurement voxel; , , and are the Gaussian profile width at the 1/e value of the maximum intensity of the point spread function (PSF) in x, y, and z directions, respectively; and 0 is the wave number of the central frequency of the transducer. The time varying ultrasound signal detected from a measurement voxel at time t can be considered as the integration of all moving point scatters within the voxel, as shown in Eq. 2.
where, sIQ is the complex ultrasound quadrature signal of the moving particles of the voxel; is the index of the ℎ scatter; is the total number of scatters within the voxel; R is the reflection factor; and ( , , ) is the position of the scatter. In Eq. 2, we assumed that all scatter points have the same reflection factor R and ( − ( ), − ( ), − ( ))= ( − ( ), − ( ), − ( )), which is reasonable when we are considering the red blood cells (RBCs, [2][3][4][5] in diameter) as the primary moving 3/32 acoustic scatter, which are much smaller than the ultrasound wavelength (usually ~100 ). The ultrasound pressure arising from a given voxel can thus be written as, 1 ( ) of a flow with uniform velocity As shown in Fig. 1a, the movement of particles will cause the detected ultrasound signal to fluctuate in both magnitude and phase. This movement can be quantified based on the dynamic analysis theory of the normalized first-order field autocorrelation function ( 1 ( )).
1 ( ) of a time varying ultrasound signal for a measurement voxel is given by, where, is the time lag; E[…] indicates the average over random initial positions; ⟨… ⟩ represents an ensemble temporal average; and * is the complex conjugate. As illustrated in Fig. 1b of the basic scenario that all scatters move at the same velocity within the ultrasound measurement voxel, the ultrasound pressure arising from the given voxel at time lag can be written as, where, , , and are the flow speed in x, y, and z directions respectively. Given equations (3-5), 1 ( ) for the basic scenario of identically flowing particles can be derived to be, where, we assumed that the positions of the flowing scattering particles are independent of each other. For details regarding the derivation, please refer to the Supplementary Information. From Fig. 1b2 we see that 1 ( ) decays faster for higher speed scattering particles. In addition, for angled flows (5 and 15 ), 1 ( ) rotates and decays to (0, 0) which is due to the axial velocity component inducing a phase shift, as shown in the complex plane (the inset figure). The similar 'rotation path' in the complex plane for the angled flow (5 and 15 ) or transverse flow (9 and 20 ) is due to that they have the same flow direction for the angled and transverse flows. Different flow angles would have different 'rotation paths' in the complex plane as shown in Supplementary Fig. 4b. It is worth noting that the basic scenario, based on the assumption that all scattering particles flow with the same velocity through the given voxel, is valid when the diameter of the vessel containing the flowing particles is much larger than the measurement point spread function (PSF) such that laminar flow speed gradients 18 can be ignored. Therefore, this basic model was used for phantom validation as the plastic tube had an inner diameter of 580 which is ~5 times larger than the US PSF resolution such that the flow velocity could be assumed uniform within the ultrasound measurement voxel.
When the diameter of the vessel is comparable to or smaller than the PSF resolution, we must consider the effect of the laminar flow speed gradients 19,20 . This is true when imaging the cerebral vasculature, as the blood vessel diameter is usually less than 100 as indicated by Fig. 1c1. In this scenario, the group velocity and velocity distribution must be taken into account as the relative movement of the scattering particles will result in additional decorrelation 20 . In Supplementary Fig. 2a, a simulation result is presented to show the effect of a flow speed distribution on the 1 ( ) decorrelation. It indicates that where, ( , , ) is the speed probability distribution function; is the group velocity; and describes the velocity distribution. We can then obtain 1 ( ) for the Gaussian speed distribution with, From our observations, the typical decorrelation time ( ) for blood flow with a speed around 10 mm/s is ~5 ms. Therefore, ↔ 2 2 < 6.25 × 10 −4 2 which is more than 8 times smaller than 4 ↔ 2 ≥ 50 × 10 −4 2 , where '↔' represents the coordinate direction (i.e., x, y or z). Thus, Eq. 9 can be further simplified to, Comparing Eq. 10 with Eq. 6, we note that the axial velocity distribution term also contributes to the signal decorrelation. For details regarding the derivation of this distributed velocity model (DV), please refer to the Supplementary Information.

The mixed ascending & descending flow model
We noticed from the in vivo data that it's common to have opposite flows present in the same measurement voxel, as shown in Fig. 1c1. In this case, 1 ( ) is a mix of dynamics of opposite flows and behaves very differently from that of the single direction flow as can be observed from Figs. 1b2 vs 1c2. Meanwhile, we observed that the majority of the mouse cerebral blood vessels contain an axial velocity component to the flow. This axial flow component causes the frequency spectrum to shift to negative values if the flow is away from the transducer, and positive if the flow is towards the transducer. Thus, we used a positive-negative frequency separation method to calculate 1 ( ) for the positive frequency and negative frequency signals, respectively. With this model, 1 ( ) for the 5/32 positive and negative frequency signals each behaves like that for the single flow direction data as shown in the Supplementary Fig. 2b1-b3, and the fitting is more accurate (FitM2) compared to the no frequency separation model (FitM1) as shown in Fig. 1c2. This frequency separation model is given by, where, and are the complex ultrasound quadrature signal of the negative frequency and positive frequency, respectively; ℱ denotes the Fourier transform; and ℱ −1 denotes the inverse Fourier transform. 1 ( ) and 1 ( ) for and are obtained using Eq. 4, respectively. / , in which 'a' and 't' indicate angled and transverse flow, respectively. The blue diffuse spot in (b1) depicts the size of the ultrasound PSF. The inset in (b2) shows 1 ( ) in the complex plane. (c1) ULM measurement shows that blood vessels are usually smaller than the ultrasound PSF and opposite flows (ascending and descending) are commonly observed within the ultrasound measurement voxel; (c2) The experimental 1 ( ) (green dots) of a voxel having both ascending and descending blood flow vessels, and fitting results using the DV model only (red dashed curve, fitting accuracy R=0.39) and the negative-positive frequency separated DV model (blue solid curve, R=0.89). The white diffuse spot in (c1) depicts the ultrasound PSF. The

6/32
inset shows 1 ( ) in the complex plane. (d) Representative blood flow velocity map of a mouse coronal plane (around Bregma: -2.18 mm) obtained using vUS. Top: vUS image obtained without ULM spatial constraint; bottom: vUS image obtained with the ULM spatial mask that constrained vUS to process pixels only within the blood vessel identified in the ULM image. Descending flow velocity map was overlapped to ascending flow velocity map. Fig. 1d shows a representative blood flow velocity map in the coronal plane of Bregma ~-2.18 mm of a mouse brain. The top figure was obtained with the vUS methods described previously without using a vasculature map measured with ULM (Online Methods). With an initially obtained ULM vascular map, we can further create a spatial vascular mask to constrain the spatial pixels used for vUS processing (Online Methods and Supplementary Fig. 1d1). As shown in the bottom figure of Fig. 1d, the higher spatial resolution was given by the ULM spatial mask, which spatially constrained the vUS to process pixels within the blood vessels identified in the ULM image. This approach has two major benefits: a higher resolution vascular image in combination with higher frame rate, and a faster data processing rate by thresholding out non-vessel pixels. It can be noted that both vUS processing methods are able to quantify blood flow velocity of the whole coronal plane and provide similar results.

Validation of vUS
We first validated vUS using numerical simulation (Online Methods), as shown in Fig.  2a and Supplementary Fig. 4. The reconstructed total velocity ( ), transverse velocity ( ) and axial velocity ( ) agree well with preset speeds and angles. Note that vUS is capable of measuring transverse flows (i.e. = 0°) which is not possible with Color Doppler velocimetry.   Fig. 2b1 shows the velocity map for a transverse flow; and the middle and bottom panels show the velocity map for two angled flows. A laminar velocity profile was observed, particularly for higher speed, as indicated in the inset Fig. 2b1. We note that the vUS measurements agree well with the preset speeds even for speeds as low as 1 mm/s for both transverse and angled flows, as shown in Fig. 2b2. For more details regarding the phantom validation results, please refer to Supplementary Fig. 5.
We further performed in vivo validation by comparing the velocity measured with ultrasound localization microscopy velocimetry (vULM, Online Methods) against vUS, as shown in Fig. 2c and Supplementary Fig. 6. For a fair comparison, both the vULM and the vUS measurements were applied with a spatial mask that ensures nonzero valued pixels for both vUS and vULM measurements. Qualitatively, we noted from the velocity maps (Fig. 2c1&c2) that the measured axial velocity (Fig2. c1) and total velocity (Fig2. c2) agree well between vUS and vULM. Quantitatively, the weighted scatter plots of all nonzero pixels between vUS and vULM in Fig. 2c3 show that the vUS measurement is highly correlated with the vULM measurement. Fig. 2c4 further shows a high correlation of the mean velocity of 50 vessels (Supplementary Fig. 6a1) between vULM and vUS measurements. We also noticed that the vUS measurements of total velocity have a larger absolute value compared to the vULM measurements. In our opinion, this is reasonable as vULM only tracks the in-plane microbubble movement ( and ) while for vUS we used 2D velocities ( and ) to represent 3D dynamics (vUS fitting algorithm, Online Methods) in which the flow across planes (in y direction) also contributes to 1 ( ) decorrelation (Eq. 10). As expected, the measured axial velocity of vUS is very close to that of vULM ( = 0.99 − 0.61) as shown in Fig. 2c3.

Blood flow velocity change evoked by whisker stimulation
It is well known that regional cerebral blood flow/velocity increases in response to taskevoked brain activation 21,22 . To test the functional imaging capability of vUS, we applied vUS to monitor the blood flow velocity in response to whisker stimulation. We developed an animal preparation protocol for chronic ultrasound imaging in awake mice (Online Methods), as shown in Fig. 3a&b. Fig. 3c presents the whisker stimulation protocol which consists with 30 s baseline followed by 10 trials of 15 s stimulation and with a 45 s interstimulus interval. The vUS images were acquired at a frame rate of 1 frame/s and processed based on the ULM spatial mask that acquired at the same coronal plane. Fig. 3d shows the correlation coefficient map between the blood flow velocity measured with vUS and the stimulation pattern. We note that in addition to the significant activation of vessels in the primary somatosensory barrel field (BF), the blood vessel flowing through the posterior complex (PO) and ventral posteromedial nucleus (VPM) of the thalamus also exhibited activation. This result agrees with previous rodent functional studies that mechanoreceptive whisker information reaches the barrel cortex via the thalamic VPM nuclei 23 and the PO is a paralemniscal pathway for whisker signal processing 24 . The velocity time courses and averaged velocity relative change of the 10 trials of vessels V1 and V2 indicate robust blood flow velocity increase in response to the stimulation as shown in Figs. 3e&f. The time course of vessel V3 on the ipsilateral cortex of the stimulation was plotted as a control region, which shows no correlation with the stimulation. Additionally, we have observed increased blood flow in the auditory cortex (AuD); lateral posterior nucleus (LP) which can be explained by the sound transmitted by the Picospritzer and the grooming of the mouse right after the stimulus. We have encountered slightly increased blood flow in the ipsilateral cortex and thalamic areas that can also be attributed to grooming. The Supplemental Video 1 shows the relative blood flow velocity changes before, during and after the stimulation.
We further compared the 1 ( ) for baseline and under stimulation of the same spatial pixel in V1, as shown in Fig. 3g. It is evident that 1 ( ) decays faster when under stimulation compared to that during the baseline, indicative of faster dynamics, i.e. elevated blood flow speed in response to whisker stimulation.

Discussion
The development of robust blood flow velocity measurement technologies has been of great importance in neuroscience research as quantifying blood flow alterations enables the assessment of brain disease [25][26][27] and interpretation of regional neural function according to neurovascular coupling 28 . In this work, we introduced vUS based on ultrasound speckle dynamic analysis to quantify cerebral blood flow velocity with a temporal resolution of 1 frame/s, free of contrast with more than 10 mm penetration and multiscale spatial resolution (with (~ 20 ) or without (~ 100 ) using ULM spatial constraint). vUS provides much deeper penetration compared to optical velocimetry methods which are usually restricted to superficial layers of less than 1 mm depth 29 while maintaining high spatial and temporal resolution compared to magnetic resonance imaging-based phase contrast velocity mapping 30 .
Conventionally, ultrasound speckle is considered as undesirable noise that impacts ultrasound image contrast negatively and a main focus of research has been on the removal of ultrasound speckle 31 . Currently, two ultrasound speckle analysis methods have been studied trying to use ultrasound speckle to resolve flow information, including the cross correlation-based ultrasound imaging velocimetry (UIV) 14,15 and the B-mode intensity-based speckle decorrelation method (SDC) 16,17 . By exploiting dynamic analysis of the normalized first order field autocorrelation function ( 1 ( ) ) of the complex ultrasound quadrature signal, for the first time to the best of our knowledge, we developed the comprehensive vUS theory to accurately measure axial velocity and transverse velocity under different scenarios. The ultrafast ultrasound plane-wave coherent compounding acquisition paves the way for vUS as high frame rate is crucial in resolving speckle decorrelation caused by moving scattering particles. The spatiotemporal singular value decomposition filtering plays an important role in rejecting the bulk motion which enables the decorrelation of 1 ( ) to represent the dynamics of the motion of red blood cells. Depending on the nature of flow within the US voxel, different vUS models are selected to reconstruct the flow velocities. We validated vUS with numerical simulation, phantom experiments, and in vivo measurements against vULM, and demonstrated its capability in brain function studies.
Indeed, vUS is an 'angle-independent' three dimensional (3D) velocimetry. It's also sensitive to the through-plane speed ( ). However, since the PSF resolution in y is more than 3 times larger than that in x (i.e., >3 ) when using the 1D geometry of the acoustic transducer array which results in a more than 9 times smaller decorrelation from compared to , in this work vUS was simplified to estimate in-plane 2D velocities (i.e., and ). This simplification, however, results in an over estimation of the transverse velocity ( ) as tends to compensate for the decorrelation caused by . We note that the over estimation of is negligible and the measured axial velocity ( ) is very accurate as indicated by Fig. 2C3.
In the future, with the development of fast 3D ultrasound imaging technology using a 2D transducer matrix 32 , vUS can be easily adopted for 3D velocimetry of the whole rodent brain. In addition, in combination with oxygenation measurements using multispectral photoacoustic tomography (mPAT) and given the blood flow (velocity) measured by vUS, the direct biomarker indicating neuron activity, i.e., the metabolic rate of oxygen, can be readily measured, which will greatly facilitate neuroscience research.

Plane wave coherent compounding-based data acquisition
The ultrasound signal was acquired with a commercial ultrafast ultrasound imaging system (Vantage 256, Verasonics Inc. Kirkland, WA, USA) and a linear ultrasonic probe (L22-14v, Verasonics Inc. Kirkland, WA, USA). The Vantage 256 system has 256 parallelized emission and receiving channels, and can acquire planar images at a frame rate up to 30 kHz when the imaging depth is ~15 mm. The L22-14v ultrasonic probe has a bandwidth of 12.4 MHz (67%, -6 dB). It has an elevation focus at z=6 mm.
To accurately resolve ultrasound speckle decorrelation, high temporal resolution is mandatory, especially for fast flows. As an example shown in Supplementary Fig. 2b, only 2 time points were captured during the 1 ( ) decay period when using a 1 kHz sampling rate. With so few time points in the autocorrelation function, it is not possible to accurately resolve the dynamics of the flow. On the other hand, to enhance the contrastto-noise ratio, coherent plane-wave compounding with several tilted emission angles 33 are needed which reduces the temporal resolution. Based on the literature 34,35 and our measurements using phase resolved Doppler Optical Coherence Tomography 36 and vULM (Online Methods) that blood flow velocity of a mouse brain is in the range of ~1 mm/s for capillaries to ~40 mm/s for large arteries and ~10-20 mm/s for the majority of middle sized vessels, we selected an image frame of 5 kHz using 5 tilt angles (−6°, −3°, 0°, 3°, 6°) for coherent plane-wave compounding. As illustrated in Supplementary Fig. 1a&b, the ultrasound pulse frame rate was 30 kHz which is mainly limited by the transmit time of the ultrasound signal in the sample through the intended imaging depth; and the 5 tilted plane waves were coherently summed to form a compounded image whose frame rate was set to be 5 kHz. In addition, the total data acquisition time for each vUS frame was selected to be 200 . Thus, each frame of the complex ultrasound quadrature signal (IQ) was obtained by beamforming and coherent compounding the RF data from the 5 emission angles, and there were 1000 frames of IQ data used for a vUS frame calculation.
Ideally, if the acquired radio frequency (RF) data can be transferred to the host computer and saved within 200 the vUS frame rate is 5 frames/s. However, limited by the data transfer rate and data processing (beamforming) time, it may take tens of seconds to save the data for a single vUS frame. To realize a higher frame rate for longer time imaging, we chose to save the RF data directly to an m2 drive with a write rate of ~1 Gbit/s and post process the data offline. With this data acquisition and processing procedure, we are able to perform continuous imaging up to 30 mins with a vUS frame rate of 1 frame/s.

Clutter rejection
Clutter rejection is needed as the background/bulk motion signal is dominant in the ultrasound images which makes it difficult to analyze the 1 ( ) dynamics. Supplementary Fig. 3a presents an example to show the importance of clutter rejection. We note that 1 ( ) calculated from raw IQ data ( 1 ( ) -RAW) has a very small decorrelation (decays form ~0.995 to 0.975) which is due to the fact that 1 ( ) characterizes the signal fluctuation of dynamic components while for ultrasound, particularly for in vivo imaging, the dynamic signal only contributes a very small portion of the whole signal. In addition, the background signal from the "static" tissue is not actually static because of motion induced by physiological fluctuations such as respiration and the cardiac cycle and bulk movement of the tissue. This motion causes decay of the correlation function which confounds proper estimation of the blood flow speeds. As a result, proper clutter rejection methods are needed to filter these confounding signals before estimating the correlation functions.
Exploiting the spatial correlation of the clutter signal, we used a spatiotemporal filtering method (singular value decomposition, SVD, Eq. 13 37 ) to remove the first two (Nc=3) highest 14/32 singular value signal components for the phantom data. To reject the bulk motion signal for the in vivo data, we used a combination of SVD and high pass filtering. The first 20 highest singular value signal components were removed (Nc=21), followed by a fourth order Butterworth high pass filtering with a cutoff frequency of 25 Hz corresponding with a 1 mm/s speed cutoff. With this clutter rejection method, the static/bulk motion signal is removed and the residual signal is the dynamic representative of the moving of red blood cells. 1 ( ) calculated from the dynamic sIQ signal ( 1 ( )-SVDHP) decays from ~0.8 to 0.15 (see Supplementary Fig. 3a), much larger than observed with 1 ( )-RAW.
where, sIQ is the dynamic signal; Nc is the cutoff rank for SVD processing; ( , ) is the spatial singular matrix; is the singular value corresponding with the i th rank; and ( ) is the temporal singular vector.

( ) calculation
A long observation time is preferred for the 1 ( )-based dynamic analysis for sufficient ensemble averaging. However, a long sampling time requires longer data acquisition and results in a large data size. To determine the appropriate observation time, we compared 1 ( ) calculated with different observation times. As shown in Supplementary Fig. 3c1, 1 ( ) has a similar decorrelation for the observation time of = 200 and 400 , but that for = 50 the decorrelation differs. Supplementary Fig. 3c2 compares the similarity R (Eq. 20, R equals to 1 minus the least squares difference of 1 ( ) obtained with different observation times compared to = 500 ) and the noise level, which was the average value of 1 ( = [3: 20 ]), for different observation times. We note that = 200 provides a low noise level and high R. Thus, in this study we used a 200 period of data (i.e. 1000 coherent compounded frames) for the 1 ( ) calculation.
Another parameter in the 1 ( ) analysis is the maximum time lag in the autocorrelation function. Generally speaking, slow flow takes a longer time for 1 ( ) to decay to the same level as for fast flow. As shown in Supplementary Fig. 3d, it takes around 17 ms for 1 ( ) to decay to <0.2 when the speed is 5 mm/s, while it just takes ~ 5 ms to decay to <0.2 when the speed is 15 mm/s. Since the flow velocity of vessels of interest is >5 mm/s, we calculated 1 ( ) with a maximum time lag of 20 .
vUS Fitting algorithm Supplementary Fig. 1d illustrates the vUS data processing algorithm. The first step is to check the signal quality so that only the data satisfying the defined criteria are analyzed to fit the 1 ( ). Specifically, the criteria include the signal-to-noise ratio (Eq. 14), the ratio of positive frequency power to negative frequency power(Eq. 15), the absolute value of 1 ( ) at the first time lag (Eq. 16), and a spatial mask if using the ultrasound localization microscopy data as a spatial constraint (Eq. 17, Supplementary Fig. 1d1) Mask ULM ( , ) > 0 (17) where, ℱ denotes the Fourier transform. These criteria enable us to skip the poor quality data, which also greatly reduces the processing time. Then the sIQ data is split into and using Eq. 12, and the fitting procedure is applied for both and . In practice, random noise results in a prompt 'jump' of 1 (1), i.e. the change of 1 (0) to 1 (1) is not a smooth transition compared to 1 (1) to the end of the decorrelation as the noise is uncorrelated between sequential measurement time points. We therefore modified the 1 ( ) equation by using an 'F' factor to account for this 'jump'. Also, it is worth noting that when using a linear transducer array that the ultrasound PSF is anisotropic in the transverse directions, i.e. ≠ . In our experimental setup, was more than 3 times larger than which results in a more than 9 times slower signal decorrelation rate from compared to that from . Therefore, we omitted the y component from the 1 ( ) fitting to simplify the data processing. Thus the theoretical 1 ( ) model used for fitting out experimental data is, where, > 0; , , and 0 = 2 0 ⁄ are known parameters. With this simplification, we noticed that the vUS measurements of total velocity have a larger absolute value compared to the vULM measurements as discussed in the Validation section. A proper initial guess of the unknown parameters (i.e., F, , , and ) is important to achieve high fitting accuracy and efficiency. The initial guess of 0 was set to be 0 =| 1 (1)|. We used the real part of 1 ( ) to determine 0 by finding the time lag when 1 ( ) reaches the first minimum. Then 0 is obtained using, We tested a mesh of and values to determine the initial guess of 0 and 0 by finding the pair of 0 and 0 that maximizes the coefficient of determination, R. R is defined in Eq. 20 and was also used in the final fitting process as the objective function for a constrained least squares regression non-linear fitting procedure to estimate the values for F, , , and based on the initial guesses.
Representative reconstructed results of axial velocity and total velocity ( = √ 2 + 2 )) are shown in Supplementary Fig. 1e (without a ULM spatial constraint) and Supplementary Fig. 1f (with a ULM spatial constraint).

Power Doppler and Color Doppler calculation
The Power Doppler image (PDI) was calculated as 3 , (21) where, N is the number of samples and sIQ is the complex ultrasound quadrature signal of the moving particles. The axial velocity based on the Color Doppler calculation is obtained with 10 , where, c is the sound speed in the medium and c= 1540 m/s was used in this study; 0 is the transducer center frequency; is the frame rate; and ℱ denotes the Fourier transform.

Ultrasound Localization Microscopy
The ultrasound localization microscopy (ULM) images and the ULM-based velocity maps (vULM) were obtained based on a microbubble tracking and accumulation method described in 12,13 . Briefly, a frame-to-frame subtraction was applied to the IQ data to get the dynamic microbubble signal. The images of the microbubble were rescaled to have a pixel size of 10 × 10 . The centroid position for each microbubble was then identified with 10 precision by deconvolving the system point spread function. By accumulating the centroid positions over time, a high resolution image of the cerebral vasculature image (ULM) is obtained. Further, by identifying and tracking the same microbubble's position, the in-plane flow velocity of the microbubble can be calculated based on the travel distance and the imaging frame rate. The final velocity for coordinates (z, x) consists of descending and ascending flows, and the speed for each direction was obtained by averaging the same directional flow speed at all time points when the absolute value was greater than 0, respectively.

Numerical Simulation
In this study, two dimensional (x-z) flow and ultrasound detection was simulated to validate vUS. Point scattering particles (5 in diameter) were randomly generated at the initialization segment which is outside the ultrasound measurement voxel. Then the flowing positions were calculated for all time points based on the preset flow speed and flow angle at a temporal rate of 5 KHz. The detected ultrasound signal (sIQ) was obtained based on Eq. 2 for each time point. Then the simulated 1 ( ) was calculated according to Eq. 4 with 1000 observation time points (i.e. 200 ms) and 100 autocorrelation calculation time lags (i.e. 20 ms). Flow velocity was then reconstructed by applying vUS processing on the simulated 1 ( ).

Phantom experiment and data processing
For the phantom validation experiment, a plastic micro tube (inner diameter 580 , Intramedic Inc.) was buried in a homemade agar phantom with an angle of ~ 30° (angled flow), and another plastic micro tube was aligned close to ~ 0° (transvers flow) in another homemade agar phantom. A blood solution was pumped through the tubes with a syringe pump (Harvard Apparatus) at speeds of 1, 3,5,7,9,11,13,15,20,25, and 30 mm/s. SVD was performed to filter the background signal clutter by removing the first two highest singular value components. Since the diameter of the tube is much larger than the ultrasound resolution, the red blood cell speed distribution can be considered uniform. Therefore, the basic model (i.e., the SV model) was used for processing the phantom data.