Retrieving Pulsatility in Ultrasound Localization Microscopy

Ultrasound localization microscopy (ULM) is a vascular imaging method that provides a 10-fold improvement in resolution compared to ultrasound Doppler imaging. Because typical ULM acquisitions accumulate large numbers of synthetic microbubble (MB) trajectories over hundreds of cardiac cycles, transient hemodynamic variations such as pulsatility get averaged out. Here we introduce two independent processing methods to retrieve pulsatile flow characteristics from MB trajectories sampled at kilohertz frame rates and demonstrate their potential on a simulated dataset. The first approach follows a Lagrangian description of the flow. We filter the MB trajectories to eliminate ULM localization grid artifacts and successfully recover the pulsatility fraction <inline-formula> <tex-math notation="LaTeX">$P_{\mathrm {f}}$ </tex-math></inline-formula> with a root mean square error (RMSE) of 3.3%. Our second approach follows a Eulerian description of the flow. It relies on the accumulation of MB velocity estimates as observed from a stationary observer. We show that pulsatile flow gives rise to a bimodal velocity distribution with peaks indicating the maximum and minimum velocities of the cardiac cycle. In this second method, we recovered the pulsatility fraction <inline-formula> <tex-math notation="LaTeX">$P_{\mathrm {f}}$ </tex-math></inline-formula> by measuring the location of these distribution peaks with a RMSE of 5.2%. We evaluated the impact of the MB localization precision <inline-formula> <tex-math notation="LaTeX">$\sigma $ </tex-math></inline-formula> on our ability to retrieve the bimodal signature of a pulsatile flow. Together, our results demonstrate that pulsatility can be retrieved from ULM acquisitions at kilohertz frame rate and that the estimation of the pulsatility fraction improves with localization precision.


I. INTRODUCTION
T HE resolution of conventional ultrasound images is limited by the wavelength (λ) dependent Rayleigh diffraction limit [1] while imaging depth is inversely proportional with λ, leading to a fundamental trade-off between image resolution and investigation depth. Recently, this trade-off has been circumvented by the introduction of ULM. This super-resolution vascular imaging technique relies on the localization of sub-wavelength contrast agent [2]. By accumulating thousands of micrometer-sized agent positions with sub-wavelength precision, an image can be reconstructed with a 10-fold resolution improvement compared to conventional ultrasound while retaining the same imaging depth.
To date, ULM relies on synthetic ultrasound contrast agents made of polydisperse gas-filled microbubbles (MBs). Currently, ULM reconstruction algorithms rely on the hypothesis that individual MBs can be localized with subwavelength precision. Individual MBs are typically captured over multiple frames at kilohertz frame rates. The MB positions in subsequent frames can thus be linked together, forming trajectories that reveal the underlying vascular architecture. ULM images are then rendered using either the density of MBs per each pixel, blood velocity estimates, or other VOLUME 2, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ metrics associated with morphological features such as tortuosity [9], vessel orientation [10], or main flow direction [8].
The ULM velocity image displays the average velocity found at each pixel location during the entire ULM acquisition. Because MB trajectories are accumulated over hundreds of cardiac cycles, fast flow fluctuations are averaged out in ULM velocity maps. By accelerating the ULM acquisition, the averaging of velocities is reduced. An intuitive approach to speed up ULM acquisition consists of increasing the MB population density, as localizing larger numbers of MBs per frame reduces the total number of frames accordingly. However, this poses a new imaging trade-off between acquisition time and MB localization precision. Huang et al. [19] proposed to increase the MB sparsity via post-processing methods. They separated MB echoes into subpopulations using a spatiotemporal Fourier filter. Each subpopulation was then individually processed by conventional ULM processing steps. Alternatively, machine learning-based approaches for MB localization were shown to be able to handle higher MB concentrations by disentangling the interference pattern of spatially overlapping MBs [20], [21], [22]. ULM can also be triggered with electrocardiograms to observe phenomenons independently from the cardiac cycle such as pulsatility [23], or avoid big motion [24].
In this study, we investigate the possibility of retrieving pulsatility in the microvasculature by post-processing a conventional ULM dataset. Microangiopathy or small vessel disease, if present in the brain, is shown to be related to stroke and dementia [25]. Shi et al. [26] showed that white matter hyperintensities, which characterize cerebral small vessel disease, are associated with increased intracranial pulsatility. In their study, pulsatility was measured in the major arteries and veins of the brain using phase-contrast MRI. Having the capability to measure pulsatility in smaller vessels will provide insight into the pathogenesis of small vessel disease.
We aim to retrieve pulsatility from conventional ULM datasets without any change in the acquisition pipeline. We hypothesize that the MB trajectories sampled at kilohertz frame rate contain enough information to retrieve pulsatility. We develop two methods capable of retrieving pulsatility from raw ULM data. The first method follows a Lagrangian description of the flow. That is, using velocity estimates of a single MB, it relies on filtering to retrieve the pulsatility fraction in the reconstructed trajectories. The second method follows a Eulerian description of the flow. Specifically, the temporal distribution of velocities found at a fixed location in space is calculated. The pulsatility fraction was extracted from the bimodality of the velocity distribution. We validate both methods on simulated data generated using a custom ULM simulator that was designed to mimic pulsatile flow in rat cortical brain vessels.

