An Anatomically Realistic Simulation Framework for 3D Ultrasound Localization Microscopy

The resolution of 3D Ultrasound Localization Microscopy (ULM) is determined by acquisition parameters such as frequency and transducer geometry but also by microbubble (MB) concentration, which is linked to the total acquisition time needed to sample the vascular tree at different scales. In this study, we introduce a novel 3D anatomically-realistic ULM simulation framework based on two-photon microscopy (2PM) and in-vivo MB perfusion dynamics. As a proof of concept, using metrics such as MB localization error, MB count and network filling, we quantify the effect of MB concentration and PSF volume by varying probe transmit frequency (3-15 MHz). We found that while low frequencies can achieve sub-wavelength resolution as predicted by theory, they are also associated with prolonged acquisition times to map smaller vessels, thus limiting effective resolution (i.e., the smallest vessel that can be reconstructed). A linear relationship was found between the maximal MB concentration and the inverse of the point spread function (PSF) volume. Since inverse PSF volume roughly scales cubically with frequency, the reconstruction of the equivalent of 10 minutes at 15 MHz would require hours at 3 MHz. We expect that these findings can be leveraged to achieve effective reconstruction and serve as a guide for choosing optimal MB concentrations in ULM.

F OR long, we were bound by the limit of diffraction stated by Rayleigh's criterion [1] in conventional ultrasound imaging. Diffraction causes the point-spread function (PSF) to have a specific size in the order of the wavelength, which increases with required penetration depth: larger depths can be achieved at the cost of larger wavelength and thus degraded resolution. Wavelength thus dictates the compromise between limiting factors of resolution and penetration depth in conventional ultrasound [2]. The concept of isolated source localization inspired from optical imaging [3] and used in ultrasound localization microscopy (ULM) [2], [4], [5] challenges this compromise.
Microbubbles (MB) injected into the bloodstream are used as small, highly echogenic unique scatterers that can be tracked through vascular networks in large vessels down to, VOLUME 3, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ in principle, the smallest capillaries. It is by finding the exact centroid of each MB in reconstructed images or by fitting a paraboloid in radiofrequency (RF) data that we can obtain micrometer-scale resolution. ULM resolution is less dependent on the wavelength but more on the ability of the system to locate isolated microbubbles in ultrasound images [6].
The gain in spatial resolution comes at the cost of temporal resolution since sampling the vascular network requires localization and tracking of millions of MB in tens of thousands of images acquired during a few minutes [7]. This is particularly important for dynamic ULM applications such as pulsatility imaging [8], cardiac imaging [9], and functional imaging [10].
In [2], Couture et al. defined four conditions that are required to perform ULM, namely 1) being sensitive to the contrast agent, 2) being able to localize that contrast agent, 3) record at a sufficient framerate and 4) maintain isolated sources in a spatiotemporal referential. Our ability to respect the latter condition depends mainly on two factors: 1) acquisition parameters such as frequency and transducer geometry [6], and 2) microbubble concentration [7].
The first ULM studies demonstrating flow separability under the diffraction limit have been performed in in-vitro phantoms with various materials, configurations, and scale. Viessmann et al. [11] conducted four experiments at 2-2.5 MHz using two touching cellulose tubes with inner diameters (ID) in the order of λ/3, where λ is the wavelength. Two sets of MB could be distinguished with a sub-diffraction limit with a measured center-to-center distance of about λ/3. Similarly, in [12], O'Reilly and Hynynen conducted an experiment at 612 kHz using a PTFE tube with 255 µm of ID spiraled around a 2.5 mm rod. The estimated localization uncertainty was approximately 1/120λ in the x-y dimensions and 1/60 λ in the z dimension. In [13], Couture et al. localized MBs in 100 µm sized microfluidic channel at 5 MHz with an uncertainty of λ/250 axially and λ/38 laterally. Similarly, in [14], Desailly et al. succeeded in localizing MB in a 100 µm-sized microfluidic channel separated by 50 to 200 µm at 1.75 MHz. However, none of the studies could assess the limits of ULM vascular reconstruction.
In Heiles et al. [15], the performance of 7 MB localization algorithms was assessed in simulation using the Verasonics Research Ultrasound Simulator. The simulation media consisted of 1) a single MB moving in λ/21 increments, and 2) 11 simulated tubes disposed in various configurations and degrees of curvature while the smallest distance in flow being about λ/2. White Gaussian noise was added to match typical in-vivo data. Localization error was quantified as 0.13 λ, 0.36 λ, and 0.46 λ depending on the localization algorithm. In McCall et al. [16], localization error in the presence of image degradation was estimated in simulation using the Fullwave simulation tool. The simulation framework included inhomogeneous tissue modelling to account for reverberation clutter, trailing clutter, and phase aberration. Localization error was found to be 39 times worse in the axial dimension and up to 2.1 times worse in the lateral dimension than what is predicted by theory [17]. However, in the above simulation media, the MB distributions did not reflect in-vivo challenges.
In mice, brain vascular networks typically exhibit larger vessels with diameters in the 50-200 µm range [18] down to very small vessels in the <10-micrometer range [19]. Vascular density varies from 161-391 sections per mm 2 , or an equivalent of 51-79 µm between vessels on average [20]. Furthermore, the size of the vessel intrinsically dictates blood flow, i.e., the number of detectable MB events and their corresponding velocities, which affects effective resolution [7]. In Hingot et al. [21], in-vivo image resolution was quantified using Fourier Ring Correlation. ULM performed on rodent brain, kidney and tumor yielded an image resolution in the order of λ/10. However, to our knowledge, there is no threedimensional validation framework for ULM image formation algorithms based on anatomically-realistic MB flow to assess ULM limitations.
To overcome this problem, a highly resolved ground truth of vasculature is needed. Pulsed laser systems [22], [23], e.g., two-photon microscopy (2PM), can image in-vivo and ex-vivo tissue at micrometer scale. This approach can be scaled by stacking ex-vivo slices to reconstruct whole brain anatomy with micrometer scale resolution [24], [25]. Cerebral vascular networks that can therefore be leveraged as virtual phantoms [26]. A similar approach was developed in 2D in [27], but the missing information in the third dimension doesn't allow a complete characterization of MB localization error, network filling and other spatiotemporal metrics. We developed a 3D particle flow simulator based on highly resolved 2PM data from mice brain vasculature to conduct 3D anatomically-realistic simulations. It was first implemented in 2D as part of a deep learning framework [28] and as part of a sparse channel sampling study for ULM [29].
In this work, we introduce a novel 3D simulation framework for the validation of ULM image formation algorithms and provide an example application. From a vascular model obtained using 2PM, we generate a MB flow with known MB positions in time and where dynamics mimic MB behavior from in-vivo data [30]. These MB positions at controlled concentrations are used to simulate ultrasound raw data, which can then be input in a ULM image formation algorithm. This setup allows to evaluate the performance of MB localization frameworks with various imaging parameters, by comparing the output MB positions with the provided reference MBs. As an example application, we studied the effect of MB concentration and transmit frequency by computing metrics such as localization error and network filling. We provide new investigations that link maximal MB concentration with the inverse of the PSF volume. Such investigations are essential in optimizing MB concentration to obtain a more complete and accurate vascular reconstruction.

