Understanding the effects of microbubble concentration on localization accuracy in super-resolution ultrasound imaging

Super-resolution ultrasound (SRUS) through localising and tracking of microbubbles (MBs) can achieve sub-wavelength resolution for imaging microvascular structure and flow dynamics in deep tissue in vivo. The technique assumes that signals from individual MBs can be isolated and localised accurately, but this assumption starts to break down when the MB concentration increases and the signals from neighbouring MBs start to interfere. The aim of this study is to gain understanding of the effect of MB–MB distance on ultrasound images and their localisation. Ultrasound images of two MBs approaching each other were synthesised by simulating both ultrasound field propagation and nonlinear MB dynamics. Besides the distance between MBs, a range of other influencing factors including MB size, ultrasound frequency, transmit pulse sequence, pulse amplitude and localisation methods were studied. The results show that as two MBs approach each other, the interference fringes can lead to significant and oscillating localisation errors, which are affected by both the MB and imaging parameters. When modelling a clinical linear array probe operating at 6 MHz, localisation errors between 20 and 30 μm (∼1/10 wavelength) can be generated when MBs are ∼500 μm (2 wavelengths or ∼1.7 times the point spread function (PSF)) away from each other. When modelling a cardiac probe operating at 1.5 MHz, the localisation errors were as high as 200 μm (∼1/5 wavelength) even when the MBs were more than 10 wavelengths apart (2.9 times the PSF). For both frequencies, at smaller separation distances, the two MBs were misinterpreted as one MB located in between the two true positions. Cross-correlation or Gaussian fitting methods were found to generate slightly smaller localisation errors than centroiding. In conclusion, caution should be taken when generating and interpreting SRUS images obtained using high agent concentration with MBs separated by less than 1.7 to 3 times the PSF, as significant localisation errors can be generated due to interference between neighbouring MBs.