II. METHODS
The ULM simulation pipeline is described in Subsection II-A and pulsatility retrieval methods are described in II-B. We provide a list of symbols in Appendix A.

A. ULM SIMULATOR DESIGN
The designed ULM simulator is presented in Fig. 1(a) and is used to study two cases called the 'MB Localization' and 'No MB Localization' scenarios. In the former, the full ULM process is simulated including B-mode images. In the latter MB localization error is directly added to the MB positions. In both scenarios the simulation was performed in 2D.

1) FIRST SIMULATION SCENARIO: LOCALIZATION
The ''Localization'' scenario mimics all steps of an experimental ULM acquisition and consists of three main modules (see Fig. 1(a)): a module simulating the MB positions (in blue), a module simulating ultrasound frames (in green) and a module performing ULM processing (in yellow).
In the first module, the ground truth flow u GT n (r) is simulated by modeling its temporal and spatial characteristics separately, with r the lateral coordinate and n the time index.
The temporal behavior of the flow at the centerline of the vessel (r = 0) is modelled following in vivo observations by Santisakultarm et al. [27] performed in mouse cortical brain vessels (see Fig. 1(b)). It is thus assumed that the shapes of the pulsatile cycles of rat and mice are similar. The pulsatility fraction is defined as whereū GT n (r) denotes the temporal average. To construct the pulsatile flow at the centerline of the vessel u GT n (0) we retrieved the temporal averageū GT n (0), P f and the shape of the pulsatile cycle from the results reported in [27]. The frequency of pulsatility was fixed at 300 bpm. Fig. 1(b) shows the simulation of u GT n (0) for three different values of P f corresponding to a steady-state flow (P f = 0) a pulsatile flow in a venule (P f = 0.2) and in an arteriole (P f = 0.4).
We calculate the flow u GT n (r) for a given temporal sampling n = t/dt, from u GT n (0) following two assumptions on the spatial characteristics of the flow (see Table 2). First, u GT n (r) is assumed to be dependent on the lateral coordinate r only, corresponding to the rigid vessel-hypothesis. Second, the flow velocity profile along r is assumed to be parabolic, corresponding to Poiseuille flow [28]. The parabolic profile describing u GT n (r) along r is constructed from the centerline velocity and enforcing zero velocity at the boundaries. The separate simulation of the temporal and spatial behaviour of the flow was assumed to be appropriate in the microvasculature due to the low Womersley number flow [29], [30].
MB positions during the acquisition were simulated by initializing and propagating a fixed number of point scatterers in the vessel. The axial MB coordinate l was uniformly sampled along the vessel length L. The lateral MB coordinate r was sampled from a distribution of parabolic shape found by normalization of the parabolic velocity profile u GT n (r) such that it represents a probability density function. The MBs were propagated to the next frame by a spatial step given by the local simulated flow velocity multiplied by the timestep dt. An outlet condition was implemented so that once a MB left the vessel, a new one was initialized at the inlet. The lateral coordinate r of this new MB was drawn from the parabolic distribution while the axial coordinate was drawn from a uniform distribution on the interval [0, u GT n (r)dt)]. In the second module, the Verasonics Research Ultrasound Simulator (VRUS, Vantage, Verasonics, Kirkland, WA, USA) was used to simulate the ultrasound B-mode frames at 1kHz [31]. The L22-14vX probe (Vermon, Tours, France) was simulated to transmit a single plane wave at 0 • with pulse duration of 2 cycles and a main frequency of 17.6 MHz. The Verasonics Reconstruction software was used to beamform the Radio-Frequency data. Additional parameters can be found in Appendix A, Table 2. To improve computational efficiency the frames were constructed in two steps: simulation of the MB PSFs and addition of a speckle pattern (see Fig. 1(d)).
The MB response was simulated as a subwavelength scatterer and was simulated as if the MBs were to flow in a homogeneous medium without speckle, resulting in a frame with a few point spread functions (PSF). A speckle pattern was randomly selected from a set of a priori simulated speckle frames and added to the PSF frame. The speckle patterns were simulated by placing 10 5 weaker scatterers of random reflectivity at a uniform random position in the field of view. The signal-to-noise ratio (SNR) of the resulting image can be controlled by setting the intensity of the speckle pattern. An average SNR of 14 dB was simulated.
In the third module, we performed ULM processing on blocks of 1000 simulated frames using the LOTUS toolbox [31]. We optimized ULM parameters based on the visual improvement of the microvessel reconstruction indicated by vessel filling, vessel separation and the presence of a parabolic velocity profile. Their values are provided in Appendix A, Table 3. We applied a SVD filter to the simulated B-mode frames [32] (see Table 3). Since no static tissue is present in our simulation, the SVD filter was solely applied to reduce speckle noise. For the localization of MBs, we used the radial symmetry algorithm, due to its low localization error and low computation time in mid-to-low-SNR scenarios [31], [33]. We calculated MB trajectories using a Kuhn-Munkres assignment [34] and computed the corresponding velocity estimates u n (r) by a backward Euler method [7]. We smoothed MB trajectories using a moving average filter to generate the filtered velocity estimates VOLUME 2, 2022 u f n (r) [31] and rendered these trajectories on a super resolved grid to reconstruct the velocity map. Additionally, a density map was constructed by finding the MB count for each pixel and a flow orientation θ map was constructed, with θ the angle between the x-axis and the direction of flow (in counter clockwise direction).