II. METHODS
In this section, we first introduce a simulation pipeline to obtain a virtual flow of reference MBs in a 3D vascular network. We then describe the generation of RF signals and a conventional ULM framework. Finally, we describe how we can leverage the reference MB flow to gain insight on vascular reconstruction.

A. IN-VITRO TWO-PHOTON MICROSCOPY
The pipeline starts with a highly resolved mouse brain vasculature acquired using 2PM from [32] (Fig. 1 A). The strain (C57BL/6) used is the most widely used in biomedical research [33]. Prior to imaging, C57BL/6 mice (25-30 g males, n = 6) were anesthetized using isoflurane (1-2% in a mixture of O 2 and air) and maintained at a temperature of 37 • C. The imaging apparatus consisted of a custom-built two-photon microscope with components listed in [34]. Imaging was performed through a cranial window with removed dura and a 150 µm-thick microscope coverslip for sealing purposes. Anesthesia was reduced to 0.7-1.2% isoflurane during the acquisition. Structural images of cortical vasculature were obtained by labelling blood plasma with dextran-conjugated fluorescein at 500 nM. Using angiograms acquired with a 20X Olympus objective (Numerical Aperture = 0.95), 6 stacks of vasculature were prepared with voxel sizes of 1.2 × 1.2×2.0 µm and total volume of 600 × 600× 662 µm.

B. GRAPH MODEL
Using a modeling framework from [31] and [35], a graph model of 2PM data was generated ( Fig. 1 B). In this framework, vascular segmentation was performed using a fully-convolutional neural network based on densely connected layers [31]. The graphing/skeletonization step was achieved by generating a 3D grid-graph model followed by a geometric graph contraction and refinement algorithms. Training the segmentation model was done using the Theano framework [36]. As for 3D modeling and geometric contraction processes [35], they are based on the VascGraph Python package [35]. Using this framework, a file containing vessel network geometry and topology (nodes and edges), with their corresponding inlet and outlet vessels is generated. The corresponding vessel radii are assigned as features to the nodes of the vascular network. Vessels ranged from 2 to 57.2 micrometers in diameter, their spatial position span was 500 × 500×650 µm and vessel volumetric density (VD) was calculated at 4%, using the ratio of vesselscontaining voxels and total voxels in segmented volumes.