Introduction
Super-resolution imaging, based on the localisation of sparsely activated fluorophores has transformed optical fluorescence microscopy (Betzig et al 2006, Hess et al 2006).More recently, its acoustic counterpart, also known as ultrasound localization microscopy (ULM), has been shown to be capable of separating two acoustically responsive targets at a resolution beyond the diffraction limit both in vitro (Desailly et al 2013, Viessmann et al 2013), and in vivo (Christensen-Jeffries et al 2015, Errico et al 2015).It has been successfully applied in a range of animal models (Siepmann et (Opacic et al 2018), in limbs (Harput et al 2018), in the brain (Demené et al 2021), in liver (Lok et al 2022), and in lymph nodes (Zhu et al 2022).A recent review on Super-resolution ultrasound (SRUS) can be found in (Christensen-Jeffries et al 2020).
For microbubble (MB) and localisation-based SRUS/ULM imaging, a significant challenge is the long data acquisition time required to accumulate a sufficient number of isolated MB signals to generate an image, which can range from seconds to minutes (Christensen-Jeffries et al 2019, 2020).To date, this has been addressed by increasing the concentration of MBs injected and hence generating more MB signals per imaging frame (Christensen-Jeffries et al 2019).Higher concentrations, however, mean that MBs will be closer to each other, which makes it more challenging to isolate them, especially as the MB signals may generate coherent interfere patterns in the image, or, at closer distance, their oscillations may start to interfere with each other (Yusefi and Helfield 2022).The radiated acoustic signals from MBs start to interfere when the MBs approach each other and their point spread function (PSF) starts to overlap in space.Unlike linear scatterers within tissue, MBs oscillate, and their oscillations can be highly nonlinear when driven at resonance, scattering harmonics of the transmitted signal (Versluis et al 2020).
This non-linearity has been utilized to generate contrast-specific images using methods such as pulse inversion (PI) and amplitude modulation (AM), which are now widely used clinically (Stanziola et al 2016, Christensen-Jeffries et al 2020).Furthermore, it is well known that MB signals can be significantly impacted by the ultrasound frequency, the acoustic pressure amplitude, and the ultrasound pulse sequence used (Tang and Eckersley 2007, Tang et al 2010, Christensen-Jeffries et al 2020).However, how these factors contribute to the ultrasound image of a MB when it approaches another MB, and how the localisation accuracy is affected in ULM is currently not well understood.
The aim of this study is to obtain a better understanding of what happens to a ULM image as the MB-MB distance decreases, using a combination of acoustic field and nonlinear MB dynamics simulations (Lerendegui et al 2022).Specifically, the impact of MB separation, MB orientation, imaging mode (Bmode or contrastspecific mode), MB sizes, ultrasound frequency, and pulse amplitude, on the images and on MB localisation accuracy in ULM are investigated.

Methods
2.1.The ultrasound imaging simulation software Bubble flow field (BUFF) (Lerendegui et al 2022), was the ultrasound imaging software used for simulation.It uses a linear ultrasound field simulation back-end (Jensen andSvendsen 1992, Jensen 1996), coupled to a MB dynamics model (Marmottant et al 2005) to generate synthetic ultrasound images, using the following scheme: the acoustic pressure signal incident upon the MBs is first calculated.This pressure is then used as an input to the MB dynamics simulation.Finally, the acoustic simulator is used to obtain the propagation of the radiated pressure from the MB back to the transducer.Ultrasound images were reconstructed using the delay and sum method.Synthetic random coloured Gaussian noise, i.e. white noise filtered by the transducer bandwidth, was added to the beamformed image data with a signal to noise ratio (SNR) ratio of 20 dB.
Two transducers with different central frequencies were modelled: one, the GE M5Sc-D probe (GE Healthcare), a low frequency (LF) clinical cardiac transducer with a centre frequency of 1.75 MHz, and the other, the GE L3-12-D probe (GE Healthcare), a high frequency (HF) transducer with a centre frequency of 6 MHz.

Bubble sizes and simulation parameters
For this study, the clinically approved agent Sonovue/Lumison (Bracco) was simulated using the fitted parameters displayed on table 1 (Gorce et al 2000, Marmottant et al 2005).
As SRUS is mainly focused on imaging microvascularture, we have used 1 atm as the ambient pressure parameter (P0).Normal pressure in the main arteries can reach 1.15 atm, but this is reduced in the microvessels to values very close to atmospheric pressure.
The contrast agent has a radius distribution shown in figure 1 (Gorce et al 2000).An average human capillary internal diameter of 10 μm is assumed, so no MBs larger than this are considered.For each simulation, two MBs were used, one referred to as the modal MB, having the modal radius of the size distribution (Modal Radius R = 1.5 μm), and the other referred to as the resonant MB, having the resonant radius.Note that the resonant size depends on the ultrasound frequency and hence is different for the LF and HF simulations (Resonant Radius either R = 2.745 μm at 1.75 MHz or R = 1.115 μm at 6 MHz).The resonant MB sizes were calculated approximately using equation (1), which calculates the linear natural frequency (Versluis et al 2020).
The scattered pressure from the bubble was calculated using equation (2).Note that this equation ignores the kinetic term that goes as r 1 4 and only the far field scattered pressure is considered this way.The inter bubble distances considered were all of the order of 1 wavelength and thus it was deemed acceptable to neglect the r In each simulation two MBs were positioned at a single imaging depth corresponding to the elevational focal depth of the transducers (77 mm for LF and 22 mm for HF respectively), and their lateral positions across the xaxis were varied (inter MB distances: 0.1−6 wavelengths for HF and 0.1-12 wavelengths for LF), to study the impact on the image and localization accuracy.The same simulation was repeated with their axial positions varied across the z-axis.Linear scatterers and isolated single MBs were also simulated as controls.
Both B-mode and AM contrast enhanced ultrasound (CEUS) images were generated.For each B-mode image, 9 plane wave of 1 cycle were transmitted with a steering angle range of [−4, 4]°with a step of 1°.Received signals were coherently compounded.For AM images one full and one-half amplitude pulse were used.Three different amplitudes were investigated at each frequency corresponding to mechanical indexes (MIs) of: 0.05, 0.1, 0.2.Random noise was beamformed separately and then added to the beamformed RF data before generating the final B-mode and AM images.The SNR ratio was 20 dB.Adding noises this way represents the same result as adding them to the RF data prior to beamforming but is computationally more efficient.
For the simulation, multiple parameters were varied to represent common imaging scenarios.More than 24 000 simulations were performed, using all possible combinations of the choices for the following parameters: • Transducer: low frequency, high frequency   • Bubble types: resonant-resonant, resonant-modal, modal-modal • Mechanical index: 0.025, 0.05, 0.1, 0.2 • Target type: microbubble, linear scatterer • Interaction type: two bubbles, isolated bubbles • Distances: 128 steps.

Bubble isolation and localization
MB images were first normalised with respect to max pixel intensity across frames.To remove noise and isolate the main lobe of MB signals for localisation, an image threshold for each sequence was determined empirically by selecting the strongest side lobe and tail of the PSF.This entailed an inevitable trade-off to minimize false positive detections since higher thresholds will have the disadvantage of weak MBs not being detected in cases with large MB mismatch (resonant-modal), while lower threshold values produce less separation between MBs.
For the HF data, 0.2 was used as threshold for all the MB pairs.For the LF data, the main lobe of modal MB image was not much stronger than the side lobe of resonant MB image.Therefore, only resonant-resonant and modalmodal cases were analysed using threshold values of 0.2 and 0.3 respectively.Isolated lobes with dimension less than 10 pixels (0.57λ for HF probe) or 20 pixels (0.88λ for LF probe) were rejected to reduce the effect of lacking signal on the precision of localisation.Three commonly used localization methods, i.e. centroiding (Christensen-Jeffries et al 2015), normalized cross-correlation (Song et al 2018), and Gaussian fitting (Heiles et al 2022) were chosen as we found these methods to work well for asymmetrical PSFs with multiple local maxima and low-SNR images.
In the centroiding method, MBs were localised by the intensity-weighted centroid of the isolated main lobes.In the normalised cross-correlation method, images of single linear scatters at the centre were regarded as the PSF of the simulated imaging system.After interpolating PSFs and MB images with an up-sampling factor of 10, achieving a pixel resolution of around 0.01λ, the PSF of each frequency was then cross-correlated with the MB main lobe and the peak of the correlation coefficient map was used to localise MBs.In the Gaussian fitting method, a 2D Gaussian function with six degrees of freedom-two mean values: μ x , and μ z , two standard deviations: σ x , and σ z , the amplitude A, and an offset C-was fitted to the main lobe and the peak of the fitted Gaussian function was used to localise the MB.
Gaussian fitting and normalized cross-correlation were done on the same MB images segmented for calculating centroiding to avoid the effect of MB isolation or image segmentation on localization error evaluation.
The locations obtained by centroiding and Gaussian fitting, were rounded by the same pixel size that is used for the interpolation of the normalised cross-correlation.

Quantification of localisation errors
The super-localisation methods were applied to all the simulations, obtaining estimated MB coordinates.The localisation error was used as the metric to evaluate which method deals better with the interference between neighbouring MBs.This metric is calculated as the distance between the localised position of the MB with and without the neighbour being present.Mean Localisation errors were calculated by averaging at three levels: across simulated distances, distances and directions, and distances directions and MB types.The first levels show the expected average localisation error for specific scenarios and were used to see if errors change with the direction and MB mismatch.During experiments, we expect to see all kinds of MBs interfering with others in any direction, which is why the third level, a more general level, was used to rank the performance of localisation algorithms to determine which one is affected the least.For this ranking, a lower number corresponds to better performance, i.e. smaller average localisation error.

Results
In this section, some example images of two MBs approaching each other generating interference patterns are presented with different imaging parameters.These parameters include frequency, imaging mode (B-mode or AM), and MI.MB images are also compared to those of linear scatterers.The influence of the interference by a neighbouring MB on localization was then quantified.
Figure 2 shows how the AM image of one MB is affected by the presence of another MB as they approach each other vertically and horizontally with both the simulated HF and LF probes.The MB signals start to interfere, even when the MBs are multiple wavelengths apart.With the HF probe, the two MBs cannot be visually separated when they are ∼2λ apart.With the LF probe, two main lobes start to overlap at ∼5.5λ.This difference is expected from the different lateral resolution of the two probes.Moreover, main lobe of the modal MBs is not much stronger than the side lobes of the resonant MB when two MBs are placed horizontally or even weaker than the tails of resonant MBs when two MBs are placed vertically, which makes it challenging to include the main lobes of modal MBs and exclude side lobes/tails of resonant MBs for MB isolation at the same time.
Figure 3 shows the B-mode images for the same pair of MBs.While there are some differences between the B-mode and AM images, the qualitative findings are very similar, i.e. that MBs can only be localized when they are many wavelengths apart, and for the LF probe the separation needs to be much larger.
Figure 4 compares the B-mode images of two MBs and two linear scatterers.It can be seen that there are some differences in their B-mode image appearance, and these differences, as shown in the images, are larger at lower frequency.
Figure 5 shows the impact of MI on the AM images.While increased MI could increase the signal intensity, it also increases the overlap between the two MBs signals making them harder to isolate.
Figure 6 shows the localisation error obtained with normalised cross-correlation for two MBs as a function of their separation distance for the HF and LF probe along x-direction.It can be seen for HF probe that when the MBs are ∼1.5λapart, the localisation error is just under 0.15λ.However, when they are closer, the simple thresholding algorithm failed to separate the two MBs, represented by the grey shaded area in the plot.It is also interesting to see that the localisation error seems to oscillate as a function of distance, likely due to the alternation between constructive and destructive interference as the MBs move towards each other.For the LF probe, the localisation for resonant-modal MB pairs is not displayed here as the main lobes of modal MBs cannot be effectively isolated with the threshold set as the strongest side lobes and tails.Lateral localization error oscillates with distance between MBs, as in the HF cases.More importantly, the grey area grows significantly, which means when using the LF probe it is much harder to separate the two MB signals and is consistent with the observation in figure 2. This is related to the small aperture of the clinical cardiac probe and the correspondingly poor lateral resolution.Figure 7 shows the equivalent plot for B-mode images with similar findings.
Figure 8 shows the equivalent plot for vertically placed MBs.Compared to figure 6, the localization error is larger, which means the interference along the axial direction is stronger than that along the lateral, even if the axial dimension of the PSF is usually considered smaller than the lateral dimension.Again, the results for the B-mode images are qualitatively similar, as shown in figure 9.
Tables 2 and 3 show the average localization errors in three different levels for the LF and HF simulations respectively.
To summarize, the results show that the MB signals start to interfere even when the MBs are multiple wavelengths apart, and the interferences increase with MI.The localization errors due to this interference oscillates as a function of distance and can be up to 1/5 of a wavelength, until the two MBs cannot be separated  anymore.Cross-correlation and Gaussian fitting seem to be able to localize the MBs slightly better, if the MBs are isolated in the first place.The signals from a resonant MB can often overwhelm those from a non-resonant MB, and the side lobes of the resonant MB can often be higher than the main lobe of the non-resonant MBs.

Discussion
In SRUS, high MB concentrations are increasingly being used to reduce the total acquisition time, and several However, there is still a lack of understanding of the impact on the image and the localisation of a MB when there is another MB nearby.This study used simulations of both ultrasound field and MB dynamics to gain some Table 2. HF Localisation errors for the three localisation methods, ranked for all types of MB mismatch, orientation and image generation method.Res corresponds to a resonant MB, and Mode corresponds to a modal MB. further understanding of this topic.Over 24 000 synthetic ultrasound images of MB at different distances were generated and processed using several ULM localization algorithms.The results were analysed to understand what happens to a SRUS image as MBs get closer.The simulation parameters were chosen to cover a large range of possible scenarios.
The results show that MB signals can start to interfere with each other and generate interference fringes even when the MBs are separated by many wavelengths, and that this negatively impacts localization accuracy.There also seems to be a certain MB separation distance below which the localisation cannot separate the two MBs.This distance limit is 1.7-3 times larger than the corresponding PSF in this study using an available intensity thresholding algorithm.
Beyond this limit, the localisation error due to a nearby MB oscillates with distance, likely due to constructive and destructive interference between MB signals, which is also affected by the direction between two MBs.In this study it was found that when using AM at 1.5 MHz, localisation errors up to 0.2λ can be generated when two MBs are 10λ apart for horizontally placed MBs and up to 0.5λ for vertically placed MBs.This error increased with decreasing MB separation.
The results also show that localisation accuracy also varies with MB size, ultrasound probe characteristics, and pulsing sequence; with larger localization errors being generated at lower frequencies, when MBs were at resonance, and when AM was used.Again, it is hoped that these findings may help to guide the development of improved localization and tracking algorithms.For example, a larger uncertainty can be allocated to a MB localization if there is another localization within a certain distance and varied with the orientation between MBs.With adaptive localisation uncertainties, the position of the MB might be updated to be closer to it real location after Kalman filtering when using flow dynamics.
It should be noted that this study focused on evaluating localization errors using the different localization methods, and not on the methods' ability to isolate individual bubbles.Bubbles are segmented through thresholding in the same way, before Gaussian fitting, centroiding, or cross-correlation methods were applied.While cross-correlation has the least localization errors, it is likely to be less able to isolate closely spaced bubbles than e.g. a Gaussian fitting based method.
Table 3. LF Localisation errors for the three localisation methods, ranked for all types of MB mismatch, orientation and image generation method.Res corresponds to a resonant MB, and Mode corresponds to a modal MB.
One key challenge in localising MBs is that there is a large variation in the MB signals depending on their size and whether they are resonant.Smaller MBs and/or MBs that are not resonant usually produce weaker signals that can be overwhelmed by the signal from a nearby larger or resonant MBs, as shown in e.g.figures 2, 3, 5.
In this case not only the weak MB can be missed, but its overlapping signal can also contribute to localisation errors of the large MB nearby.The weak MB signals also pose challenges in the presence of high noise levels.Modelling of such weak MB signals could be helpful for algorithms to distinguish such weak MB signals from side lobes of strong signals and detect these weak signals by using algorithms better than the used simple thresholding.
Furthermore, the B-mode and AM images generated using linear scatterers and MBs show some small differences.These should be considered when using data e.g. to train a learning algorithm, since most clinical systems that image MBs in contrast-specific mode often use either AM, PI or a combination of both.
The echoes from microbubbles are heavily affected by their shell properties (Versluis et al 2020).In this study we used a set of bubble shell parameters as in a previous reference (Marmottant et al 2005).However, a change to the shell properties may impact the resulting images.E.g. a more elastic shell could increase the amplitude of the echoes and increase the difference in appearance between a bubble and a linear scatterer, which in turn could require a farther separation of bubbles for them to be individually detected.
Acoustic coupling between MBs was not taken into account in the simulations.Previous studies have suggested that such coupling is only significant when MBs are less than 20 diameters apart, which is about 0.5-0.2λ in this case (Harkin et al 2001).The results above, however, show that the image fringes caused by the coherent interference have a significant impact on localization accuracy even when the two MBs are relatively far apart (∼100 diameters).
It has been shown that echoes from microbubbles are affected by its surrounding tissue structures.E.g. a bubble in a microvessel could have non-symmetrical and reduced oscillation (Doinikov and Bouakaz 2015), which could change their image appearance.This will increase the variations and pose additional challenges to the SRUS bubble identification algorithms.
In this study, a simple empirical thresholding algorithm was used for segmenting MB signals from noise and side lobes.There are more advanced algorithms available, such as deconvolution on single frame (Yan et al 2022) or deconvolution exploiting flow dynamics (Solomon et al 2019), which might improve the detection and separation of MBs at smaller distance.Furthermore, this study only looked at separating MBs in lateral and axial directions, which are the more challenging cases given the asymmetry of the PSF, and it would be expected that the trends identified in this study would still apply to MBs approaching each other in other orientations.

Conclusion
Over 24 000 simulated ultrasound images of two neighbour MBs were generated, with parameters chosen to cover a range of possible scenarios: 128 different distances, two orientations, three types of MB mismatch, two MB radii, three MI, two probes: HF and LF, and two image modalities: AM and B-Mode.Images were processed using 3 ULM localization algorithms: centroiding, normalized cross-correlation, and Gaussian fitting.The results were analysed to understand the effect that MBs getting closer and closer will have on ULM images.
The results show that: • MBs give different images than linear scatters, and this difference is greater on the LF simulation.
• MB radii affect the MB images, and the resonant MBs generate brighter images.This effect is more significant at the lower frequency.
• Brightness of MB images increases with MI for all the MB sizes and both frequencies.
• Wave interference happens when MBs are close, which is reflected in the B-mode and AM images and localisation errors on these two kinds of images.
• Localisation errors can also be affected by MB sizes and the orientation between two MBs.
• Interference from neighbouring MBs can produce localisation errors of up to 0.2λ, even if the distance between them is large (∼ 10λ).
• Normalised cross-correlation performs better than centroiding and Gaussian fitting when localising MBs using a thresholding method to discard noise and side lobes.
Considering factors that affect localisation errors can be helpful when developing further better localisation algorithms.The findings of this research also provide a solid foundation for defining appropriate MB concentrations for ULM imaging.

Figure 2 .
Figure 2. AM images showing the visual artifacts and fringes generated by two different MBs (one resonant and one modal) at different separations.MB distance decreases from top to bottom.Red crosses indicate ground truth MB locations.All images are displayed at 50 dB dynamic range.The first two columns show the HF simulation with MBs approaching each other in the x and z directions, with a MI = 0.1.The two columns on the right show the LF simulation with MBs approaching each other in the x and z directions, with a MI = 0.05.Note that interference patterns occur even when the two MBs are many wavelengths (WL) apart, and the non-resonant MB signal is barely visible when a resonant MB is nearby, as shown in the low frequency cases.

Figure 3 .
Figure 3. B-Mode images showing the visual artifacts and fringes generated at different MB distances, with the same configuration as figure 2.

Figure 6 .
Figure 6.Localisation error as a function of MB horizontal separation distance at high frequency (MI = 0.1) and low frequency (MI = 0.05) using AM for different MB size combinations.Only one sample was collected per data point.Top row: lateral error; middle row: axial error; bottom row: absolute error.Types of simulated MBs are labelled at bottom row of each column.res is the abbreviation of resonant MB.Grey shaded areas indicate that at this distance the two MBs cannot be separated anymore by the use thresholding method, i.e.only one MB is detected, or any localisation error is over half a wavelength.The percentage of minimum separation distance to either the lateral or axial full width of half maximum of PSF is presented as the value above the black arrow in the bottom lefthand corner of each panel.The main lobe of mode MB is not stronger than the strongest side lobe of resonant MB under low frequency, and thus the resonant-modal column is not displayed here.

Figure 5 .Figure 8 .
Figure 5. Impact of MI on the AM images of two resonant MBs: images of high and low frequencies are normalised by the maximum intensity in corresponding frequency respective.Note all images are displayed at 50 dB dynamic range.

Figure 9 .
Figure 9. Same with figure 8, except for that two MBs are imaged by B-mode.

Figure 7 .
Figure 7. Same with figure 6, except for two MBs imaged by B-mode.
algorithms have been developed to localize and track MBs at high concentrations.Linear filtering approaches have been proposed by Elder (Bar-Zion et al 2018), (Harput et al 2017) and (Huang et al 2020) to separate individual MB signals based on either their frequency responses or spatial-temporal features.More recently, deep learning has shown promise in localising MBs at higher concentrations (van Sloun et al 2019, 2021) provided the training datasets are representative of the imaging tasks.While it is difficult to control the concentration of MBs in vivo after their injection, MB signal concentration can be much better controlled through acoustic switching of phase change contrast agents in the so-called AWSALM approach (Zhang et al 2018, 2019a), which has been demonstrated to achieve real-time super-resolution (Riemer et al 2022), on a cross-tube micro-flow phantom (Zhang et al 2018), and on a rabbit kidney model (Zhang et al 2019b).