2) SECOND SIMULATION SCENARIO: NO-LOCALIZATION
In the No-Localization simulation scenario a spatial offset mimicking a MB localization error is added to the ground truth MB positions simulated by the first module (orange section in Fig. 1(a)). The MB positions with added localization errors are then fed to the tracking module ( Fig. 1(e)). The localization errors in x and z direction were assumed independent and identically distributed (i.i.d.) as N (0, σ 2 ), with σ the localization precision. We assume that σ = σ x = σ z . This simulation scenario allows us to investigate the effect of localization precision on retrieval of pulsatility since σ can be varied in simulation. It should be noted that in practice, different localization algorithms lead to different shapes and mean values of localization error distributions [31].

3) SIMULATION CONFIGURATIONS
We defined two distinct microvascular configurations as illustrated in Fig. 1(f), with specific simulation parameters provided in Table 1. For both configurations the simulated flow characteristics were matched to that of in vivo observations [27].
The double vessel configuration mimics the in vivo configuration of cortical brain vessels. Here, the vascular architecture consisted of a straight 30 µm wide descending arteriole and a straight 100 µm wide ascending venule that converge towards each other. We used this configuration in combination with the Localization simulation scenario to evaluate the performance of our ULM simulator.
The second configuration consisted of a single 100 µm wide straight arteriole and was used in combination with the No-Localization simulation scenario. This configuration enables us to assess the performance of the pulsatility retrieval methods for different values of the localization precision σ .

B. METHODS TO RETRIEVE PULSATILITY IN ULM DATA
In this section, we introduce two methods for pulsatility retrieval in ULM data following respectively a Lagrangian and Eulerian description of the blood flow.

1) SINGLE-TRAJECTORY VELOCITY FILTERING
In this Lagrangian method to recover pulsatility, the velocity estimates of individual MB trajectories were filtered. The velocity estimates u n (r) are found by dividing the traveled distance over one frame by the time-step, also known as the backward Euler method [7]. To reduce the effect of MB localization errors on the velocity estimates, we filtered trajectories using a moving average filter defined as where s is the span of the moving average filter and · is the floor operation. Note that the span is expressed in number of frames and it was chosen to be an odd number such that an equal number of velocity estimates before and after frame n are included. It dictates the time window over which the u n (r) is averaged. Its value is therefore dependent on the imaging frame rate. Here, s should be compared to the 200 frames covering one simulated pulsatile cycle (5 Hz, dt = 1 ms). Grid-based artifacts appear in the velocity estimates when localization is dependent on the MB position within the beamformed pixel of size dx. A MB of constant steady-state velocity u GTss (r) passes a beamformed pixel every dx/u GTss (r), i.e. with frequency f dx = u GTss (r)/dx. Trajectory velocity estimates should inhibit this same frequency in case of grid-dependent MB localization, since they are computed directly from MB positions. Here, the simulated MBs move along with the pulsatile flow u GT n (r). Therefore f dx corresponds to a range of frequencies as We inspected the frequency spectrum of u n (r) to assess the presence of these frequency-dependent grid-based artifacts.
To filter out the artifacts which could not be eliminated by using the moving average filter, we applied a bandstop filter in the frequency domain as a pre-filtering step. We set the stopband corresponding to the derived f dx range. Then we applied a moving average filter with span of 51 frames to improve both the retrieval of the velocity extremes as well as the visual appearance of the trajectories. We defined the filtered velocity estimates obtained from our designed filtering process as u f n (r). We implemented this pulsatility retrieval method on trajectories measured in the double vessel simulation for the Localization scenario. The pulsatility fraction was calculated on the filtered single-trajectory velocity estimates. To retrieve the maximum and minimum velocity measured during a trajectory the average of three points centered around respectively the local maximum and minimum was taken.  Table 2

2) ACCUMULATION OF TRAJECTORY VELOCITIES AS OBSERVED BY A STATIONARY OBSERVER
We introduce a second method of retrieving pulsatility that relies on the accumulation of velocity estimates in a fixed lateral position r in the vessel, following an Eulerian flow description. As multiple MBs go through the same super-resolution pixel during a ULM acquisition, multiple velocity estimates are found in each pixel forming a velocity estimate distribution. The average of these estimates is reported in a conventional ULM velocity map. This method retrieves a P f estimate from the shape of the velocity estimate distribution.
To design this method we start by considering the expected shape of the velocity distribution at a fixed r for a fully one-dimensional flow along the axial coordinate l. Note that a one-dimensional flow is considered that is defined in a two-dimensional space (l, r). Since the flow is assumed to be fully developed, it is independent of l. We hypothesize a steady-state flow, i.e. u GT n (r) = u GTss (r), after which the case of a pulsatile flow is considered, see Fig. 2