C. PARTICLE TRAJECTORY GENERATION
The graph model that we obtain from the vascular data has nodes that are arbitrarily spaced and are located at the centerline of each vessel. However, we need MBs to be flowing at specific velocities and at different radial positions in the vascular network. We also need to measure the impact of MB concentration, therefore it must be kept constant throughout the simulation.
Before generating a steady state flow simulation where MB concentration is always kept constant, MB trajectories were computed using the graph model source nodes, target nodes and vessels radii with the shortestpathtree MATLAB (MathWorks) function.
For each trajectory, the path between a source node (i.e., has no parent node) and an extremity node (i.e., has no target node) was selected (Fig. 2). Since the graph model was a tree graph, only one trajectory was possible from a source node to an extremity node. A cubic spline fitting was then applied to allow for the generation of MB positions with a controlled time-dependent spacing. In Hingot et al. [7], MBs were followed for 5 minutes in a rat model with a continuous VOLUME 3, 2023 where parameters p 1 = 1.9 mm µm −1 s −1 and p 2 = −6 mm.s −1 . Another log-log fit (2) ln(c) = p 3 ln (d) + p 4 (2) was computed on MB count (c) as a function of diameter (d) with parameters p 3 = 3.7 MB µm −1 , p 4 = −8.2 MBs, which is used in the next section (II-D) for trajectory selection. We modelled the MB flow in a vessel with assumptions of a Poiseuille flow in a tube as an approximation of the normal physiological state [37]. Such flow exhibits a parabolic axial velocity profile denoted as: where v max is the axial velocity at the centerline of the tube, r is the radial distance from the centerline, R is the radius of the tube, and v is the axial velocity at position r. In our case, R corresponds to the vessel radius of the closest original node, r is the radial distance of the new MB position from the vessel centerline, and v max was calculated using its relationship with the vessel diameter from equation (1). A minimum r/R ratio was set at 0.3 to avoid quasi-static MBs (an r/R ratio equal to zero would require an infinite number of samples in a given trajectory). A MB trajectory is terminated if the trajectory reaches an end node or if the MB can no longer flow in the vessel because it is larger than the vessel. In this case, we modelled the MBs to be 2 micrometers in diameter. With the centerline trajectory, an infinite number of trajectories that run parallel to it can be generated. To do so, two orthogonal sets of vectors were calculated for each trajectory (Fig. 2, green and blue arrows). Using a linear combination, we can span all possible trajectories to fill the vessel. This step allowed for the generation of multiple trajectories at different radial positions in a specific vessel according to their r/R ratio. The radial boundary of each MB position was effectively the vessel diameter of the closest node. Lastly, the size of the MB was also considered where its radial span was the vessel diameter minus the MB diameter. For example, a MB of 2 µm in diameter travelling in a vessel of 2 µm in diameter [38] is limited to flowing only through the vessel centerline.

D. STEADY-STATE FLOW
Throughout this study, MB concentration is reported as the number of MBs per millimeter of beamformed volume (MB/mm 3 ). This allows for direct comparison between results shown here and one's ULM processing when visualizing data after MB extraction. Using computed MB trajectories, a steady-state flow was generated using a constant concentration of MBs ( Fig. 1 C). In practice, it meant achieving a constant number of MBs in the field of view at each time step. However, MB trajectories can't be randomly selected since the in-vivo distribution is dependent on the size of the vessel in which each MB is flowing in. To match in-vivo data from [7], MB trajectories were sorted according to their diameter. Then, MB count rate (c i ) was used to generate a probability density function by normalizing equation (3) as where X is the random variable representing the MB count, and where i is a trajectory index. We then used the pdfrnd MATLAB (MathWorks) function to select a random trajectory given the calculated probability obtained using equation (4). To populate the first set of MBs in the steady-state flow, a first population of particles was generated at random positions in their respective trajectory to match a desired concentration. Then, a new microbubble was added to the flow each time a microbubble reached the end of its trajectory. After this process, we were left with sets of constant number of microbubbles positions (X i , Y i , Z i ) at each time frame.