(a).
For steady-state flow, we assume that the spread in measured velocities is solely caused by error in MB localization. Consider two linked MBs at frame n and n − 1 as illustrated in Fig. 2 (b). Each estimated MB positionl n along axial coordinate l is given by the ground truth MB position l n with an added localization error l n , i.e.l n = l n + l n . The velocity estimate at frame n is then given as where u n is the velocity error caused by error in MB localization. We assume l n and l n−1 to be independent and N (0, σ 2 ) distributed, with σ the localization precision. The distribution of the velocity error u n can be found by a combination of the two independent normal distributions, resulting in U ∼ N 0, 2 dt 2 σ 2 . The velocity estimates of a trajectory are smoothed by a moving average filter. Substitution of (5) in (2) for general odd-numbered s results in the filtered velocity estimate as where u GTf n (r) is the filtered ground truth flow and u f n the filtered velocity error resulting from localization error. The distribution of the filtered velocity error is found to be In the steady-state scenario, where the first term of (7) can be replaced by u GTss (r), we find the filtered velocity estimates at a fixed lateral location r to be distributed as U f ∼ N (u GTss (r), σ 2 u ). In a pulsatile flow, the filtered ground truth flow u GTf n (r) does not equal the steady state flow u GTss (r). The distribution of u f n (r) for pulsatile flow is then found to be U f ∼ N u GTf n (r), σ 2 v . Two effects cause the spread in this distribution. Similar to the steady-state case, the variance σ u is caused by the localization error. In contrast to the steady-state case, an additional spread due to a varying mean u GTf n (r) is caused by the pulsatile flow behaviour. The distribution U GTf of the varying mean can be approximated by creating a histogram of u GTf n (r) over one pulsatility cycle, as in Fig. 2 (c). Assuming u GTf n (r) to be independent of u f n we find the distribution of u f n (r) in pulsatile flow to be where represents the convolution operation. In Fig. 2 (c) the derivation is performed for three different values of P f for a fixed σ and s. By employing the moving average filter, it is implicitly assumed that u GTf n (r) ≈ u GT n (r). The span s VOLUME 2, 2022 should be small enough compared to the pulsatile period for this assumption to hold. The minimum and maximum values of u GTf n (r) are given by the peaks in the histogram U GTf . After convolution, these peaks appear in U f as well. In ULM acquisitions U f can be found by the accumulation of the filtered velocities found at a fixed lateral location r. The pulsatility fraction can be estimated by substitution of the found peaks in (1).

III. RESULTS
The results presented here are divided into three sections. We first report the performance of the simulator, then the simulation results of the Lagrangian pulsatility retrieval method that relies on the filtering of single MB trajectories, and finally theoretical and simulation results found with the Eulerian method.

A. THE FULL ULM PROCESS IS SIMULATED AND SUPER-RESOLUTION IS ACHIEVED
An average velocity map was reconstructed from the simulation of the flow in the double vessel configuration in Fig. 3. This average velocity map serves as the ground truth with which the ULM velocity reconstruction should be compared. Fig. 3 (b) displays a simulated B-mode image containing three MBs.
Based on the acquired B-mode images, a contrastenhanced Power Doppler image is computed (Fig. 3 (c)), and the ULM pipeline for localization and tracking is implemented resulting in three super-resolved maps displaying density (Fig. 3 (d)), velocity (Fig. 3 (e)) and orientation ( Fig. 3 (f)). The cross-sections of the ULM density map and Power Doppler image taken at 3 different positions along the axial direction are plotted in Fig. 3(g). The two vessels can be separated as close as 15 µm (∼ λ/6), well below the wavelength of λ = 86.2 µm. Separation of the vessels is also visible in the velocity based renderings (Fig. 3 (h)). Finally, the orientation map displays a clear separation of the two vessels of opposite flow using a color code for the main orientation of the velocity vector, the venule appearing in blue and the arteriole in red. The root mean square error (RMSE) of the average blood flow velocity and orientation are 1.3 mm/s and 0.23 rad respectively.
Histograms of the found localization errors in x and z are plotted in Fig. 3(i). The localization precisions σ x and σ z were 8.5 µm and 12.7 µm respectively and are within the same range as reported in [31]

B. PULSATILITY FRACTION CAN BE RECOVERED FROM SINGLE TRAJECTORIES BY FILTERING OUT GRID-BASED ARTIFACTS
Three lateral vessel locations r 1 , r 2 and r 3 were defined, corresponding to the center of the venule, the side of the venule, and the center of the arteriole in the double vessel configuration. Trajectories at these lateral locations were found by controlling the inlet position of simulated MBs, see Fig. 4 (a) and (b). ULM parameters (see Table 3) were adjusted to recover longest trajectories.
The presence of a periodic grid-based artifact was noticed upon post-processing via a moving average filter with a span s = 51. At the position r 2 , the grid-based artifact is clearly visible by the sawtooth-like behavior (see Fig. 4(c)). At r 1 and r 3 this artifact was also noticed. Increasing the span s did not result in elimination of this artifact which can be seen in Fig. 7 of Appendix B. Therefore, additional filtering is needed to eliminate the artifact.
The frequency spectrum of the trajectory at r 2 was computed with a Fast Fourier Transform (Fig. 4 (d)). A frequency peak is visible in the frequency range f dx at which a simulated MB passes a beamformed pixel, computed from (3) and indicated as the gray shaded region. This peak was also found for different beamformed pixel sizes and is always in the f dx frequency range (Fig. 8 of Appendix B). This artifact is referred to as a beamforming grid-based artifact in the rest of this manuscript.
The filtered velocity estimates found by pre-filtering with the bandstop filter and smoothing with the moving average filter are shown in Fig. 4 (e). The pulsatility fraction estimates calculated from the found maximum and minimum velocities are given in Fig. 4 (f). The RMSE of the pulsatility fraction estimates at r 1 , r 2 and r 3 are 0.022, 0.056 and 0.013, which corresponds to 11%, 28% and 3.3% of the simulated pulsatility fraction respectively.