E. SCATTERERS' ARRANGEMENT
The dataset that we had access to spanned only about 0.5 × 0.5 × 0.65 mm in space, which corresponds to about 0.1625 mm 3 . It provides a realistically dense vascular network, but it does not span the entire field of view. Since the PSF varies based on its position in the field of view, we needed MB scatterers in most of the field of view to assess the performance of the system. We therefore decided to distribute unique MB trajectories of a vascular network throughout the field of view by translating their positions in a rigid manner in all 3 dimensions of space. MB trajectories were randomly distributed in a 10 × 10 × 10 pattern in the x, y, and z dimensions to span 5 × 5 × 6.5 mm (Fig. 3). This arrangement also allowed to span a usable range of MB concentrations. The original network spanned 0.1625 mm 3 , therefore if we take the minimal theoretical concentration of 2 MBs flowing in that volume span, it corresponds to about 12 MB/mm 3 , which is considerably high for lower transmit frequencies. Therefore, by redistributing MBs in space, we have in this case a 1000-fold discretization on the original MB concentration range. This way, MBs can be farther apart to test for very low concentrations. After the MB rearrangement, the minimal theoretical concentration of 2 MBs in the field of view can be as low as 0.012 MB/mm 3 since the new volume span is much larger. The upper limit on the MB concentration can be as high as needed simply by increasing the number of MBs in the simulation. After ultrasound simulation, reconstruction, localization, and metrics computation, MB positions from the 1000 sub-regions were added into one original small volume to emulate a longer acquisition time only for visualization purposes (Fig. 3). This step facilitated the qualitative comparison between density maps shown in Fig. 8.
In summary, to generate our ground truth MB positions in time, we performed the following steps: 1) acquire a 3D volume of brain tissue using 2PM 2) segment blood vessels from surrounding tissue 3) generate a labelled graph model of the vessels 4) calculate possible MB trajectories using diameterdependent MB flow rate and velocity 5) generate a constant set of MBs per time frame 6) rearrange groups of MBs to obtain a large field of view F. ULTRASOUND SIMULATION 3D probe parameters were set to match 32 × 32 matrix arrays 1024 EL 3 MHz [38], 1024 EL 8 MHz [39], [40] and the custom 1024-element 15 MHz as used by Brunner et al. [41]. The theoretical 12 MHz probe parameters were set to match the 3-and 8-MHz probes' specifications. Bandwidths and transmit frequencies are depicted in TABLE 1. Pitch and element width were set to 0.3 and 0.275 mm, respectively. The transmit pulse was a 0-degree plane wave, and the effective framerate was equal to the PRF of 1 kHz. The pulse length was chosen to be 4 cycles to enhance contrast. Since 2D imaging is widely used in ULM, we incorporated the study of a linear array. Probe parameters matched a L22-14 array with 128 elements, central frequency of 15.6 MHz, pitch of 0.1 mm and element width of 0.08 mm (Vermon, France).
To simulate the ultrasound images, an in-house GPU implementation of the SIMUS simulation software was performed as described in [42] and [43]. Microbubbles were modeled as linear point scatterers. White gaussian noise was added to the acoustic response to obtain a signal-to-noise ratio (SNR) of 15 dB in 2D and 10 dB in 3D to mimic typical SNR we find in our in-vivo data.

G. ULM IMAGE FORMATION ALGORITHM
In-phase-quadrature complex (IQ) volumes were generated from the RF data using a fast GPU-based delay-and-sum beamformer [44]. Individual MB were then identified as local maxima on correlation maps, resulting from the correlation of the reconstructed in-quadrature (IQ) volumes with the IQ spatial PSF of the imaging system. MB were located using a 3D subpixel gaussian fitting on IQ magnitude local maxima as in [8], [9], [28], and [45]. Localized MB were sorted as a function of their correlation from the correlation map. Only MB exceeding a certain correlation threshold were kept. Thresholding using correlation allows to screen interacting MBs. This is essential to meet the condition stated by [2], which requires isolated sources. To establish an optimal threshold, tests were conducted on various concentrations at specific transmit frequencies. After finding concentrations that yielded the most localized MBs, localization error was plotted against correlation threshold (Fig. 4). The threshold was selected to maximize MB count while preventing excessive localization error (Defined in Section II-I). In the example provided in Fig. 4, the selected threshold value was 0.6.
For visualization purposes and for reconstruction times calculations, density maps were generated. MB subpixel positions were accumulated over time on a 3D sampled volume, VOLUME 3, 2023  where each voxel is approximately one fiftieth of the wavelength in each dimension to generate a micro vasculature density map. On a density map, each voxel intensity represents the number of microbubbles in that specific voxel.

H. MB SIMULATOR METRICS
To validate MB behavior, 5000 MB trajectories were simulated to replicate similar MB count as in [7]. The two dependencies from in-vivo data are 1) the MB count as a function of vessel diameter and 2) MB velocity as a function of vessel diameter. To assess the first dependency, the occurrence of MBs flowing in specific vessel diameters was calculated. Natural logarithm (ln) was applied on both dependencies to generate a linear fit as in [7]. As for the second dependency, MB velocity was simply plotted against its corresponding vessel diameter.