C. THE ABILITY TO RETRIEVE PULSATILITY FROM THE MEASURED VELOCITY DISTRIBUTION IS DICTATED BY σ AND s 1) THEORETICAL RESULTS
The theoretical velocity distributions corresponding to the single vessel configuration are given in Fig. 5 (a). The calculation was performed for different values of σ and a fixed moving average span of s = 21. In Fig. 5 (b) and 5 (c) the calculation was repeated for pulsatility fraction of P f = 0.2 and P f = 0.
Pulsatility introduces a bimodal velocity distribution for sufficiently low σ . For P f = 0.4 two peaks could be retrieved for σ ≤ 20 µm, while for P f = 0.2 this was only possible for σ ≤ 10 µm. The maximum and minimum velocity during the pulsatile cycle were retrieved from the location of these peaks and a pulsatility fraction estimateP f was found by their application in (1), see Fig. 5 (a)-(c). No bimodality was found in the distribution of the steady-state flow case. The location of the centroid of the distribution (shown with asterisks in Fig. 5) indicates the meanū GTf n (r). When rendering the MB trajectories to a ULM reconstruction this mean value is obtained and displayed in the velocity map.ū GTf n (r) is the same for all distributions given in Fig. 5(a)-(d).
From the derivation in II-B.2 it was already found that increasing the span s narrows U f . This effect was also found in U f when comparing Fig. 5 (a) with (d) with respectively s = 21 and s = 51. We find that increasing the span aids the

(g) Cross-sections (1)-(3) of the ULM density rendering and Power Doppler image. (h) Cross-sections (4)-(6) of the ULM velocity rendering (i) Histograms for measured localization errors x and z for a total of 29417 localized MBs.
retrieval of the bimodality. In Fig. 5 (d) two peaks could still be retrieved for σ = 30 µm.
The localization precision σ influences the shape of the distributions. For increasing σ the peaks move inwards and the distribution flattens until no two peaks can be detected any longer at σ = 20 µm for s = 21 and σ > 30 µm for s = 51. To further study this, Fig. 5(e) plots the pulsatility fraction estimate as a function of σ for a fixed value of the VOLUME 2, 2022 span s = 21. The ground truth pulsatility fraction is indicated as a dashed line. We find that P f is biased and always underestimated. A deteriorated localization precision results in an increased underestimation of P f , well in accordance to what is observed in Fig. 5(a)-(d). For σ > 20 µm only a single peak was detected and thus, noP f was found.
Increasing the span of the moving average filter s leads to increasing of the height of the velocity distribution peaks. The velocity distributions are plotted in Fig. 5(a) and (d). The peaks also move inwards for higher s. To further study this, Fig. 5(f) plots the pulsatility fraction estimate as a function of the span s for a fixed value of σ up to s = 51. In accordance to what was previously described,P f is underestimated for all the scenarios. The smallest theoretical attainable estimation error is found at the maxima of the curves for each value of σ . With this approach an optimal setting of s can be derived for any obtained localization precision σ , as is shown in Fig. 5(f). The curves show that a lower localization precision (higher σ ) requires a larger span s to enable retrieval of the pulsatility fraction.

2) SIMULATION RESULTS
The velocity distributions as derived in Fig. 5(a) were acquired in simulation for σ = 5, 10 and 20 µm (single vessel No-Localization). The resulting histograms are given in Fig. 6. The location of the peaks of the distributions found in simulation are in good accordance with the theoretical peak location. The simulations for σ = 5 µm and for σ = 10 µm result in the same peak locations for the bin size of the histogram used here.
Additionally, the method was applied to the same data set as used in Fig. 4, where the double vessel configuration was simulated with the Localization scenario. Histograms were obtained by accumulation of velocity estimates at the lateral coordinates r 1 , r 2 and r 3 (see Fig. 6 (b)). In the center of the venule (r 1 ), we foundP f = 0.143 which is an underestimation of 28.6%. It was needed to increase the moving average span to 31 frames for retrieval of bimodality. In the histogram found for r 2 no bimodality could be retrieved despite further increasing the span. In the center of the arteriole (r 3 ) a pulsatility fraction estimate of 0.379 corresponding to 5.2% underestimation was found. Bimodality was already achieved at s = 21. Using the Lagrangian method we achieved a RMSE of 3% and 11% for respectively the arteriole and the venule. We find that pulsatility retrieval is best performed at the centerline of the vessel.