I. ULM METRICS
To quantify the effect of concentration on the localization itself, localization error was computed using highly correlated MB positions. First, each localized MB was paired with the closest reference MB using the Euclidian distance.
For a specific concentration, localization error was calculated as the average distance between localized MBs and their corresponding closest reference MB. Only the in-plane distances were considered when computing 2D localization error of the linear array, i.e., the error in the elevation direction was ignored. This process was performed in 2D at 15.625 MHz and in 3D for different transmit frequencies spanning current usable range: 3, 3.5, 6, 7, 9 12 and 15 MHz.
Since the MB count is directly linked to acquisition time [7], it was used as a metric to quantify the latter. After correlation thresholding, the MB count was measured at the mentioned transmit frequencies and at various concentrations. Concentrations were chosen iteratively to span a range where the MB count naturally increases with an increase in concentration and decreases when MB interactions become more important.
To assess effective reconstruction, density maps were generated using accumulated MB positions over different acquisition times. We chose a transmit frequency of 6 MHz at four different concentrations: 0.16, 0.4, 0.8 and 1.6 MB/mm 3 . They correspond to 20, 50, 100 and 200 MB in the field of view and to 10%, 25%, 50% and 100% of the maximal MB concentration derived from Fig. 6. The four concentrations were chosen iteratively, using localization error and number of localized MBs as guidelines to span a usable concentration range and to depict differences in the resulting vasculature reconstruction. A reference density map with a fully populated network was generated by simulating all possible MB trajectories without a constraint on MB count (N) vs diameter (d). A quantitative comparison with the reference density map was made at different acquisition times using the Sørensen-Dice coefficient [46], [47]: a similarity index quantifying the network filling as in [28].
To establish the link between PSF size and MB concentration, the PSF volume was calculated for each of the different transmit frequencies and concentration at which MB count is maximal. PSF was modelled as an ellipsoid with semi-axes calculated using the full width at half maximum (FWHM) in each dimension. For 2D imaging, ellipsoid semi-axis in the elevation plane was calculated as half the focal elevation width.  Fig. 5 displays the two dependencies that were originally extracted from in-vivo data [7]. The black lines represent the reference dependencies, colored dots represent the statistics of simulated MBs, and the blue dashed lines represent the linear fits on the simulated data. Each dot color corresponds to one of six vascular networks. Results show that the MB count increases with vessel diameter with a slope of 4.15 MB µm −1 compared to the reference 3.7 MB µm −1 . MB velocity also increases with vessel diameter with a slope of 1.77 mm µm −1 s −1 compared to the reference 1.9 mm µm −1 s −1 . Fig. 7 shows the different MB flow simulations from six vascular networks. Each of the simulations contains an accumulated 1 million MB positions. In Fig. 7, the top section depicts simulations where MB flow was constrained to follow the dependencies mentioned in section II-C. The bottom section shows MB flow simulations without the constraint on MB count vs the vessel diameter. The latter simulates a random MB distribution.

B. LOCALIZATION ERROR
In Fig. 6, localization error is depicted as a function of microbubble concentration for different transmit frequencies. Fig. 6 also displays the count of localized MBs as a function of concentration. At lower frequencies, MB count increases until reaching a maximum and decreases at higher concentrations. However, from 9 MHz, MB count reaches a plateau at higher concentrations. We see that at 3 MHz, the number of detected MB is approximately 10 000 MB/s or 10 MBs per frame. If we double the transmit frequency (6 MHz), the MB count increases with an order magnitude with approximately 100 000 to 200 000 MB/s. The rate at which the number of MB events increases seems to decrease with higher frequencies. See Supplementary Fig. 1 for localization error and MB count as a function of MB concentration in units of MB/λ 3 .

FIGURE 6. (Top left) 2D localization error is depicted as a function of microbubble concentration for the 15.625 MHz linear array. (Bottom left) 3D localization error is depicted as a function of microbubble concentration for matrix arrays at various transmit frequencies. (Top right). MB count resulting from localization as a function of concentration is depicted for the 15.625 MHz linear array. (Bottom right) MB count resulting from localization as a function of concentration is depicted for different transmit frequencies using matrix arrays. Bold curves
represent the average and shaded areas cover the average ± standard deviation of the 6 vascular networks. Fig. 8 shows the effect of concentration on vascular reconstruction. Top row shows density maps acquired at different transmit frequencies, increasing from left to right. MB concentration was selected as 50% of the maximal concentration (0.12 MB/mm 3 ) for a 3-MHz transmit frequency. Bottom row shows the density maps of sub-regions of reference and localized microbubbles at increasing concentrations from left to right. Fig. 9 shows the link between maximal MB concentration and PSF volume, where maximal MB concentration is proportional to the inverse of PSF volume. From left to right in blue, inverse PSF volumes correspond to transmit frequencies of 3, 3.5, 6, 7, 9, 12 and 15 MHz of the matrix arrays. The dashed orange line represents a linear fit as y = 0.043x + 0.51, R 2 = 0.94. The inverse PSF volume x is in mm −3 and maximal MB concentration y is in MB/mm 3 . The uncertainty bars span the mean ± the standard deviation of the MB concentration of maximum MB count from the 6 vascular networks.   a framerate of 1000 frames per second, where each frame corresponds to an imaging volume. MB concentration in each simulation corresponds to 10% of the concentration where MB count is maximal for each transmit frequency. 10% was the concentration that yielded what was qualitatively the best vascular reconstruction.

G. RUNTIME PERFORMANCE
The computer used for this study was an Intel Core i7-7820X @ 3.6 GHz, 128 GB RAM, NVIDIA GeForce RTX 2080 Ti. The generation of 4.5 million MB trajectories sampled at 1000 Hz was performed in approximately 19 hours. Radio-frequency signals of 1000 volumes, 1024 transducer channels, at 1000 Hz were generated in 58 seconds. Beamforming RF signals at λ/2 into 1000 IQ volumes took 54 seconds. Generating correlation maps between the PSF and IQ data took 54 seconds with a kernel size of 9. Microbubble localization was performed in 51 seconds.

IV. DISCUSSION
In this study, we took advantage of a novel 3D anatomicallyrealistic framework, where mice brain vasculature obtained from 2PM was used to propagate virtual MBs with velocities and MB distribution replicating behavior from in-vivo data. Having access to known MB positions at all times allowed for quantitative measurements to determine ULM capabilities as well as reference qualitative imaging. 2D and 3D localization error could be computed as well as MB count for different MB concentrations and transmit frequencies. A link between maximal MB concentration and PSF volume was found in matrix array imaging and network filling as a function of time was calculated. Results allow us to answer fundamental questions related to in-vivo imaging such as whether the vasculature is accurately or fully depicted as imaging parameters and MB concentration are varied.
Simulations were conducted on the six vascular networks separately. The variety of vascular morphologies allowed to generalize the link between MB concentration and acquisition parameters. Fig. 5 shows that our MB simulation framework was successful at replicating observed MB behavior from in-vivo dependencies, namely the MB count and velocity vs vessel diameter. Due to these dependencies, realistic MB flow is possible in subsequent ultrasound simulations. Fig. 7 shows that constraints on MB flow favor MBs flowing in larger vessels. Visually, fewer small vessels are mapped for a set amount of MB events, but large vessels appear fuller. This means that for a set acquisition time and a set MB concentration, some small vessels will not or will not be fully mapped in the final density map regardless of imaging parameters.