IV. DISCUSSION
In this study, we introduced two methods to retrieve flow pulsatility with ULM. The Lagrangian method retrieves pulsatility from single-trajectory velocity estimates after a twostep filter. The second method takes on a Eulerian approach of flow modeling, and relies on the accumulation of velocity estimates at a fixed location in the microvessel of interest. With this approach, the pulsatility fraction could be retrieved from the location of the two peaks in the bimodal velocity distribution. To validate these methods, we modelled pulsatility in arterioles and venules from experimental data and used the Verasonics Research Ultrasound Simulator to generate ultrasound B-mode frames. The low computation time 290 VOLUME 2, 2022 of this simulator (as reported in Appendix A) compared to k-space-based simulators makes it suitable for large dataset generation. Our results show that ULM datasets contain rich temporal information that can be processed to retrieve pulsatility in addition to conventional ULM image reconstructions. By relying on kilohertz frame rates, we were able to retrieve the temporal dynamics induced by pulsatility without changing the localization, concentration, or tracking of the ULM processing pipeline. In application to our simulated data set, the Lagrangian method outperforms the Eulerian method based on the reported errors in P f estimate. Potentially, the estimates of both methods could be combined to create a P f map of the vasculature.
Our derived velocity distributions were validated in simulation. The peak locations in the acquired histograms (Fig. 6) closely match those derived theoretically (Fig. 5). The larger presence of high velocities in the histograms acquired from the simulation is hypothesized to be a consequence of the spatial interpolation of the trajectories. A MB of higher velocity travels a longer distance on a frame-to-frame basis, causing its velocity to be accumulated in the velocity distribution of multiple pixels. This effect is beneficial to the retrieval of pulsatility since it causes the peak corresponding to the maximum velocity to be more prominent. The derivation of U f should be updated to include this effect. Pulsatility induces bimodality in the distribution of filtered velocity estimates U f (Fig. 5). Localization precision σ dictated both the ability to find a P f estimate as well as the quality of that estimate. In the derivation of U f the MB localization error was assumed to be normally distributed as ∼ N (0, σ 2 ), which does not always hold for all localization algorithms [31].
In the Lagrangian method, both the pre-filtering by the bandstop filter and the smoothing by the moving average filter were found critical to retrieve pulsatility from single trajectories. In the trajectory velocity estimates, we observed a large presence of noise resulting from MB localization error. In the simulation of the double vessel configuration, we obtained σ values higher than the average distance traveled by the MBs over one frame, 12.7 µm versus 9.0, 3.3 and 8.6 µm for r 1 , r 2 and r 3 respectively. Additionally, we discovered the presence of beamforming grid-based artifacts in single-trajectory velocity estimates caused by grid-dependent localization, as shown in Fig. 4(c). The presence of these artifacts is consistent for different beamformed pixel sizes dx and for localization using both a radial symmetry and a Gaussian fitting algorithm (Fig. 8). Applying the Lagrangian method to a dataset without speckle for which no SVD filter was applied, validated that the artifacts result solely from beamforming (Fig. 9). To eliminate the grid-based artifacts, we designed a bandstop frequency filter with a stopband f dx based on available ground truth MB velocity. In practical applications this ground truth information is not available. In this case, an iterative process based on the measured MB velocities and its frequency spectrum should be applied to determine an appropriate f dx . Additionally, further research could focus on the replacement of the moving average filter with an estimator that utilizes all intermediate MB localizations within the assigned span. Applying the minimum The Lagrangian method relies on the temporal sampling of a trajectory. For each trajectory, it will recover a pulsatility fraction using (1). In the case of slow MBs -at the side of the venule for example -trajectories spanning hundreds of consecutive frames could be found (Fig. 4(e)). The main factor limiting the trajectory length in our simulation was the time during which the MB was present in the vessel. Acquisition of long MB trajectories is therefore crucial in the Lagrangian method. The trajectories in Fig. 4(e) succeed each other every ∼100 frames. Due to the sparsity of in vivo MB data the intervals between two succeeding trajectories are expected to be longer. The reported RMSE for the different lateral locations scale inversely with the absolute velocity deviation during the pulsatile cycle. Therefore, pulsatility retrieval can best be performed at the center of a vessel.
The Eulerian method explores the distribution of the velocity estimates at a fixed spatial sample, rather than inspecting the velocity estimates from individual MB trajectories. Because it relies on the full acquisition time, we expect less dependency on the frame-rate and on the retrieval of long MB trajectories than the Lagrangian method. In the calculation of U f we noticed two possible effects a change in frame rate might cause. First, the temporal discretization causes an additional numerical error in (6). This numerical error will increase for lower frame rates with O(dt 2 ). Second, the time step dt, experimentally defined by the frame-rate, influences the variance of the filtered velocity error through σ u = √ 2 sdt . When the span of the moving average filter is increased to match the original temporal window, this effect is compensated for.
The moving average filter was found to narrow the velocity error distribution U f by dividing the standard deviation by the value of s. It therefore reduces the adverse effect of localization error on the velocity estimates. For estimation of a steady-state flow it is advised to use a moving average filter that spans the full trajectory length. However, beyond a certain point increasing the span has adverse effects on pulsatility retrieval (see Fig. 5(f)). The heavy smoothing leads to an increased underestimation of the pulsatility fraction due to the flattening of the pulsatile cycle (see Appendix C). This underestimation effect of the moving average filter is apparent in Fig. 5(e) where s = 21 was applied. Even for perfect localization (σ = 0 µm) pulsatility fraction was underestimated. Unfortunately, the span is currently not reported in ULM studies. The theory introduced here can be used to determine an optimal setting of the moving average span for a specific flow scenario with given localization precision.
The occurrence of false MB pairing by the tracking algorithm will influence both methods. It will cause unexpected behavior of the filtered trajectories (Fig. 4(e)) and introduce additional contributions to the measured velocity histogram (Fig. 6). To limit the adverse effect of false MB pairings, an appropriate tuning of the maximum linking distance in ULM tracking is necessary. The theoretical derivation of 292 VOLUME 2, 2022 the velocity distribution performed in this study aids the development of a systematic method to determine the optimal maximum linking distance (Appendix D).
Pulsatility was simulated by our ULM simulator by modelling the spatial and temporal characteristics of the microvascular blood flow separately corresponding to a straight vessel with rigid wall of constant diameter. The ULM simulator could be improved by adapting a more complex vascular architecture corresponding to in vivo observations. Additionally, introducing a third dimension in acoustic simulations would describe more faithfully the physics at play in this problem. The derivations of the velocity distribution should be extended to take the elevational projection into account.
By investigating the effect of MB localization precision on the retrieval of pulsatility, we focused on the spatial quality of MB trajectories. Future research could focus on the temporal sampling of the MB trajectories in relation to pulsatility reconstruction. We expect the minimum frame rate needed for resolving hemodynamics to depend mainly on the time scale of the targeted hemodynamic phenomena, which are usually in the tens of milliseconds to seconds timescale for hemodynamic events. The methods reported here will need to be validated on experimental data. We anticipate specific challenges in an in vivo context, such as shorter MB trajectories, sparser MB data and incomplete filling of vessels for short-time acquisitions. In-vivo it is expected that the vascular architecture contains crossings that will lead to false MB linking. Since the Eulerian method does not rely on capturing long MB trajectories, we expect it to be more reliable for use on experimental data. However, for higher sampling rate of the pulsatility, the Lagrangian method will be more appropriate as each trajectory can give a pulsatility measurement and thus will lower acquisition time.