B. MB LOCALIZATION
In Fig. 8, we see that as concentration increases, even though most vessels are preserved, they appear blurred, are missing, or are inaccurately depicted from the reference density map. This is in accordance with theory since ULM requires isolated MB signals. The higher the concentration, the more likely MBs will interact with each other. However, relative concentration around 10% of the concentration at maximal MB count allows to depict most of the original vasculature. Furthermore, reference density map shows that because of realistic MB count vs diameter, some small vessels are not present. Probability of a MB flowing in vessels of only a few micrometers is so low that it would require a longer acquisition time than the equivalent of 10 minutes at 15 MHz, according to results in Fig. 10. This means that in current in-vivo imaging, some vessels of only a few micrometers are most likely not mapped if the transmit frequency is not high enough or the acquisition time not long enough. The visual assessment of the vascular reconstruction (Fig. 8) can also show the behavior of the localization processing in challenging areas, such as vessel bifurcations, where MB interactions are more probable.
In Fig. 6, we can see that MB count increases as a function of concentration as similarly observed in [7]. On the other hand, the decrease in MB count at high MB concentrations shows that thresholding MBs using their correlation with beamformed PSF allows to screen for overlapping MB PSFs. Moreover, we can observe that the usable MB concentration range increases radically with higher transmit frequency. This would mean that the margin of error is much lower for low transmit frequencies than for high transmit frequencies.
Localization error increases with concentration as predicted by theory and is in conjunction with the blurring or the non-detection of vessels. Higher transmit frequencies are associated with lower localization errors as predicted by theory [48]. Higher frequencies are associated with smaller PSFs -hence a higher concentration is possible while maintaining isolated sources.
In Fig. 6, the rate at which maximum MB count increases eventually slows down with higher frequencies. This is mainly due to the size of the PSF only decreasing in the axial direction because of constant aperture. Moreover, since the pitch is maintained constant while increasing the transmit frequency, grating lobes can be introduced. For example, a pitch of 0.3 mm in matrix probes operated at 3, 6, and 15 MHz have a respective normalized pitch of about 0.6λ, 1.2λ, and 2.9λ. We know that the pitch influences the position of grating lobes, and that the element width dictates the angular response weighting [49], [50]. To avoid grating lobes at any steering angle, the pitch must be less or equal to 0.5λ. If the pitch is equal to λ, then the grating lobes are located at ± 90 degrees of steering angle. When the pitch is equal to 2λ, grating lobes are located at ± 30 degrees. Therefore, at higher frequencies, we can expect more contribution of grating lobes to the PSF. Thus, a large pitch can limit our ability to separate MBs at higher concentrations.
The dependency in Fig. 9 shows that PSF volume can be used to predict maximal MB concentration. In theory, to reduce MB interactions, MB concentration must be decreased. We see here that reducing PSF size also reduces MB interactions. Since MB concentration should be as high as possible to reduce acquisition time, this dependency can be used for choosing optimal concentration for a specific probe's acquisition parameters.
In Fig. 10, the network filling represents the intersection between reference density maps and density maps resulting from MB localization at different transmit frequencies. The higher the network filling index, the more reference vessels are present in the mapped vasculature from the imaging setup. Results in Fig. 10 show that acquisition time is drastically reduced using smaller PSF volumes. We already know that a higher transmit frequency should be prioritized VOLUME 3, 2023 over a lower transmit frequency if the required penetration depth allows for it. It is also important to keep in mind that for a set acquisition time, depending on MB concentration and PSF volume, only a portion of the vascular network is mapped.
As for the runtime performance, it took about one day to simulate all MB trajectories and a steady-state flow simulation. Since that step was only performed once, runtime was of less importance. The bottleneck in terms of runtime in this study was the simulation pipeline after the generation of a MB flow, i.e., simulating RF signals, beamforming, and the running the ULM framework to localize individual MBs. The ultrasound simulation framework was already optimized, but since we varied MB concentration as well as transmit frequency, we had to simulate about a hundred thousand volumes of ultrasound data for each of the concentration-frequency set of parameters.