V. CONCLUSION
We have introduced two methods to extract pulsatility from raw ULM datasets. First, by filtering out grid-based artifacts from single trajectory velocity estimates, second, by retrieving the distribution of velocities found at a fixed location in a vessel of interest. Both methods have been validated on simulated ultrasound beamformed data featuring pulsatility induced flow variations. Our study shows that ULM datasets contain more information on the hemodynamics of the blood flow than that provided by conventional ULM image reconstructions. By looking at reconstructed trajectories independently or at the distribution of velocities, more insight can be gained into the hemodynamics experienced by the microbubbles. By reporting vascular anatomy and hemodynamic function, ULM has the potential of becoming a full-fledged ultrasound diagnostic method of unprecedented resolution.

DATA AND CODE AVAILABILITY
The simulated data is available through the 4TU.Resear-chData portal at https://doi.org/10.4121/21517878. The code  that incorporates the two pulsatility retrieval methods is available at https://github.com/qnano/ulm-pulsatility.

APPENDIX A
All symbols used in this study are given in Table 2 along with their definition.
Total computation time of the ULM simulator was on average 24 ms/B-mode frame excluding the time needed for computing the speckle frames. Therefore, simulating a ULM reconstruction of 10s acquisition time at 1kHz requires VOLUME 2, 2022 FIGURE 7. The periodical grid-based artifact found in single-trajectory velocity estimates is not eliminated by increasing the span s (in number of frames). The velocity estimates filtered with a moving average filter of s = 111 still show the artifact while their shape is heavily altered by the heavy smoothing, resembling the effect shown in Fig. 10(a). FIGURE 8. The frequency spectra of the single-trajectory velocity estimates u n (r) are shown for simulations in which beamforming was perfomed on different pixel sizes dx. f dx corresponds to the range of frequencies at which a MB passes a pixel. In the top row with dx = 100 µm, peaks lie within the frequency range f dx for localization performed with a radial symmetry algorithm and with Gaussian fitting. In the bottom row a peak is visible in the f dx range for dx = 50 µm. For dx = 10 µm ∼ 1/10λ the frequency range was not found to align with a distinct peak in the spectrum. On the right an illustration of the grid-based artifact is included. The distance between the estimated MB positions (x,ẑ) is dependent on the ground truth MB location within (x, z) a beamformed pixel. a total average computation time of 240s. The hardware specifications of the computer used in this study are given in Table 5.