C. LIMITATIONS
This study set forth some limitations of the proposed simulation framework. First, a graph model was generated from 2PM, through automated segmentation, surface modeling and contraction. From that, MB trajectories were generated using MB velocities and number of MB as a function of diameter from two models derived from in-vivo data. This approach allows for anatomically-realistic modelling since vasculature is derived from an ex-vivo model. Another avenue would consist of generating a fully synthetic fluid flow-driven graph model as in [51], where vasculature is generated using constraining parameters, namely vascular resistance, diameter, pressure, and bifurcation position. Whole cortical circulation can be generated while distinguishing between arterial and venal flow. That framework combines anatomical data with artificial construction laws to overcome limitations in coverage and resolution in reference anatomical data.
As for variability in the MB count vs vessel diameter in Fig. 5, it is due to the specificity of the vascular network. The smaller the vascular network in terms of field of view, the more significant this variability is. Variability in MB velocities vs vessel diameter can be due to numerical error between the coarse sampling of the original labeled nodes compared to the finely sampled new MB positions.
Moreover, the reference MB velocity and MB count versus the vessel diameter dependencies were sourced from a study with a 150-vessel sample in a rat model, as opposed to a mouse model. A recent study by Demeulenaere et al. [52] on a whole mouse brain showed MB velocity and vessel diameter measurements. They are in the same order of magnitude as in a rat model but appear to span higher velocities for an equivalent vessel diameter range of 10-100 µm. Therefore, using reference data from a mouse model would be preferable to avoid underestimating MB flow velocity. In this case, if the MB velocity is underestimated, it doesn't affect the MB concentration but can underestimate the network filling and optimal MB concentration. A faster MB covers more distance per time unit, therefore populates more vascular space per time unit. Underestimating optimal MB concentration is on the safer side in the sense that it reduces possible MB interactions.
As for ultrasound signals, they were generated using either readily available or realistic ultrasound probes acquisition parameters. The rationale behind this was to assess current limitations and seek optimal ULM parameters. White Gaussian noise was added to radio frequency data to improve realism, but no tissue was modeled. Since PSF size and shape greatly influences our ability to precisely detect single MB events, any alteration to the PSF from higher noise or lower contrast from neighboring tissue could have an impact on MB localization and ultimately maximal MB concentration. Other studies with tissue modelling are therefore to be conducted with tissue modelling as in [53].
The computation burden of the generation of RF signals was proportional to the number of simulated MBs. Combined with the burden of 3D beamforming increasing with the frequency cubed, this made high-frequency and high-concentration simulations significantly longer.
In terms of imaging parameters, common planewave transmit sequences include planewave compounding of 5-12 tilted planewaves and a pulse length of 1-10 cycles [39], [55], [56]. A single planewave was simulated in this work to reduce computation time, and a pulse length of 4 cycles was selected to send more energy in the medium to increase contrast. We worked under the assumption that most of the impact of coherent compounding could be accounted for by the PSF size. A more complete study on the impact of compounding on ULM performance would be required To generate a correlation map, a cross-correlation was performed between IQ data and the PSF. The PSF was generated by a point-scatterer in the middle of the field of view. Since the PSF of a system is non-uniform throughout the field of view, a more accurate approach would consist of using different PSFs for each voxel position.
In Fig. 6, we can see that thresholding does not completely prevent MB interactions, as depicted in the increasing localization error as a function of concentration. Therefore, correlation threshold must be selected wisely to maximize MB detection events for a fuller network while minimizing strongly interacting MBs which results in blurred vasculature in density maps. Localization error can reach a plateau or decrease at higher concentrations. The latter seems to happen when the MB concentration is so high that the pairing assumption is violated. The pairing assumption was that the error of a MB's position can be estimated using the closest neighboring MB. In that case, the calculated MB error isn't reflective of the true error. At very high concentrations of MBs, the closest neighbor of a localized MB is most probably not the reference MB's position. Therefore, the Sørensen-Dice coefficient could be a more robust metric to evaluate vascular reconstruction for high concentrations. Moreover, for frequencies of 9, 12, and 15 MHz, more simulations at higher concentrations would be interesting to see the same behavior in lower frequencies, where the MB count declines at higher MB concentrations. Due to the computation burden, simulations were stopped at those frequencies when MB count reached a maximum.
All simulations were conducted on networks of mouse brain vasculature where MB flow reflected behavior from in-vivo data. However, the vascular networks were relatively small in terms of field of view (about 500 µm). This means that MB trajectories were limited to a small range of vessel sizes (mainly small vessels). This reduced vessel size range and increased inter-network variability, which can be seen through results from Fig. 9, where MB count variability is not proportional to the inverse of PSF volume. This also implies that the network size in wavelengths is about 5λ at 15 MHz and about 1λ at 3 MHz. At 3 MHz, the resolution limit is about λ/10 or 50 µm, which is the size of the largest vessel in the network. Therefore, at that frequency, only one vessel can be reconstructed. However, even with the small dataset we had access to, we still wanted to test the limits of the ULM processing with this validation framework. The proposed simulation framework provides access to both quantitative metrics and a visual assessment of the vascular reconstruction compared to the reference network. This provides insight on aspects that can be harder to capture with metrics only. Future work with limited access to large datasets could mitigate the issue by artificially dilating small networks isotropically to span vessel diameters that are representative of a whole organ's distribution.
In the medium, MBs were spread uniformly to artificially cover a large portion of space, but realism could be improved by having access to a larger network in terms of field of view and ultimately a whole brain or whole organ vasculature. A vascular network that exhibits whole organ vascular size distribution will increase metrics repeatability. Since lowfrequency probes are generally used to allow tissue imaging at great depths (order of 10 cm), a larger dataset would also allow to account for the compromise between imaging depth and resolution.
In subsequent studies, acquisition parameters such as number of compounded angles [44], number of pulse cycles, pulse frequency modulation (chirp) [57], SNR, probe elements geometry, MB size distribution, MB nonlinear response, could be interesting to better understand the underlying limitations of ULM and progress towards finding optimal imaging parameters.

V. CONCLUSION
In conclusion, we have developed a framework capable of characterizing ULM image formation algorithms in 3D. We described how to obtain virtual anatomically-realistic MB positions from 2PM data. We have shown and quantified the impact of concentration on MB localization and acquisition time. We quantified the link between maximal MB concentration and PSF volume. Moving forward, standardized anatomically-realistic MB simulations could become a useful tool in the validation of an ULM imaging setup and the PSF volume a straightforward metric for a specific probe's acquisition parameters since it would define the appropriate MB concentration. In practice, we hope that these results can be used to achieve optimal vasculature reconstruction at lowest acquisition time for a specific imaging setup.

VI. CODE AVAILABILITY
All code used for the generation of the MB flow is available on GitHub at https://github.com/provostultrasoundlab/ hatimb-particle_flow_simulator. Example code to simulate and reconstruct RF signals is available in the MATLAB MUST toolbox https://www.biomecardio.com/MUST/.

VII. DATA AVAILABILITY
The main data supporting the results in this study are comprised within the article. Example steady-state MB flow datasets in six networks are available from the FRDR repository at https://www.frdr-dfdr.ca/repo/dataset/ae5706dc-22fb-41c2-be22-082d893fcb9a.