APPENDIX B
Increasing the moving average span did not eliminate the beamforming grid-based artifacts (Fig. 7). Even at a span of 111 frames, which is more than half of the duration of the pulsatile cycle, the artifact is still observed.
The artifact is likely to be a result of MB localization that is dependent on the location in the beamformed pixel. In Fig. 8 this is illustrated for a MB that moves at constant speed with localization biased towards the center of the pixel. The velocity estimates based on the estimated MB  positions (red) will fluctuate between a high value measured for the MB crossing a pixel border and a low value measured for a MB at a central pixel location (see the gray arrows in the illustration). To conclude that the artifact indeed results from this grid dependent localization, the simulation was performed on three different beamformed pixel sizes dx: 100 µm, 50 µm and 10 µm. The resulting spectra are shown in Fig. 8 with the f dx range given by the shaded area.
For dx = 100 µm and 50 µm a peak is found within the f dx range. In the case of a 10 µm pixel size, this is not observed. Note that this pixel size is approximately ∼ 1/10λ and typical localization precisions that would be obtained are often larger than 10 µm [31].
To validate that the grid-based artifact results solely from beamforming, a simulation on PSF frames without any speckle was performed. In processing of this data, no SVD filter was applied. The results in Fig. 9 show that the grid-based artifact is also present in this data set.

APPENDIX C
When increasing the span of the moving average filter, two effects cause the resulting velocity distribution U f to change.
Firstly, as introduced in section II-B.2, by employing a moving average filter one implicitly assumes that u GTf n (r) ≈ u GT n (r). This assumption is violated when increasing the span s, as can be seen in Fig. 10(a). Due to the flattening caused by the moving average filter, the filtered pulsatile cycle does not resemble the original pulsatile cylce (s = 1) any longer. As a result, the shape of the histograms of U GTf changes and their peaks move inwards.
Secondly, a more beneficial effect of increasing the span is the fact that U f narrows. After convolution with U GTf this gives rise to narrower peaks in U f . For lower localization precision a higher span is needed to be able to resolve two peaks. For example, in Fig. 10(d) for a localization precision of σ = 20 µm, the two peaks are more distinct for s = 51 than s = 21.

APPENDIX D
The ULM simulator was also used to inspect the influence of the pulsatility fraction on the rendered ULM velocity map. The quality of the reconstruction was assessed by the RMSE of the average velocity prediction at each super resolved pixel compared to the ground truth average velocity ( Fig. 11(a)). Additionally, the ULM velocity reconstructions were inspected along the cross-section (Fig. 12).
In Fig. 11(a) the velocity RMSE and the fraction of tracked MBs is plotted against the localization precision for different values of P f ranging from 0 to 1 and two maximum linking distance (MLD) settings. In tracking, the linking of two MBs is not permitted if the distance between them exceeds the MLD.
Generally the RMSE increases with σ . The reported RMSE values become unreliable, as seen by the enlarged 95% confidence interval, when only a low number of MBs is tracked. The fraction of MBs that is tracked drops for increasing values of σ , since the distance between two MBs becomes more likely to exceed the MLD.
The violins in the cross-sections of Fig. 12 report the spread inū f n (r) reported in the ULM velocity rendering at the 296 VOLUME 2, 2022  Fig. 11 and 12 show that P f has a small effect on the ULM velocity reconstructions through two mechanisms. Firstly, for an insufficient maximum linking distance the truncation of the velocity distributions of different P f results in different values reported in the ULM reconstruction, see MLD = 25 µm of Fig. 11(a) and 12(b). We observe that the fraction of MBs tracked varies for different P f . Fig. 11(b) displays the non-filtered velocity distributions U as found from U GT U for P f = 0, 0.2 and 0.4 and a fixed σ = 5 µm. For a MLD of 25 µm velocities over 25 mm/s can not be measured and the velocity distributions are truncated at that location. The truncation of these distributions by the insufficient setting of the MLD results in different locations of the new centroids, which represent the average velocityū f n (r) as found in a ULM reconstruction. Note that the non-filtered velocity distribution is applicable here since the MLD is applied before the moving average filter. Due to the truncation, the centerline velocity is found to be more underestimated for larger P f . The theoretical explanation of Fig. 11 corresponds to the simulation results of Fig. 12, where the underestimation at the center of the vessel is clearly visible.
Additionally, the effect of insufficient MLD is visible in the average density reconstruction over the cross-section. For MLD = 25 µm the number of MBs tracked drops for the central pixels.
A second effect of P f on the reported RMSE values is found in the high localization precision range (σ < 10 µm) even for sufficient MLD setting of 50 µm. This can not be caused by the truncation of the velocity distribution as described above since no loss in tracked number of MBs is found. When inspecting the cross-sections of Fig. 11 we only observe a slight increase in violin size for larger P f , which indicates a higher spread of reconstructed velocities in the ULM velocity rendering. We expect this effect to fade for longer acquisitions due to the averaging that is performed in rendering of the reconstruction. To construct these cross-sections, the single vessel configuration was simulated for a 10 s acquisition, while typical in vivo acquisition length is > 100 s.
To prevent underestimation of the centerline velocities, the MLD should ideally be set such that it captures the full non filtered velocity distribution. It should therefore be set significantly higher than the actual maximum blood flow velocity present. Based on our velocity distribution theory, new systematic ways of setting the MLD can be formed replacing the trial and error tuning that is currently common.