Super-resolution ultrasound imaging

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


INTRODUCTION
Our circulatory system is so vital that the loss of blood flow is one of the key parameters defining death. This vessel network created by nature comprises billions of vessels that carry fundamental nutrients, hormones and gases at distances longer than simple diffusion in large living beings (Pugsley and Tabrizchi 2000). Would one lay out all of the arteries and veins and the 40 billion capillaries in one human adult, they would reach more than 100,000 km or two times the circumference of the earth. The tiniest components of our vasculature, the capillaries, are less than 10 mm in diameter (Lenasi 2016) or about a tenth of the diameter of a human hair. Some capillaries are even smaller in diameter than blood cells, forcing cells to distort their shapes to pass through.
From a biomechanical point of view, the circulatory system is also a piece of extraordinary machinery ensuring rapid transport and complete distribution of blood at meters per second in our largest arteries down to less than 1 mm/s in the capillaries feeding the vast territory of our organs at microscopic scales. To achieve this amazing feat, the 3-D geometry of our vasculature and the rigidity of each arterial segment are carefully optimized. The arterial stiffness also adapts itself transiently after a load or arterial pressure changes. This highly nonlinear elastic nature of the arterial walls is essential to effectively damp the large oscillations in blood flow coming from the heart. It provides for a better flow homogeneity in tiny blood vessels distal in the arterial vascular tree. As a consequence, pathologic changes of the mechanical properties of arteries strongly affect the transmission of blood to tiny vessels (Webb et al. 2012). Today, although the field of regenerative medicine and biomaterials is rapidly progressing, mimicking the complete vascular system with optimal structural and functional properties remains challenging.
In some noble organs such as the brain, the complexity extends to an even higher level, as tiny vessels are intimately connected and communicating with the neuronal system via the glial system and particularly astrocytes and pericytes (Iadecola, 2004). Such communication allows for precise coupling between neuronal activity and blood flow, which ensures that activated brain regions are properly nourished in oxygen, glucose and other nutrients. This neurovascular coupling is the basis for functional magnetic resonance imaging (MRI) (Kim and Ogawa 2012) and functional ultrasound (Deffieux et al. 2018). In the brain, as in other organs, the microvasculature is a dynamic system that adapts to the constantly changing metabolism of surrounding cells.
At such a microscopic level, large territories of our knowledge remain unexplored mainly because of the lack of imaging methods providing non-invasiveness, microscopic resolutions, a macroscopic field of view and sufficient temporal resolution for dynamic imaging. From a fundamental point of view, it is, for example, striking to note that scientists do not fully understand the functioning of the human placental vascular system and exchanges between maternal and fetal blood systems (Mayo 2018).
Many discoveries have also shed new light on the major importance of our vascular system in various disease processes, ranging from cancer, to diabetes, to neurodegenerative diseases such as Alzheimer or Parkinson diseases (Zlokovic, 2011;Stanimirovic and Friedman 2012). For example, it has been known for more than 40 y that angiogenesis, the development of new blood vessels, is a hallmark of solid tumors (Folkman 2006). Angiogenesis is driven by tumors that outgrow their host tissue's native blood supply, with the result of the release of pro-angiogenic factors that locally stimulate increased microvascular development to feed the growing malignancy. Pathologic angiogenesis is differentiated from the normal microvascular structure by the lack of hierarchical branching, the presence of tortuous and erratically shaped vessels and immature and leaky vessels (Jain 1988).
Dementia also has a microvascular component, as it is more likely to be present when vascular and Alzheimer disease lesions co-exist (Jellinger 2008). In elderly individuals, the association between stroke and Alzheimer disease increases in patients with vascular risk factors (Iadecola et al. 2010;Gorelick et al. 2011). Vascular endothelial growth factor, one of the most potent mediators of angiogenesis, can be envisioned as a potential therapeutic for neurodegenerative disorders (Storkebaum et al. 2004).
The diagnosis of these diseases would benefit greatly if the microcirculation could be characterized in each tissue, organ and patient. Indeed, pathologies leave first their signature in tiny vessels before becoming detectable much later in larger vessels by a dramatic domino game.
Unfortunately, the naked eye cannot resolve the vessels smaller than 100 mm forming the microcirculation. Moreover, most of these vessels lie beyond the penetration depth of coherent light in tissue. Histopathology can be performed after a biopsy or a surgical resection and remains the gold standard for cancer diagnosis. However, such approaches are limited by their invasiveness in the clinical setting or pre-clinical research. Microvascular parameters linked to angiogenesis, such as microvascular density and intercapillary distance, are obtained by observing and measuring stained vessels on highly magnified thin slices under a microscope. Several types of tissue staining can help reveal microvessels specifically such as anti-CD31, CD34 and von Willebrand factor (Weidner 1995;Marien et al. 2016).
Within a few hundred microns depth, microscopy can also be applied directly on the skin or mucosa to observe its microcirculation. Techniques such as orthogonal polarization spectral imaging and sidestream darkfield imaging can extract specifically the light from blood and provide a map of blood vessels under the surface (Leahy 2012). Flow can also be assessed with laser Doppler either at a single point or on an entire map. Additionally, retinography can assess the evolution of the microcirculation, in diabetic patients for example, by exploiting the clear window provided by the eye (Pieczynski and Grzybowski 2015).
Various modalities are able to reach microscopic resolutions such as two-photon imaging (Soeller and Cannell 1999) or optical coherence tomography (Jia et al. 2012) at the cost of a limited field of view and penetration depth. Other approaches based on tissue clearing coupled to light-sheet microscopy lead to high-resolution volumetric imaging of the microvasculature in organs but are limited to dead tissues (Ertuk et al. 2012;Ragan et al. 2012). Additional approaches such as photo-acoustics (Wang and Hu 2012) and functional ultrasound (Deffieux et al. 2018) based on ultrafast Doppler (Bercoff et al. 2011;Demene et al. 2016) recently improved our ability to image small vessels (of the order of 100 mm in diameter) but fail to reach microscopic resolution scales.
Moreover, for human diagnostic or animal imaging in-depth, it is necessary to exploit modalities that can explore entire organs, at a depth beyond 10 cm. Wellknown medical techniques, such as MRI (Williams et al. 1992), computed tomography (CT) (Miles 1991), nuclear imaging (Underwood et al. 2004) and ultrasound (Cosgrove and Lassau 2010) all have versions that are sensitive to blood perfusion. Perfusion CT works by following a bolus of injected iodinated material in different parts of the organ. Parameters such as blood flow, blood volume, time to peak and mean transit time can be extracted from the bolus curves. In a similar fashion, perfusion MRI can also be performed with a contrast agent such as gadolinium chelate complex. Furthermore, using arterial spin labeling, MRI can even provide information on perfusion without contrast agent injection (Peterson et al. 2006). With the use of injected radionuclide, single-photon emission CT (SPECT) can be made sensitive to perfusion, which is clinically applied to the cardiac muscle for instance. However, these techniques provide only broad generalities on the microcirculation in each imaging voxel. None of them can define the microvascular architecture itself because these macroscopic modalities are limited in resolution to the submillimeter and millimeter scales.
In particular, ultrasound imaging is limited in resolution by diffraction to the scale of its wavelength (wavelength = speed of sound/frequency, a 5-MHz ultrasound wave in tissue has a 300-mm wavelength). It relies on the echo of tissue owing to variation in compressibility and density to reconstruct anatomic images. In Doppler mode, it is sensitive to blood vessels through the motion of red blood cells serving as scatterers. However, small vessels have a limited number of weak scatterers, which are also moving slowly, making it particularly difficult to distinguish vessels from tissue motion. Generally, vessels with blood velocities below 1 cm/s are difficult to distinguish, making Doppler ultrasound a poor imaging method for the microvasculature. Even with recent advances exploiting ultrafast plane-wave imaging and spatiotemporal filters (Bercoff et al. 2011;Demene et al. 2016), which improved drastically the sensitivity of Doppler ultrasound, micro-arterioles and microvenules remain invisible to Doppler ultrasound.
As in other medical imaging techniques, ultrasound imaging can be made sensitive to unresolved microvessels by the introduction of contrast agents (Cosgrove and Lassau 2010). These agents are microbubbles, smaller than capillaries, which are injected intravenously and flow within the entire vasculature for about 3 min (Ferrara et al. 2000;Ferrara et al. 2007;Blomley et al. 2001;Burns and Wilson 2006). Microbubbles act as resonators with a resonance frequency in the range 1À15 MHz, vastly increasing their scattering coefficient in the clinical frequency range. Moreover, microbubbles re-emit ultrasound in a non-linear fashion, providing a tool to separate them from tissue (Frinking et al. 2000;Stride and Saffari 2003;Dollet et al. 2008). The presence of these contrast agents highlights the vasculature, including the capillaries, as the ultrasound scanner is also sensitive to slowly moving microbubbles. One great advantage of microbubbles for perfusion imaging is that they are entirely intravascular because of their micrometer size. Hence, after injection, the only compartment to be taken into account for the calculation of parameters such as the mean transit time is the vasculature. Crucially, in contrast to optical agents, microbubbles can be detected deep within the body, making them advantageous as a clinical contrast modality. Furthermore, contrast-enhanced ultrasound (CEUS) scans are already an established and routine clinical procedure in many clinics around the world, making fast clinical translation a real possibility.
Unfortunately, conventional ultrasound remains limited by resolution in the same way as MRI, CT or SPECT. The extracted parameters are linked only indirectly to modifications in the microcirculation. If a medical imaging technique directly maps microvessels, it would provide a revolutionary wealth of information, bridging the gap with histopathology. For instance, such a technique could measure directly vessel density, interdistance, size, unique flow pattern, tortuosity or fractal factor.
The extensive work on microbubble imaging has recently inspired a new technique that has caused an important rupture in a fundamental characteristic of ultrasound: its resolution. Introduced 10 y ago, super-resolution ultrasound can improve the resolving power of ultrasound imaging by a factor of 10 with respect to the diffraction limit (wavelength/2).
This review describes super-resolution ultrasound imaging as it is conceived by several groups, which introduced several of the precursor works in the field. It will first detail its origin and its technical aspects, as well as define its key concepts and technical aspects. It will detail both ultrasound localization microscopy and other approaches based on fluctuations imaging. The latter section will discuss the current and future applications for oncology and neurology. 1992; Couture et al. 2018). The goal of super-resolution is to separate echoes coming from sources closer than the classic diffraction limit. Such a quest was performed in parallel to the improvement of resolution through the increase in acquisition frequency (Lockwood et al. 1996).
Approaches such as near-field imaging (Shekhawat and Dravid 2005) were found to differentiate subwavelength sources as the resolution close to a probe is proportional to the distance with respect to the object rather than the wavelength . However, in the body, organs are several centimeters deep, which could be a hundred wavelengths away. A far-field approach is thus required for medical imaging.
In the far field, refocusing on close individual sources could be performed when a limited number of them were present (Blomgren et al. 2002;Lehman and Devaney 2003;Prada and Thomas 2003). A precise knowledge of the source could also allow subwavelength imaging (Clement et al. 2005). However, a limited number of sources or strong a priori knowledge is not applicable to conventional B-mode imaging, which observes tissue formed by a multitude of scatterers at various scale: cells, organelles, fibers, and so forth.
Further rupture of the half-wavelength limit in ultrasonic imaging was inspired by new developments in optical microscopy. In 2006, fluorescence photoactivated localization microscopy, photoactivated localization microscopy (PALM) and stochastic optical reconstruction microscopy (STORM) were introduced, breaking the diffraction limit in optics by at an order of magnitude or more (Betzig et al. 2006;Hess et al. 2006;Rust et al. 2006). It relies on photoswitchable fluorescence sources and fast cameras, which take sequential images where only a subset of the sources is lit in each image. By isolating the sources closer to the wavelength, the interference of the wave they emitted could be avoided. Moreover, knowledge of the point-spread function (PSF) of the system leads to an extremely precise localization of isolated sources from their intensity map. By collecting thousands of subwavelength localizations, a picture with a resolution in the tens of nanometer could be obtained with a microscope using visible light. These developments were so revolutionary that they led to the attribution of the 2014 Chemistry Nobel Prize to Eric Betzig, Stefan Hell and William E. Moerner.
In 2010, an ultrasonic version of FPALM, now called ultrasound localization microscopy (ULM), was proposed (Couture et al. 2010). The fluorescent beacons were replaced by ultrasound contrast agents, and the cameras, by an ultrafast programmable ultrasonic scanner (Couture et al. 2009(Couture et al. , 2012. Beyond that, the same principle applied: the interference between different microbubbles was avoided by observing them sequentially so that isolated sources could be detected in each image. When the PSF on the radiofrequency channel data or on the beamformed image is known, a localization with a micrometric precision can be obtained for each microbubble. As these contrast agents are purely intravascular, the accumulation of these subwavelength localizations would yield a super-resolved map of the microvasculature. The ULM approach ( Fig. 1) was rapidly illustrated in vitro through imaging of a single micro-channel containing flow microbubbles (Couture et al. 2011). In parallel, the first in vivo application was reported by Siepmann et al. (2011), who described a technique to improve maximum intensity projection images of dilute microbubbles by implementing centroid detection. By 2013, four of our teams were already exploring super-resolution imaging. In vitro, Viessman et al. (2013) reported for the first time that ULM can distinguish two vessels separated by less than half a wavelength. Two 3-D super-resolution approaches were proposed, one with a 1.5-D array (Desailly et al. 2013) and another with a hemispherical array through human skull bone (O'Reilly and Hynynen 2013). Further in vivo applications were rapidly implemented afterward: Christensen-Jeffries et al. (2015) illustrated the application in the mouse ear, introducing super-resolved velocity mapping (Fig. 2); Errico et al. (2015) in the rat brain; and Lin et al. (2017a) in a cancer model. First-in-human demonstrations of the techniques with clinical scanners were provided for breast cancer by Opacic et al. (2018) and imaging in the lower limb by Harput et al. (2018).
These applications will be detailed in the latter sections of this review. At this point, it is important to explain precisely the principle of ultrasound localization microscopy. Non-ULM approaches will also be introduced.

GENERAL TECHNICAL ASPECTS OF ULM
Like optics, ultrasound faces a limit inherent to all wave-based imaging processes, where diffraction of the transmitted and received waves mean point sources become indistinguishable from one another when closer than approximately half the transmitted wavelength. Beyond this, interference of scattered sound results in acoustic speckle. After the revolutionary developments seen within the optical field, analogous approaches were proposed to exploit these same principles, but in the ultrasound field. Here, instead of utilizing molecules to provide the individual signal sources required, ultrasound contrast agents, called microbubbles, were proposed.
Thus, the super-resolution ultrasound process requires the introduction of a contrast agent into the body. Akin to its optical counterpart, it also requires the acquisition of a sequence of frames. A crucial principle within localization microscopy techniques is that by limiting the number of sources detected in each image, the responses do not interfere with each other. Under this constraint, the location of the underlying scatterers, in this case, microbubbles, can be estimated to a precision far higher than the diffraction-limited resolution of the system. Here, this is exploited to accumulate the localization of flowing microbubbles over thousands of images to re-create a super-resolved image of vascular structures.
Increased worldwide attention within this area of research means that the contributions of many international groups are facilitating rapid progression, diversity and innovation in this field, which will inevitably encourage and accelerate clinical implementation. Despite methodological differences, there are a number of common steps that form the basis of the technique throughout the literature. These are (note: post-processing steps are visualized in Fig. 1):

Microbubble introduction
All current single-bubble localization methods for super-resolution involve the injection of contrast agents as a bolus or infusion. For localization-based methods, the concentration of the agent is required to be low enough that bubbles can be spatially separated by the image system's diffraction-limited point spread function after post-treatment filtering. Recent sparsity-based approaches (Bar-Zion et al. 2018a) and deep learningbased methods (van Sloun et al. 2018a(van Sloun et al. , 2018b can alleviate these requirements and permit higher concentrations

Video acquisition
An ultrasound pulse is emitted into a medium containing microbubbles. Then, a video of microbubble flow is acquired (Fig. 1a). This could be B-mode with or without contrast-specific pulse sequences, and at conventional or ultrafast frame rates. The received data can be collected as a matrix of radiofrequency (RF) data acquired by each channel or as beamformed image data. Differences may depend on equipment availability, restrictions with data accessibility or specific clinical protocols.

Motion correction
Long video acquisitions are often required to observe the smallest vessels. As super-resolved images are created from the superposition of many localizations gathered over time, motion between frames will significantly affect their visualization. Crucially, the scale of motion present in clinical scanning is often orders of magnitude larger than the super-resolved vessels themselves. Adequate motion correction techniques are therefore vital.

Microbubble detection
A crucial part of the super-resolution process is distinguishing microbubbles from the surrounding tissue (Fig. 1b). This step creates candidate bubble regions for subsequent localization. Inadequate bubble detection will mean that long acquisition times are needed to obtain sufficient localizations for the final rendering. Too many erroneous signals may create a challenge for subsequent filtering processes and result in noise in the super-resolved image. Techniques to extract microbubble signals from tissue have been thoroughly studied in the field of CEUS (Cosgrove and Lassau 2010), while new techniques are also being introduced and developed within the super-resolution ultrasound field itself.

Microbubble isolation
Once detected, a filtering step can identify isolated microbubbles in each frame (Fig. 1c). Indeed, the echoes from two microbubbles that are closer than a few wavelengths in an image will interfere, making the corresponding localizations of each microbubble inaccurate. These signals should, therefore, be rejected at this stage or at the tracking stage. Consequently, only a limited number of microbubbles can be detected in each image to avoid such overlapping.

Localization
A key step in the process is the localization of these isolated signals (Fig. 1d). Ultrasound waves are coherent when propagating into tissue. Furthermore, the response of ultrasound waves to a single isolated point scatterer is generally well defined by the point spread function of the imaging system. Under the assumption that these signals are isolated, this procedure can estimate the location of the underlying scatterers to a precision far higher than the diffraction-limited resolution of the system. Localization of the microbubbles can be performed with a number of techniques in the RF domain or on beamformed images. In reality, the ultimate limit on the resolution in a super-resolved image is determined by the precision of this localization step.

Tracking
By determining the displacement of a microbubble between two images, a velocity vector can be created (Fig. 1e). With sufficient detections, features of the local vascular velocity can be determined through the creation of detailed velocity maps. The tracking of these localized points to define the paths and velocity of microbubbles in microvessels can vastly improve the quality and information content of images and modify their interpretation. Such microbubble tracking has distinct advantages over Doppler in spatial resolution and is largely independent of the direction of the flow displacements.

Visualization
The visualization of localizations, their density, or their calculated velocities is performed by creating a map from the microbubble positions after correcting for motion to reveal details of the vascular network (Fig. 1f). This can be performed either by accumulating bubble detections in fine grids with pixel sizes that reflect the scale of the super-resolved detail or by plotting each localization as a Gaussian distribution, where its size accounts for the uncertainty in the localization procedure.

DEFINITIONS
Super-resolution ultrasound imaging is a new technique with much promise, but with terminology that could become confusing. One of the goals of this review is to establish a common language and provide several definitions. This can be particularly important to unite applications that can differ widely, such as the simple improvement of vascular imaging in deep-seated tissue in humans or the direct observation of capillaries in the animal brain.

Resolution
In super-resolution ultrasound imaging, as in any imaging modality, resolution is the capacity to separate two objects to be imaged. However, in our case, these objects are the vessels to be imaged. The separation of two vessels or the capacity to separate two distinct flows patterns within the same vessel is the ultimate measurement of resolution in our field. Other approaches from the optical field should also be exploited such as the Fourier shell correlation threshold (Van Heel et al. 2005).

Super-resolution
Super-resolution is the capacity to distinguish two objects, here vessels, beyond the classic limit. Several classic limits have been proposed, and their relevance depends ultimately on the application. The lower limit for resolution is usually considered to be the diffraction barrier at half the wavelength. One can also use the Rayleigh resolution criterion of 1.22λ (focal length/aperture) (Cobbold 2006). The latter definition is more lenient in the vast majority of cases, but providing a 150-mm resolution at a 15-cm depth with a 5-cm-aperture, 5-MHz transducer remains an interesting exploit for medical applications, especially for larger organs. In both cases, the particular limit should be clearly stated in the article to allow a better comparison between studies.
It should be noted that post-filtering techniques on diffraction-limited images are not considered super-resolution as they assume a priori knowledge of the vessels to be separated. Indeed, prior information on the ultrasound echo sources has to be exploited before their accumulation, as it is their physical interference in each frame that has to be avoided to achieve super-resolution.

Super-resolution ultrasound imaging
Super-resolution ultrasound imaging (SRUS) refers to a field of study comprising several techniques-using high-frequency acoustic waves-that can distinguish objects or structures closer than the diffraction limit, which is about a half-wavelength. These techniques include ultrasound localization microscopy, super-resolution fluctuation imaging or structured illumination.

Ultrasound localization microscopy
Ultrasound localization microscopy (ULM) refers to a super-resolution ultrasound technique. It exploits the accumulation of sub-wavelength localizations of many separate sources over a great number of frames to reconstruct a super-resolved composite image. In each frame, the localization sources are sufficiently sparse so that their interference can be disregarded or taken into account to yield a precise positioning.

Localization precision
Ultrasound localization microscopy relies on the accumulation of the centroid's localization of isolated microbubbles to reconstruct a super-resolved image of the microvasculature in which they flow. One of the key aspects is hence the capacity to localize a single source with a good precision either on the radiofrequency channel data or on the beamformed image. The localization precision is the closeness between several measures of the localized position of a single isolated microbubble. It depends, for instance, on the signal-to-noise ratio (SNR) . This precision is often close to the size of the microbubbles and that of the capillaries. This precision is not the resolution of the imaging system as the microbubbles are not the structure to be imaged, but they are rather the probes that highlight that structure.

Isolated microbubbles
A source is deemed isolated if, in a single image, no other sources can bias significantly the localization of its center. This depends on the techniques that have been developed to distinguish non-isolated sources, along with the contrast-to-tissue ratio (CTR), the SNR and the PSF of the system.

Temporal resolution
Temporal resolution corresponds to the total time to display a super-resolved image. It has to be distinguished from the frame rate of the acquisition. Temporal resolution and spatial resolution are closely linked because the accumulation of microbubble localization will yield an improved image. After injection, microbubbles remain rare compared with red blood cells (about 1/100,000) and it takes time for these agents to fill every single capillary.

MICROBUBBLE DIFFERENTIATION VERSUS TISSUE
The first stage for ULM is the separation of microbubbles from tissue. Fortunately, because the first ultrasonic detection of microbubbles in vivo (Gramiak and Shah 1968), several techniques have been proposed to highlight the signal of contrast agents. With respect to tissue or blood, microbubbles have several distinct characteristics.
First, these intravascular microbubbles move with the blood, allowing them to be separated from static or slowly moving tissue. The first microbubble detection technique was hence ultrasound Doppler. For instance, on continuous Doppler, isolated microbubbles in the blood make a very distinctive broadband "ploc" as they pass through the acoustic beam.
As discussed before, microbubbles are clearly distinct from the blood by the intensity of their echo. Their scattering cross-section can be three orders of magnitude greater than their physical size (Wheatley et al. 1990). This is partly owing to the important acoustic mismatch between air and an aqueous solution. But it is also owing to a very beneficial coincidence of nature. Indeed, microbubbles act as a resonator, where the compressible air is the spring and the surrounding water is the mass. The resonant frequency of microbubbles small enough to move in the capillary falls exactly in the megahertz range within which clinical ultrasound imaging is performed. This resonance increases the echo of microbubbles by several orders of magnitude.
However, beyond moving faster than tissue and scattering more than red blood cells, microbubbles can also be detected through their non-linear echoes (Schrope et al. 1992;Shi and Forsberg 2000). Because of their strong oscillations around their resonance, the acoustic behavior of microbubbles goes beyond simply following the incident wave. In the scattered spectrum, this non-linearity can be detected as a rich harmonic, subharmonic, superharmonic and ultraharmonic content and microbubbles can be detected as such (Burns et al. 1994;Forsberg et al. 2000;Basude and Wheatley, 2001). However, multipulse sequences are even more sensitive to contrast agents, and commercial systems have implemented techniques such as pulse inversion (PI) (Simpson et al. 1999;Eckersley et al. 2005), amplitude modulation or a combination to highlight perfusion after microbubble injection.
Finally, beyond a few hundred kilopascals, microbubbles oscillate so violently that they can disrupt, their shells rupturing and their gas diffusing into the surrounding fluid (Postema et al. 2004). The disappearance of the microbubbles after a short but intense pulse (Chomas et al. 2001) can be used as a contrast mechanism through differential imaging (DI), but can also be applied for flash-replenishment techniques (Johnson and Powers 1995;Greis 2009).
All these approaches have already been exploited to detect microbubbles for ultrasound localization microscopy. Specific techniques depend generally on the frequency used, which itself depends on the depth at which the imaging is performed. At higher frequencies, microbubbles are poorly resonant and scatter little harmonics, so techniques based on the motion or disruption of microbubbles are preferable. For instance, spatiotemporal filtering with singular value decomposition (SVD) has, importantly, increased Doppler sensitivity and has been exploited to separate microbubbles from tissue (Errico et al. 2015;Desailly et al. 2017;Song et al. 2017). In addition, deep learning techniques can be used to solve the microbubble tissue separation problem through an unfolded robust principal component scheme (Cohen et al. 2019). More simply, DI can also highlight the motion or disruption of microbubbles (Desailly et al. 2013).
Closer to resonance, it is highly beneficial to use non-linear techniques. Pulse inversion and amplitude modulation have been used to extract the signal from microbubbles (Viessman et al. 2013;Christensen-Jeffries et al. 2015, 2017a, 2017b. These various approaches were compared in a recent article by Brown et al. (2019).
For cardiac, abdominal or transcranial applications, depths on the order of 10 cm are required. However, detecting separable microbubbles in microvessels at these depths is challenging because the intensity of the signal backscattered from microbubbles is low, and increasing the transmit energy is not a viable option because of microbubble destruction. Previous studies have proposed the use of larger microbubbles to increase the amplitude of the backscattered signals (Lin et al. 2017b); however further improvements are needed to achieve super-resolution at clinically relevant depths. Furthermore, the accurate evaluation of microvascular structure requires imaging in three dimensions where the small aperture sizes of existing matrix arrays are challenged in terms of sensitivity and contrast. However, the volume fraction of tissue that is occupied by vessels is typically less than 10%; therefore, imaging sequences do not need to interrogate the entire volume.
A simultaneous multifocus beamforming approach was proposed to simultaneously sonicate two or more foci with a single emission (Espindola et al. 2018). This has the advantage of retaining a high frame rate, yet achieving improved sensitivity to microbubbles. In the limit of one target, this beam reduces to a conventional focused transmission; and for an infinite number of targets, it converges to plane wave imaging. By interleaving targeting sequences with imaging sequences, only volumes of tissue that contain visible bubbles can be targeted, thus reducing the number of transmitÀreceive events required to construct an image while improving the sensitivity. Experimental results have indicated that the adaptive multifocus sequence successfully detects 744 microbubble events at 60 mm, while these are undetectable by the plane wave sequence under the same imaging conditions. At a shallower depth of 44 mm, the adaptive multifocus method detected 6.9 times more bubbles than plane wave imaging (1763 vs. 257 bubble events).

MICROBUBBLE ISOLATION
A current research area is the method by which the microbubble signals are separated. Indeed, to make localization work precisely, scattered signals from microbubbles need to be isolated in each individual image. Neighboring echoes that are closer than the Rayleigh criterion in a single frame cannot be distinguished from one another and interference will make the corresponding localizations inaccurate. Notably, in van Sloun et al. (2019), a deep neural network was specifically trained to learn the interference patterns of closely spaced microbubbles and directly infer their locations. Yet, in most methods, these signals are typically rejected at this stage.
Super-resolution ultrasound has been illustrated using relatively high microbubble concentrations in the blood, where signal separability is achieved by significantly limiting the number of microbubbles detected in each frame. As early as 2011, Siepman et al. and Couture et al. achieved microbubble isolation after bolus injection using pairwise frame subtraction, where preceding frames were subtracted to detect only the signals from moving or disrupted microbubbles in each image. In Siepman et al., this was performed using a high frequency of 50 MHz, meaning a considerably smaller diffraction-limited signal size, as well as the inherent compromise of limited penetration depth. With a lower frequency suitable for deeper imaging, 872 Ultrasound in Medicine & Biology Volume 46, Number 4, 2020 Couture et al. (2011) and later Desailly et al. (2013) implemented pairwise frame subtraction, termed differential imaging (DI), using ultrafast imaging frame rates (up to 20,000 Hz for shallow tissue). This allowed microbubble isolation by detecting the decorrelation between successive ultrasound echoes caused by the fast movement or disruption of bubbles. Although applicable, the disruption of microbubbles remains undesirable as it restricts microbubble tracking capabilities, as well as the visualization of microbubbles flowing through small microvessels. The most direct way to achieve isolated signals is to reduce the concentration of microbubbles in the blood. By use of an infusion, a constant microbubble concentration can be maintained that minimizes the likelihood of overlapping signals occurring. This method was implemented successfully in some of the earliest demonstrations of ultrasound super-resolution, including those of Viessman et al. (2013), O'Reilly and Hynynen (2013) and Christensen-Jeffries et al. (2015), and has been used extensively in the field since. Notably, detection and isolation techniques are often intertwined, and some detection procedures may perform both of these steps to some degree. For example, either B-mode, contrast-specific imaging techniques (e.g., PI or contrast pulse sequencing [CPS]) or SVD (Desailly et al. 2017) has been used primarily for detection; however, these also help to isolate signals because they reduce the number of bubbles that are visible in each frame. This can be owing to, for example, some bubbles being off-resonance or weakly responding in the case of contrast-specific methods, or slow-moving or stationary bubbles in the case of SVD.
Despite the advantages of infusions for this technique, standard clinical contrast imaging protocols often involve a bolus injection. As such, soon after injection, the bubble concentration in the blood reaches its peak, and isolated detections are often nearly impossible (unless detection strategies are such that only a small proportion of bubbles are visualized). In these instances, the selection of suitable video segments during both the inflow and later stages of the acquisition can improve the separability between individual microbubbles. This method was implemented in vivo by Ackermann and Schmitz (2016). Alternatively, multiple lower-volume bolus injections can be performed to improve the duration of the desired bubble concentration as implemented in Errico et al. (2015) and Zhu et al. (2019). Additionally, destructionÀreplenishment sequences can be performed after contrast injection to capture the reperfusion of microbubbles into a region before their signals overlap, as in Opacic et al. (2018).
Although these methods aim to extract mainly isolated microbubble signals, an additional step is often required to reject any remaining interfering signals. For example, in O'Reilly and Hynenen (2013), all in vitro frames containing too high a sidelobe intensity (50% global maximum in frame) were rejected. Most other methods have involved searching for or rejecting signals based on their similarity to expected microbubble signatures. In Christensen-Jeffries et al. (2015), the size and intensity of any connected signal regions were compared with those of the expected PSF (point-spread function) based on characterization experiments. In Errico et al. (2015), after spatiotemporal filtering procedures, a Gaussian model based on typical microbubble signal properties was used to perform a 2-D normalized crosscorrelation with each interpolated frame. Similarly, in Ackermann and Schmitz, (2016), after background subtraction, a 2-D convolution of the foreground frame was performed with a Gaussian kernel to identify single microbubble signals for subsequent localization.
These methods can be performed at both conventional clinical imaging frame rates (100 Hz

LOCALIZATION
The prior knowledge that the diffraction-limited image of a microbubble originates from a single source allows estimation of its location with a precision well beyond the diffraction limit. Within RF data, a single echo from a microbubble appears as a hyperbola, determined by the time-of-flight to arrive at each transducer element. In beamformed data, which are often more readily available, particularly in a clinical setting, microbubbles appear as individual "blurs" in an image. Localizations from each domain are similar as they both relate to the fitting of the time-of-flight and, as distinct sources, can both be determined very precisely. Nevertheless, the exact equivalence of the localizations in each domain is yet to be fully examined and may depend on factors such as the respective SNR, the benefit of phase information in the RF space and any non-linear processes used to transfer into the beamformed space.
The localization step requires identification of some feature of the returned signal that corresponds to the bubble's position. In most existing publications regarding ultrasound super-resolution, the position is expected to coincide with a measure related to the "center" of the echo because of the assumption that the image is a convolution of a point source with the system PSF. In early work using RF data, Couture et al. (2011) and Desailly et al. (2013) found a local axial maximum from the traveling hyperboloid in RF data lines by finding the summit of a parabola fit. In some of the earliest studies using post-beamformed data, localization was performed by calculating the intensity-weighted center of mass (Siepmann et al. 2011;Viessmann et al. 2013;Christensen-Jeffries et al. 2015), fitting either a 3-D Gaussian function to the original backscatter signal (O'Reilly and Hynynen 2013) or a 2-D Gaussian after deconvolving with a predicted PSF (Errico et al. 2015), or cross-correlation of signals with an expected response (Ackermann and Schmitz, (2016)).
In optical localization microscopy, provided a sufficient number of emitted photons are sampled, the signal received from a molecule is well defined, and the intensity centroid is generally used to define its position to nanometer scale. Likewise, if each microbubble from the insonated population behaves similarly, generating a signal close to the system PSF, the aforementioned methods adequately represent the relative position of the contrast agents. However, more recently. Christensen-Jeffries et al. (2017b) demonstrated in vitro and in simulations that the considerable variability of bubble responses makes it challenging to predict both their linear and nonlinear responses and that this variability should be taken into account during localization. It indicated that centroiding, peak detection and 2-D Gaussian fitting methods introduced hundreds of micrometers in error when using non-linear imaging at low clinical imaging frequencies. As a result, a new "onset" localization method was proposed to identify the beginning of the signal in the axial direction while ignoring ringing or elongated signal duration. This resulted in considerable improvement in localization error compared with existing methods and is particularly important in the resonant regime of microbubbles.
A recent article by Song et al. (2018) described the influence of the spatial sampling of the beamforming on the localization error. This study found that the Fourier analysis of an oversampled spatial profile of the microbubble signal could guide the choice of beamforming spatial sampling frequency. Furthermore, parametric Gaussian fitting and centroid-based localization on upsampled data had better localization performance and were less susceptible to quantization error than peak intensity-based localization methods. When spatial sampling resolution was low, parametric Gaussian fittingbased localization had the best performance in suppressing quantization error, and could produce acceptable microvessel imaging with no significant quantization artifacts.
The theoretical limits of ultrasound super-resolution, and under what conditions these could be achieved, are of vital importance. Many parameters, including physiologic ones, influence the separability of microvessels or other structures in super-resolution imaging.
However, the best attainable resolution is limited by the localization precision of each microbubble. This lower bound is linked to the minimal time delay, which can be estimated between similar echoes, and is known as the CramerÀRao lower bound (Walker and Trahey 1995). In 2015, Desailly et al. investigated the theoretical resolution limit to which a single, linear, point scatterer can be localized using hyperboloid fitting and validated this experimentally. In the far field, Desailly et al. found that in this situation, the achievable resolution was dependent on the standard deviation of arrival time estimates on each channel RF line. These estimates depend on the pulse bandwidth, the pulse center frequency, the SNR and the temporal jitter between electronic channels of the scanner. Both the number of transducer elements and the speed of sound also influence the resolution, while the lateral resolution was also dependent on the depth of the target and the size of the array aperture.
Alongside the extraction and isolation of individual microbubbles from tissue, the SNR of each microbubble signal is crucial to the performance of the localization algorithm. In Foiret et al. (2017), the coherence factor, the ratio of coherent intensity to incoherent intensity (Mallart and Fink 1994), was applied to CPS images to improve the SNR of echoes from individual microbubbles. Ghosh et al. (2017) and Lin et al. (2017b) have reported that selecting a population of larger microbubbles improves the SNR from individual microbubbles. Finally, because of the challenge of separating bubble signals from noise in images using just spatial information, Song et al. (2017) proposed to denoise the microbubble signal in a spatiotemporal domain using a nonlocal means filter to improve the separation of microbubble "tracks" from random background noise that does not resemble any feature or pattern. More robust denoising such as this is beneficial to improve the achievable isolation and localization precision.
Microbubble localization can also be improved by modifying the beamforming process, from a delay-andsum to a maximum variance approach, for example (Diamantis et al. 2018).

TRACKING
Current ULM techniques usually go beyond the simple accumulation of microbubbles subwavelength positions to create an image. The fact that microbubbles remain exclusively intravascular allows the displacement of the microbubbles to be interpolated from the various localizations along a track. Christensen-Jeffries et al. (2015) and Errico et al. (2015) used such an approach to establish the direction and the velocity of each microbubble at the micrometer scale. An important advantage of tracking is the exclusion of artifactual microbubbles.
Indeed, by removing very short tracks (a few images), only microbubbles with a coherent path are retained, forming an image with lower noise. Early on, tracking algorithms were based on the maximum intensity cross correlation within a search window (Christensen-Jeffries et al. 2015) or closest-neighbor detector, which identified from the localized position of a microbubble on image i the most probable localization on the image i + 1 (Errico et al. 2015). Algorithms have to define several thresholds such as the maximum distance to the next microbubbles for starting and ending a track.
Currently, more advanced algorithms can be exploited such as Monte Carlo based on Markov chain (Ackermann and Schmitz 2016;Opacic et al. 2018) or joint tracking and detection in a Kalman framework that is regularized through optical flow estimates . Hungarian assignment, already exploited in transportation analysis (Kuhn 1955;Munkres 1957), can also be used to improve tracking (Song et al. 2018). New approaches to bypassing the localization process through a dynamic method have also been discussed (Alberti et al. 2018). Tracking was also performed in combination with motion correction, which is discussed in the next section.
The images are often reconstructed by projecting the detected tracks on an image grid where the pixels are much smaller than the wavelength, 5 £ 5 mm, for example. However, it is important to note that the size of the vessel where a single microbubble is detected is misrepresented using this approach. For instance, a 5-mm "vessel" could simply be a larger vessel, such as an arteriole, where a single microbubble has passed. In fact, even the aorta would look like a capillary if only a single track is detected inside its lumen. The true size of a vessel should only be established with a large number of microbubble tracks, at a minimum sufficient to cover the vessel diameter. Consequently, the minimal number of tracks required to reconstruct a vessel was proposed by Hingot et al. (2019) as the width of the vessel to be imaged divided by the super-resolved pixel size.

MOTION CORRECTION
Super-resolution ultrasound imaging can achieve resolution down to 5 mm. However, motion in the body or from the ultrasound transducer can be several orders of magnitude beyond this level, especially because microscopic resolution requires long acquisition times Christensen-Jeffries et al. 2019)sometimes tens of minutes to reconstruct capillaries. A brute force approach is to fix the organ to be imaged with a stereotactic frame for the brain (Errico et al. 2015) . However, this approach is limited to organs not affected by breathing motion, and its implementation is difficult to imagine in humans apart from neuroimaging. Christensen-Jeffries et al. (2015) and Lin et al. (2017a) used respiratory gating to exclude images of tumors that were affected by breathing. Such an approach requires a steady motion of the organ, which has to return exactly to its original position.
Several techniques are now being implemented to correct the motion based on the ultrasound images themselves. The phase correlation of speckle motion can provide an idea of the collective rigid motion of tissue (Hingot et al. 2016). It is, however, necessary to remove the effect of the microbubbles, which can bias the motion measurement. This can be done directly on Bmode as has been reported by Foiret et al. (2017). More advanced non-rigid motion correction can also be implemented to take into account complex motions in the image . Indeed, because of its phase sensitivity, IQ data provide ample information to detect micrometric motion as reported in the fields of speckle tracking (Walker and Trahey 1995) and transient elastography (Sandrin et al. 1999). Motion correction packages already implemented on MATLAB can also be used directly on the IQ data.
However, because it cannot be tracked on the ultrasound data, out-of-plane motion cannot be easily corrected. Hence, motion correction is one of the principal arguments in favor of 3-D ultrasound localization microscopy, as displacements in all directions could be accounted for.

Three dimensions
Super-resolution performed on a 2-D plane will always be limited for several reasons. First, tissue motion is inherently 3-D, and its correction necessitates sampling in all directions, including the elevation axis. Second, vessels follow tortuous paths in three dimensions, and super-resolved imaging of small sections of numerous vessels intercepted by the imaging plane would give a restricted amount of information. Third, microbubble positions need to be accumulated for a long time, up to several minutes, to attain the capillary level. Reconstructing a volume through plane-by-plane imaging as in Errico et al. (2015), Lin et al. (2017b) or Zhu et al. (2019) is difficult to transfer successfully to the clinic. Indeed, it could make ultrasound imaging even more user dependent, as the plane will have to be chosen and preserved for an extended time to avoid trading off the underlying resolution.
Because of these 2-D imaging drawbacks, superresolution ultrasound imaging was rapidly implemented with matrix (Desailly et al. 2013) or hemispherical (O'Reilly and Hynynen 2013) probes to perform ultrasound localization over volumes. Both these groups were able to map, in three dimensions, microbubbles flowing in microvessels with subwavelength resolution. However, these approaches were limited either in time resolution or in the reduced imaging volume. The latter issue also limited another approach exploiting two perpendicular probes for 3-D reconstruction (Christensen-Jeffries et al. 2017a, 2017b. Appropriate 3-D imaging can be performed with matrix arrays, where each transducer can be addressed either through a beamformer in the probe or by connecting each channel to the ultrasound system. However, matrix arrays are costly and present some technological hurdles at higher frequencies, making the technology less available. Heiles et al. (2019) described in vitro a 3-D super-resolution imaging device in which each of the 1024 channels (32 £ 32) was connected to a programmable ultrafast scanner (Provost et al. 2014). In vivo, this apparatus was able to reconstruct entire super-resolved volumes of the rat brain, alleviating the previously described drawbacks with 3-D tracking and motion correction. Harput et al. (2019b) reduced the number of connected channels by implementing super-resolution on a 2-D sparse array. 3-D super-resolution could even be performed with a clinical scanner as described by Christensen-Jeffries et al. (2018).

TRADE-OFFS AND LIMITATIONS
With the potential for spatial resolution to be improved by orders of magnitude in ultrasound, the performance of the technique now depends on many new factors, including vascular flow rates, contrast-signal density, time resolution, SNR, probe geometry, postprocessing algorithms and many more parameters. Although each of the methodological steps has been investigated, adapted and improved by various groups, several new trade-offs have been discovered, and numerous strategies have been introduced to address them.
First, much discussion has arisen in recent years regarding the requirement for sufficient measurement times to adequately visualize vascular structures. For all the aforementioned localization methods, microbubble signals need to be spatially isolated for accurate localization. Thus, there is an upper limit on the number of isolated signals attainable in each frame. The temporal domain is, therefore, used to gather sufficient localizations to adequately sample the structure of interest. Yet, because microbubbles travel within the bloodstream, slow flow rates associated with the smallest vessels limit our ability to sample these structures over short periods.  Dencks et al. (2017) provided an expression for the acquisition time needed for the localization coverage in an image to saturate to a value assumed to be proportional to relative blood volume (rBV). The reliability of rBV estimates from shortened measurement times is then examined experimentally for bolus scans of mouse tumors at high frequency and conventional frame rates (40 MHz, 50 Hz), where approximately 45À170 s was found to be necessary for a 90% mapping of the microvessels. In Dencks et al. (2017) the saturation model was applied to clinical data on breast tumors (10 MHz, 15 Hz), and 28%À50% filling was found for a 90-s acquisition time. Hingot et al. (2019) aimed to establish the link between imaging time and microvascular scale in super-resolved images and found that larger vessels (>100 mm) are fully reconstructed within 10 s, while mapping the entire capillary network would take tens of minutes. This is owing to the limited blood flow in capillaries, along with their density in tissue. Moreover, the characterization of flow profiles was shown to be bound by even more stringent temporal limitations. Christensen-Jeffries et al. (2019) developed a generalized model to predict the optimal bubble signal density and acquisition time for user-defined imaging conditions, parameters and target vasculature. Similarly, this study estimated acquisition times longer than 10 min for microvascular targets at a 5-cm depth and λ/20 super-resolution. Both studies highlight that super-resolved vascular imaging is inherently dependent on the physiology of the vasculature. Nevertheless, the distinct advantages provided by ultrasound super-resolution, along with the strengths of this technique over comparable imaging modalities, mean that several minutes of acquisition time may be acceptable in the clinical world.
Super-resolution ultrasound techniques have a particular advantage over high-frequency techniques to improve resolution without the conventional penetration trade-off associated with higher frequencies. Nevertheless, increased depth still has the capability to degrade the performance of ultrasound super-resolution. First, imaging at greater depths will increase the minimum super-resolution acquisition time. This is because greater depths require lower transmit frequencies as in conventional imaging. For instance, super-resolution ultrasound imaging at greater depths (>10 cm) or through the skull would require 1-to 3-MHz transducers, and at depths of a few centimeters, this would be achieved around 4À12 MHz, whereas small animal imaging would be performed at frequencies beyond 15 MHz. Each time, super-resolution imaging would have an enhanced resolution with respect to conventional imaging, but it would still be degraded as the frequency is lowered to penetrate deeper. The resulting larger PSF size decreases the maximum number of separable signals that can be detected in each frame. Therefore, the acquisition time required to create a super-resolved image increases with depth, as discussed in Christensen-Jeffries et al. (2019). This, along with the inherent time-of-flight versus depth relationship, therefore represents a trade-off in minimum possible super-resolution acquisition time with depth.
Another trade-off is with depth is reduced sensitivity, particularly using plane-wave imaging. Conventional focused imaging has a lower frame rate, but because of the well-established contrast and SNR, it has a higher sensitivity to microbubble detection than plane waves. For high-frame-rate sonications, the wide planar beam used to illuminate the field of view produces poor contrast, and the transmitted energy is limited to minimize microbubble destruction over long imaging times. This means detecting microbubbles in microvessels at depth can be challenging. Strategies to improve the sensitivity of contrast agents at depth are being investigated for super-resolution imaging, for example, the use of larger microbubbles (Ghosh et al. 2017;Lin et al. 2017b) and the use of adaptive multifocus methods, which combine the high frame rates of plane wave imaging with the increased bubble detection sensitivity of focused beams (Espindola et al. 2018).
More generally, the performance of post-processing techniques and their influence on factors such as SNR and localization error present potential trade-offs to be considered when making methodological choices. So far, Brown et al. (2019) have studied the performance of various microbubble detection methods using simulations, namely, PI, DI and SVD filtering. They found that PI is the most appropriate for some applications, but also the most dependent on transducer bandwidth. PI was found to be largely independent of flow direction and speed compared with SVD and DI and, thus, more appropriate for visualizing the slowest flows and tortuous vasculature. SVD was found to be preferable for highfrequency acquisitions, where localization precision on the order of a few microns is possible, but can introduce localization errors on the order of hundreds of microns over the speed range 0À2 mm/s, while DI is suitable only for flow rates >0.5 mm/s or as flow becomes more axial. Such investigations highlight the need for an increased understanding of how different techniques can affect the super-resolution process.
Because of the nature of research progression in this field, a number of other issues requires increased attention in the near future. There is an inherent compromise between the SNR of the signal and the localization precision. This has been highlighted in studies such as Desailly et al. (2015), Ghosh et al. (2017), Lin et al. (2017a) and Espindola et al. (2018). Exploration of the impact of SNR, specifically in relation to depth and aberration effects, is greatly needed. Additionally, with increased acquisition time, and with the field moving increasingly into the 3-D domain, the management and analysis of large data will become more and more important in coming years. However, it must be stated that ultrasound localization is a powerful data-compression technique. Indeed, radiofrequency channel dataterabits in size for long 3-D acquisitions-can be converted to an accumulation of hundreds of thousands of small vectors-a few bits each-representing the superresolved path of microbubbles.
Another important trade-off in ultrasound localization microscopy-based methods is the concentration of the injected microbubbles. Too low a concentration would certainly further prolong the acquisition time, but increasing the concentration would increase the probability of having more than one microbubbles in the same imaging resolution cell, and too high a concentration would cause many detected signals to be not useful for super-localization. Such multiple bubble events, if not rejected, would contribute to the noise of the super-resolution image. Christensen-Jeffries et al. (2019) has reported that an optimal bubble concentration exists for super-resolution ultrasound imaging.
Super-resolution ultrasound imaging needs to demonstrate that it can provide superior microvascular and microstructural information for clinical applications to be relevant. In particular, it should be compared with the clinical diagnostic benefit of perfusion ultrasound imaging. Is there a benefit to displaying super-resolved microbubble tracks compared with intensity images or maximum intensity projection? For instance, if we reconstruct a single microbubble track, it is impossible to determine if the microbubble was moving inside a large or a small vessel. Several microbubbles are hence required in each vessel to determine its true size and bifurcation. To assess its accuracy for larger vessels, for mapping subwavelength flow patterns, super-resolution ultrasound imaging should be compared with other vascular imaging modalities such as ultrafast Doppler (Bercoff et al. 2011) or optical techniques such as microscopy. The field of ultrasound should explore these concepts of accuracy more deeply and move beyond the concerns over the smallest observable vessel or the distinction between two vessels.

SUPER-RESOLUTION ULTRASOUND USING PHASE-SHIFT AGENTS
A fundamental difference between the localizationbased optical super-resolution methods and the microbubble-based ultrasound counterparts is that there is no switching on and off of the contrast agents for the ultrasound approach. In optical super-resolution such as PALM (Betzig et al. 2006) and STORM (Rust et al. 2006), some individual contrast agents (i.e., fluorescence molecules) from a high-concentration population are stochastically switched on and off, creating "blinking" signals based on which these individual agents can be localized. The accumulation of the localizations depends on the stochastic on and off switching of different subpopulations of the agents. On the other hand, microbubbles, once injected, are already in an "on" state. Therefore, for localization-based super-resolution ultrasound approaches to work, low concentration of the bubbles are required to ensure the spatial distribution of the bubbles are sparse, and physiologic flow during imaging is also required to displace the bubbles so that different localizations can be accumulated over space and time.
While localizing the same contrast agent over time in the presence of flow offers the opportunity to generate super-resolved velocity maps of microflow (Errico et al. 2015;Christensen-Jeffries et al. 2015), a limitation of the approach is that they heavily depend on the agent concentration and presence of flow and, consequently, may take much longer acquisition time. Regardless of the imaging frame rate, time has to be spent in waiting for spatially isolated individual bubbles to flow past the microvessels of interest (Dencks et al. 2017;Christensen-Jeffries 2017;Hingot et al. 2019). Furthermore, the in vivo microbubble concentration is difficult to control and it can change quickly over time, which can cause either too few localizations in small microvessels or, in larger vessels with high agent concentrations, too many incorrect localizations of multiple bubbles that, if not discarded, can contribute to image noise.
Phase-change droplets, typically made from perfluorocarbons, can change phase from liquid to gas (i.e., be vaporized into bubbles) by ultrasound (Kripfgans et al. 2000;Shpak et al. 2014). Such phase-change agents have been studied in the past couple of decades for both therapeutic and diagnostic applications (Kripfgans 2000;Kawabata et al. 2005;Rapoport et al. 2009;Reznik et al. 2011;Sheeran et al. 2012) and recently have been reported to offer subwavelength ultrasound drug delivery and super-resolution imaging (Hingot et al. 2016). Most current research in phase-change agents is focused on nanometric agents, which allow longer circulation and tumor deposition through the enhanced permeability and retention effect, while reducing the size of the resulting microbubbles.
As an alternative ultrasound imaging contrast agent, nanodroplets offer three potential advantages over microbubbles. First, the nanodroplets can be selectively activated, both spatially and temporally, to provide an ultrasound contrast signal on demand. Second, the nanosize of nanodroplets potentially allows extravasation, for example, into cancerous tissue because of its leaky vasculature and enhanced permeability and retention effects, offering contrast beyond vasculature. Third, it has been found that nanodroplets can last longer during in vivo circulation than microbubbles (Sheeran et al. 2012;Sheeran et al. 2015). In the past decade it has been well reported that the vaporization of such phase-change nanodroplets into microbubbles can be controlled by acoustic pressure (Sheeran et al. 2012). This phase transition also depends on nanodroplet characteristics (e.g., size distribution and boiling point of the gas) and external factors (e.g., microvessel confinement ] and attachment to cells [Zhang et al. 2019b]).
The capability of phase change nanodroplets to be switched on and off by ultrasound on demand allows super-resolution ultrasound imaging to mimic STORM or PALM in optical super-resolution. The switching-off mechanism can be either the recondensation of the resulting microbubble back into a nano-agent or its disappearance through disruption or dissolution. This has recently been observed by Zhang et al. (2018) through acoustic wave sparsely activated localization microscopy (AWSALM) in two proof-of-concept studies (Zhang et al. , 2019a. One of the fundamental distinctions between bubble-and droplet-based super-resolution imaging techniques is that, with bubble-based super-resolution approaches, accumulated signals are limited by low bubble concentration and depend on flow decorrelation, which is slow in the smaller vessels. However, with the droplets being activated and deactivated at the ultrasound pulse repetition frequency, signals can be accumulated as quickly as the imaging pulses are sent, being able to achieve much faster super-resolution imaging. In the droplet-based super-resolution ultrasound approach, much higher concentrations of droplets can be used to pre-populate the vessels (as droplets are invisible to ultrasound before activation) and localizations can then be accumulated quickly by switching them on and off even with very slow or no flow.
In such approaches there are two key factors: droplet composition and the transmitted ultrasound wave parameters. The transmitted ultrasound amplitude has to match the vaporization threshold of the droplets to activate them, and the threshold depends on the boiling point of the perflurocarbon gas, their size distribution and whether the droplets are encapsulated by a shell. In its initial realization, Zhang et al. (2018) used shelled decafluorobutane (boiling point: À1.7˚C) nanodroplets and focused waves with relatively high ultrasound amplitude (mechanical index [MI] = 1.3) were swept through the imaging area to switch on new droplets, while switching off previously activated droplets at the same time, followed by loweramplitude (MI = 0.25) imaging plane waves. Such an approach of separate activating and imaging pulses is necessary in this case as the activation of the decafluorobutane droplets requires acoustic pressure that is difficult for plane-wave to achieve. Consequently, significant time is spent on sweeping the focused activating pulses and this slows down the imaging.
More recently the same authors have developed a fast version of the AWSALM (Zhang 2019a, Fig. 3), in which shelled octafluoropropane (boiling point: À37˚C) nanodroplets are used, which can be activated by ultrasound with a relatively low MI that plane wave transmission can achieve. The study has demonstrated that octafluoruopropane nanodroplets can be simultaneously imaged, activated and deactivated at a very high frame rate (5000 fps) at an intermediate acoustic amplitude (MI = 0.76) to perform super-resolution ultrasound imaging. Two tubes, which contain high-concentration droplets without flow and are separated by a distance below the ultrasound diffraction limit, are well resolved within 200 ms, achieving a two-orders-of-magnitude improvement in temporal resolution when compared with existing ULM techniques.
Photoacoustic super-resolution has been developed by taking advantage of the non-linear photoacoustic signal dependence of certain molecules or nanoparticles on the laser energy (Danielli et al. 2014;Nedosekin et al. 2014). Instead of using ultrasound, perforocarbon nanodroplets can also undergo phase change under laser excitation, allowing photoacoustic super-resolution. Highboiling-point phase-change nanodroplets for super-resolution photoacoustic imaging (Luke et al. 2016) have been used. A unique feature of this approach is that such droplets can condense back to droplets after being vaporized into bubbles under relatively high intensity laser pulses, creating "blinking" signals for ultrasound localization. An advantage of this approach is that the same droplets can be used repeatedly. Flowing optically absorbing particles have also been used for photoacoustic super-resolution (Vilov et al. 2017;Dean-Ben and Razansky 2018). Compared with acoustic activation, such photoacoustic techniques use lasers, and the penetration depth is restricted to regions that can be illuminated optically (typically up to a few centimeters). Furthermore, as the pulse repetition frequency of photoacoustic lasers is usually low (10 Hz), longer acquisition times are required.
Furthermore, because of their size, nanodroplets may be able to leak out of the vasculature, and it is possible to offer in vivo super-resolution beyond the vascular space (e.g., to visualize greater details of tumor tissue structure by using targeted nanodroplets).

SUPER-RESOLUTION ULTRASOUND WITHOUT LOCALIZATION
Standard localization-based techniques for ULM rely on sufficient separability of microbubbles, placing an upper bound on the allowable concentration. Depending on the desired image fidelity, this in turn places a substantial lower bound on the acquisition time, typically being on the order of minutes Dencks et al. 2018Dencks et al. , 2019Hingot et al. 2019). By introducing other trade-offs, super-resolution ultrasound can be performed faster through different strategies exploiting structured illumination (Ilovitsh et al. 2018) or the signal structure.
In this section, several non-localization-based techniques are described that push this lower bound down significantly by permitting higher concentrations without compromising precision. These methods make effective use of structural signal priors to overcome the limitations of localization-based techniques for scenarios with overlapping PSFs. Such priors are obtained explicitly through signal models or implicitly by learning from data. An overview of these approaches is given in Figure 4.

EXPLOITING SIGNAL STRUCTURE
We first assume that the tissue clutter has been removed. This can be achieved by using high-pass filtering, singular value decomposition or deep learning methods (Cohen et al. 2019). The image obtained can then be written as a convolution between the PSF of the system and the microbubbles. We can express this relation in matrix-vector form by vectorizing the low-resolution image frame into a vector y. As reported in Bar-Zion et al. (2017, 2018b the resulting vector follows the measurement model where A is the PSF matrix consisting of shifts of the PSF which relates the vectorized microbubble distribution on a high-resolution grid x to the image, and n is a noise vector. The goal is to recover the high-resolution vectorized image x. In its most basic form, signal structure can be exploited by realizing that the microbubble distribution within a frame x is highly sparse on such a high-resolution grid (van Sloun et al. 2017). This implies that the vector x contains only a few non-zero values and can therefore be recovered by relying on the compressedsensing literature (Eldar and Kutyniok 2012), which provides methods for recovering sparse vectors from underdetermined linear systems such as that in eqn (1) using computationally efficient methods. In particular, one can formulate a convex minimization problem that balances a data fidelity term and a structural sparsity-promoting term, which results in sparse solutions to the problem of eqn (1) on a high-resolution grid.
To allow for increased overlap between microbubbles, one can extend beyond single frames by processing multiple frames of the form (1), leading to what is known as a multiple measurement vector model (MMV) (Eldar and Mishali 2009). Such models permit exploitation of signal structure not only within a frame but also across frames, by relying on the prior that the non-zero elements in each of the frames are located in the same positions. In line with this, Bar-Zion et al. (2017) proposed computing high-order statistics of the temporal microbubble signals at each pixel, achieving high-resolution imaging in a fashion that is similar to super-resolution optical fluctuation imaging (SOFI) in fluorescence microscopy (Dertinger et al. 2009). This ultrasound variant of SOFI makes effective use of the statistical independence between the fluctuations of CEUS signals originating from different vessels. The attainable resolution gain was found to scale with the order of the statistics. In practice, the authors propose to use second-order moments, as estimation of high-order moments requires an exponentially increasing number of image frames to retain the same SNR level.
Sparsity-based super-resolution ultrasound hemodynamic imaging (SUSHI) extends this SOFI-inspired solution by also exploiting the sparsity of the underlying vascular architecture, thereby further improving spatial resolution (Bar-Zion et al. 2018a). This prior can be incorporated by solving a minimization problem in the correlation domain with a sparsity-promoting regularization term on the high-resolution correlation image, thus accounting for both the joint sparsity and the statistical properties of the image. Other than exploiting the sparsity of the vasculature itself, sparse representations may also be used in a transformed domain, for example, through a wavelet or total variation transformation. SUSHI was found to effectively exploit the distinct temporal fluctuations of microbubbles in different vessels when imaging at a sufficiently high frame rate, that is, when leveraging ultrafast imaging. In this scenario, spatial resolutions beyond the diffraction limit were attained at an unprecedented temporal resolution of 25 Hz.
When frame rates are low compared with the vascular flow velocities, the differences between temporal fluctuations among vessels and their high-order statistics exploited in CEUS-SOFI and SUSHI are less significant. To address this, another method that also resides within the MMV framework, simultaneous sparsity-based super-resolution and tracking (triple-SAT), has recently been proposed . Triple-SAT combines frame-by-frame weighted sparse recovery for detection and simultaneous tracking of the microbubbles, thereby explicitly exploiting flow as a prior in the detection process. In doing so it leverages both a structural spatial sparsity prior and a flow-based before resolve vascular architectures with super-resolution. Although the tracking and data association methods used in triple-SAT share similarities with those presented in Fig. 4. Overview of several non-localization methods that leverage priors through signal structure or by learning from data. An input contrast-enhanced ultrasound (CEUS) sequence (a) can be modeled to derive estimators that exploit signal structure (b) or to generate realistic data to train data-driven estimators in the form of deep neural networks (c). MMV = multiple measurement vector. Ackermann and Schmitz, (2016), the method distinguishes itself by the incorporation of the aforementioned tracking-based flow prior in the detection stage, in conjunction with crude optical-flow-based velocity estimates within a Kalman framework. From an optimization point of view, the methods presented above can all be cast in a similar form, as solving regularized least-squares problems of the form where the first term represents a data fidelity term, and R (x) denotes a regularization term that embeds the structural priors. The fidelity term may represent the image directly, as in eqn (1), or may represent the image correlation function, as in SUSHI and SOFI. In those settings, x effectively represents the variance of the pixels in the underlying high-resolution image. By taking into account the structure of A, problem (2) may be solved in the Fourier domain, leading to computationally efficient methods ).
In the case of frame-by-frame sparse recovery, one may choose RðxÞ ¼ k x k 1 , promoting sparse solutions. For triple-SAT, RðxÞ ¼ k diagðwÞ ¢ x k 1 , with w being a flow-prior-based weighting vector, promoting solutions that are both sparse and close to the tracking predictions. In SUSHI, RðxÞ ¼ k Cx k 1 , where C denotes a basis transformation that maps to a domain in which the vascular architecture is sparse (e.g., through the wavelet transform), and x represents the image variance. In van Sloun et al. (2017) and Bar-Zion et al. (2018b), eqn (2) was solved through an efficient Fourier-domain implementation of the fast iterative shrinkage and thresholding algorithm (FISTA).

LEVERAGING DEEP LEARNING
Exploiting structural signal priors in some domain, either at a single-frame level or in an MMV model, has provided an effective set of tools that can alleviate the harsh constraints on allowable microbubble concentration. However, the potential it carries for reducing the acquisition time comes at the cost of high reconstruction time; even highly optimized Fourier-domain implementations of FISTA rely on a time-consuming iterative procedure. Appropriate optimization settings and the optimal regularization (thresholding) parameter λ moreover depend on the actual sparsity of the data, which in practice varies within and across acquisitions. Furthermore, exact knowledge of the PSF is needed to form the measurement matrix.
The aforementioned challenges motivate the development of a solution that is fast, has fixed complexity and is robust across a wide variety of imaging conditions. To that end, van Sloun et al. (2018avan Sloun et al. ( , 2018bvan Sloun et al. ( , 2019 propose deep-ULM, a method that leverages recent advances in deep learning to fulfill these desirable properties. Deep-ULM uses a deep fully convolutional encoderÀdecoder neural network to map lowresolution input frames to sparse localizations on a high-resolution grid. It is trained using simulations that reflect a wide variety of imaging conditions, covering varying microbubble concentrations, backscatter amplitudes, background clutter and noise and variations in the PSF. Once trained, inference by Deep-ULM is very fast, recovering more than a thousand high-resolution patches of 128 £ 128 pixels per second using GPU acceleration. As such, it is about four orders of magnitude faster than the iterative FISTA scheme. Beyond speedup, Deep-ULM also yields improved performance (in particular for high concentrations), which may be attributed to its ability to learn the relation between specific interference patterns of ultrasound waves reflecting off closely spaced microbubbles and their locations .

DEVELOPING APPLICATIONS
The previous section discussed mostly technical aspects of super-resolution techniques. We focused on ULM as it is currently the most common approach for ultrasound super-resolution. However, the choice of the appropriate technique depends mostly on the requested spatial resolution and temporal resolution. ULM still achieves the highest spatial resolution, but an acquisition with sufficient microbubble tracks takes an extended time. Fluctuation-based approaches or structured illumination could provide an interesting compromise between resolution and frame-rate-a middle way between ultrafast Doppler and ULM. Acoustic switching of phasechange contrast agents has shown initial promise to achieve real-time temporal resolution without any significant compromise to spatial resolution in the localization-based approach but further evaluation in vivo is required.
We now discuss two main fields of applications that could benefit strongly from super-resolution ultrasound imaging: oncology and neurology.

ONCOLOGY
Medical imaging is widely recognized as one of the most significant advances in the war on cancer. Early detection of lesions through screening, or diagnosis of suspicious lesions found by other screening methods, while tumors are still small and before metastasis, drastically increases the chances of successful treatment (World Health Organization 2019) Imaging during treatment provides crucial information regarding efficacy of Super-resolution Ultrasound Imaging K. CHRISTENSEN-JEFFRIES et al. the therapeutic approach. Because of its low cost, widespread accessibility and safety, ultrasound has the potential to be a preferred modality for cancer imaging (Kiessling et al. 2014). Yet, traditional gray-scale ultrasound falls behind other modalities for certain cancer imaging applications. Breast imaging is one such application. Breast cancer is the most common cancer in women and second most common cancer overall (World Cancer Research Fund International, 2018). Mammography is the first choice for screening for breast cancer, with a sensitivity typically reported to be between 70% and 90% (albeit reported with varying ranges) (Heijblom et al. 2011). Ultrasound can be combined with mammography to improve this sensitivity, although the specificity of this combination is poor, as low as 40% in women with dense breast tissue (Heijblom et al. 2011). Mammography is even more challenging in women with saline or silicone implants, which drastically attenuate X-rays and substantially limit imaging efficacy. The result is that many imaging results are indeterminate and biopsy is recommended for pathologic confirmation. Unfortunately, the biopsy process is associated with emotional stress, risk of complications and additional time and cost (Rey et al. 2005;Montgomery and McCrone 2010;Montgomery 2010). Improved ultrasound imaging methods, which provide improved specificity and sensitivity with respect to cancer, would provide clear advantages in screening and diagnosis.
These abnormalities provide an opportunity, however, for a diagnostic imaging technique that can visualize microvascular features.
Although traditional gray-scale or Doppler ultrasound techniques are challenged to visualize microvascular features on the scale of early tumor angiogenesis, it has been found that biomarkers of angiogenesis can be detected with the use of high-resolution contrast ultrasound. Recently, investigators have visualized cancer-associated microvascular abnormalities through acoustic angiography, a high-resolution superharmonic imaging technique (Gessner et al. 2012). Acoustic angiography, which requires dual-frequency transducers to excite ultrasound contrast agents near resonance and receive higher-order harmonics, can provide 3-D maps of microvasculature in shallow tissues with a resolution on the order of 150 mm at 2 cm (Gessner et al. 2010). Using this microvascular imaging method, investigators have observed clear indicators of cancer-associated angiogenesis in both xenograft and spontaneous tumors (Fig. 5) (Rao et al. 2015;, and preliminary data suggest that angiogenesis as a biomarker can improve the sensitivity and specificity of ultrasound in cancer substantially over gray-scale imaging. Acoustic angiography is fundamentally limited by the diffraction limit, however, and hence resolution is ultimately limited by the frequency of the imaging system with its inherent trade-off in frequency-dependent depth attenuation, which will prohibit acoustic angiography from ever being a technique that can resolve small vessels involved in angiogenesis in most tissues deeper than a few centimeters. Super-resolution imaging using ultrasound localization microscopy, however, is not bound to diffraction-limited resolution, making this technique a potentially powerful tool to identify malignancies via their angiogenic fingerprint, even in deep tissues. Lin et al. (2017a) applied super-resolution imaging to visualize 3-D microvascular patterns in fibrosarcoma tumors in rats and compared microvascular features with control healthy tissue (Fig. 6). Similar to prior data acquired with acoustic angiography, super-resolution imaging revealed that tumors exhibited characteristic signatures of angiogenesis such as a high degree of tortuosity. Although these early studies were performed in rodents at shallow depths, data were encouraging that this technique would readily translate to deeper tissues in humans as well. The primary challenge with microvascular imaging analysis is that it requires 3-D volume imaging to be most effective, and motion correction can easily corrupt the accurate mapping of microvessels on the order of tens of microns.
The lymphatic system, including lymph nodes and lymphatic vessels, plays a crucial role in oncology. In patients with cancer, tumor cells often traffic through lymphatic channels to regional lymph nodes and beyond to distant sites. Accurate identification and quantification of metastatic lymph nodes remain essential for prognosis and treatment planning in cancer patients (Grills et al. 2003). Recently, super-resolution of the micro-vasculature of rabbit popliteal lymph nodes has been performed   (Fig. 7). Some technical aspects of imaging the slow lymph flow in lymphatic vessels to track the sentinel lymph nodes were also studied. A particular challenge for super-resolved imaging of the tumor vasculature is the complexity of the vessel morphology with a less organized and sometimes immature vessel tree. Adding to this is the higher slice thickness in clinical imaging, which projects more intersecting tracks into a single image, and the lower frame rates, which make the association to tracks more ambiguous. To resolve this, tracking algorithms incorporating linear motion models based on Kalman filters have been proposed (Ackermann and Schmitz, 2016). These methods were applied to pre-clinical tumor xenografts for successful differentiation of tumor types in a radiomics analysis based on quantitative parameters of vessel morphology . Additionally, the clinical feasibility of the method was illustrated for chemotherapy monitoring of breast cancer Dencks et al. 2018Dencks et al. , 2019). An example of microbubble tracking analysis of a contrast-enhanced ultrasound scan of a patient with breast cancer is illustrated in Figure 8. Here, the vascular architecture, displayed in detail, consists of a highly vascularized rim from which vascular assemblies spoke-like proceed toward the tumor center. These are connected to vascular hotspots that exhibit high blood velocities and distribute the blood to the overall tumor tissue.
There still remain several challenges before microbubble tracking analysis is used in clinical routine. The different human organs have distinct functional and anatomic vascular properties that need to be considered to avoid an overload with microbubbles that makes it impossible to assign individual tracks. Furthermore, the resolution in clinical ultrasound scanners is usually lower than in systems used for small animal imaging, as lower frequencies need to be applied to sufficiently penetrate the tissue. Because of the larger voxel sizes and particularly the thicker image slices, it is more difficult to calculate the correct position of the microbubbles, and voxels may even contain several microbubbles at the same time; both make it more difficult to determine correct microbubble tracks. Thus, CEUS protocols and  microbubble do need to be carefully adapted, and consideration of pre-knowledge of the vascular anatomy or the application of deep learning algorithms may help to correctly assign tracks. Another challenge for clinical translation is organ motion. This is less critical for superficial lesions in the skin or the muscle as well as for the brain but becomes highly relevant for abdominal organs such as liver, kidney, spleen and pancreas. Motion correction in the image plane has been implemented (Hingot et al. 2017;Harput et al. 2018;Dencks et al. 2018Dencks et al. , 2019 and needs to be further refined in future studies. However, motion that occurs orthogonal to the image plane remains a challenge. In this context, advancing CEUS imaging from 2-D to true 3-D using matrix transducers, which has been done in vitro (Harput et al. 2019a;Heiles et al. 2019), may not only allow display of the vascular network in much greater comprehensiveness and detail, but also holds promise for a precise motion correction in all directions (Harput et al. 2019b).
In conclusion, clinical application of microbubble tracking analyses is feasible in cancer and provides encouraging preliminary data. Some vendors of ultrasound devices also consider implementing such analysis methods in their scanners' software. In addition, in agreement with research on other non-invasive imaging modalities, there is high interest in the extraction of multiple quantitative image features and their radiomic analysis to receive a disease-specific signature. Recently, it was reported that the radiomic analyses of CEUS data reliably allows differentiation of tumor types in mice . In this context, approaches to microbubble localization and tracking have several advantages over the established methods based on intensityÀtime curves: (i) the imaging protocols are considerably more simple and robust; (ii) there is no need to apply complex pharmacokinetic models based on assumptions that may not be fulfilled (e.g., stable blood concentration of microbubbles or one sole feeding vessel); and (iii) significantly more functional parameters can be extracted than by other CEUS analysis methods because individual tracks and vessels are studied. However, there is still significant room for further research on hardware, algorithms, software and scan protocols to render the methods robust and reliable enough to improve diagnostic precision in clinical routine imaging of cancer.

Neurology
The ability to resolve small vessels deep in the brain has implications in neurology for both diagnostics and therapeutics. Small vessel pathology is increasingly recognized as a contributor in disorders such as stroke and Alzheimer disease (Pantoni 2010;Gorelick et al. 2011;Wardlaw et al. 2013), while tumors are characterized by irregular, tortuous vasculature. Anatomic and functional imaging on a capillary level scale could provide powerful new insight into the workings of the healthy and pathologic brain.
The skull bone stands as the major hurdle to highresolution ultrasound brain imaging. Bone is highly attenuating and aberrating to sound, particularly at higher frequencies (Fry and Barger 1978;Pichardo et al. 2011). Thus, ultrasound imaging and therapy must either be conducted at lower frequencies that better penetrate bone or be performed through acoustic windows, such as the temporal and suboccipital windows (Lindsey et al. 2011(Lindsey et al. , 2014. O'Reilly and Hynynen (2013) used a hemispherical array designed for transcranial ultrasound therapy to image single bubbles in a tube phantom through an ex vivo human skull cap, using the emissions to noninvasively correct image aberration and fitting the expected PSF to generate a super-resolution image (Fig. 9a). This proof-of-concept study illustrated the ability to image through the intact skull bone and the potential to image a large volume of the brain with a spatial resolution <50 um. However, the current approach suffers from long acquisition times because of the low bubble concentration used and the fact that the focal volume required steering to generate the larger image. Employing larger and/or multiple excitation volumes, as well as simultaneously localizing multiple bubbles through fitting or deconvolution (Foroozan et al. 2018) could reduce total acquisition time. In contrast, in their landmark publication, Errico et al. (2015) produced stunning in vivo images of rat brain vasculature at 15 MHz (Fig. 10), both through craniotomy windows and through thinned skull bone.
Although not performed at a frequency amenable to penetrating the human skull bone, this approach could be applied in its current form to image the human brain intra-operatively or in premature infants, where the anterior and posterior fontanelles are still open. This study leveraged ultrafast plane wave imaging and the movement of bubbles between frames to obtain isolated bubbles without having to use very dilute concentrations. Errico et al. (2015) also noted that in theory a reduction of their imaging frequency to 2.5 MHz could still yield 6-mm isotropic resolution. At this frequency, imaging though the thin temporal bone of the human skull could be possible. In fact, multiple matrix arrays imaging simultaneously through the two temporal bones have been used by Lindsey et al. (2011Lindsey et al. ( , 2014 to produce 3-D images of the circle of Willis in human volunteers at normal resolution. Additionally, Soulioti et al. (2020) recently performed transcranial super-resolution imaging of a tube phantom through human temporal bone with a 2.5-MHz diagnostic transducer, achieving an accuracy to that of O'Reilly and Hynynen (2013) and further supporting the feasibility of this approach. The implementation of super-resolution methods on such a system could create a powerful tool for imaging stroke.
A perhaps non-intuitive challenge to be faced in brain imaging, beyond distortion and signal loss caused by the skull bone, is organ motion. Human brain motion, as a result of the cardiac cycle, has been measured to be as high as 0.5 mm, and varies by structure (Enzmann and Pelc 1992). Breathing motion may cause additional lower-frequency shifts in the echo locations. Such displacements are significant relative the structures that super-resolution ultrasound aims to resolve. Motion correction algorithms (Hingot et al. 2017) will thus be necessary to image smaller vessels in the brain.
Beyond diagnostic imaging, the development of super-resolution ultrasound methods has great potential for guiding and monitoring interventions in the brain, such as ultrasound-mediated opening of the bloodÀbrain barrier (Hynynen et al. 2001;McDannold et al. 2006;Carpentier et al. 2016;Lipsman et al. 2018;Mainprize et al. 2019) and sonothrombolysis (Culp et al. 2011;Pajek et al. 2014). Super-resolution imaging thus has immense potential to drastically improve our understanding of the brain and our ability to treat it.   10. Super-resolution image (15 MHz) of rat brain vasculature through thinned skull bone using an ultrafast scanner at 500 frames/s. An acquisition time of 150 s was required to localize and track approximatively one million separable sources per hemisphere. Reprinted with permission from Errico et al. (2015).

OUTSTANDING CHALLENGES AND PERSPECTIVES
As we have seen, super-resolution ultrasound imaging can explore animal and human vasculature in vivo and in depth at a scale unattainable until recently. Both oncological and neurologic applications have been discussed here, leaving many promising applications untouched, such as kidney imaging, which has been explored by several authors (Foiret et al. 2017;Song et al. 2017), and lower-limb imaging, which opens the field of microvascular imaging to diabetic patients .
Despite its numerous advantages, super-resolution ultrasound imaging is still limited by several aspects: acquisition time, signal-to-noise ratio, dependence on exogeneous contrast agents, motion, data overdose, lack of a gold standard, partial reconstruction, exploitation of ultrafast scanners uncommon in the clinic, restriction to the vascular network, lack of 3-D reconstruction, and so forth.
An important challenge to tackle in the field of super-resolution ultrasound imaging is determination of its accuracy. In depth and in vivo, no other imaging modality can provide a similar resolution for vascular imaging. Hence, it is difficult to determine if the observed vessels are appropriately reconstructed, especially with respect to their diameter and velocimetry. One approach is to observe the same plane of a shallow microvascular network with ULM and optical means with an orthogonal setup (Christensen-Jeffries et al. 2015). For this purpose, further studies are needed, potentially with the help of confocal microscopy, multiphoton microscopy or optical coherence tomography to assess the accuracy of ULM at least within the first millimeter. For deep organs, ex vivo 3-D tissue clearing (Renier et al. 2014) could provide further confirmation.
Long imaging times can be partly alleviated through deep learning, but further development is required in ULM to separate closer microbubbles so that injected concentrations can be increased. More localization per frame would mean faster and more resolved images. This is especially true when conventional scanners are used in the clinic. Rapid switching of phasechange contrast agents could also potentially address this challenge.
Vast improvement in signal-to-noise and contrastto-tissue ratios has been achieved in the last 30 y for perfusion imaging with contrast agents. Specific sequences have been developed, but new ones are probably required specifically to extract individual microbubbles. Any gain in SNR and CTR improves drastically the precision of localization and makes microbubbles more distinguishable from background.
One of the key advances will probably come from 3-D super-resolution ultrasound. The lack of information in the elevation direction is particularly limiting in ultrasound localization microscopy, making velocimetry inaccurate and motion correction incomplete and providing only a single plane per acquisition. Fortunately, new results ) on 3-D reconstruction capabilities are encouraging (Fig. 11). However, matrix transducers with numerous parallel acquisition channels are still challenging to fabricate and such technology remains a bottleneck for 3-D ULM, especially at higher frequencies.
Should the field restrict itself to microbubbles? We have seen that nanodroplets or photoacoustic sources could be exploited as localization targets. However, other targets such as cells expressing gas vesicles could be conceived (Bourdeau et al. 2018;Farhadi et al. 2019), allowing the localization of single cells in the body and possibly opening the field of imaging metastatic cells in depth. Fig. 11. Volumetric ultrasound localization microscopy implemented on an anesthetized rat brain using a 2-D matrix array with center frequency of 9 MHz . The voxel size is 10 mm; the region of imaging is 13.28 mm deep, 10.5 mm wide and 11.4 mm long. The "glow" color map was used to encode microbubble density. The brighter the color of the voxel, the larger is the amount of microbubbles have passed through that voxel.
Should the vasculature remain the only playing field for super-resolution ultrasound imaging? The lymphatic system is already being explored . But the extravascular space could be the next frontier using the natural permeability of tumor vasculature to nanoagents, which could become localization targets after acoustic droplet localization.
Finally, the largest challenge remaining is probably definition of the optimal application of super-resolution ultrasound imaging, the application that would become irreplaceable by other imaging modalities. We will only rest-briefly-when the prognoses of patients have vastly benefited from the use of super-resolution ultrasound imaging.

CONCLUSIONS
In the last 10 y, super-resolution ultrasound has disrupted the diffraction limit of ultrasound imaging through techniques such as ultrasound localization microscopy. By observing directly the microvasculature in the brain, kidney, skin, lymph nodes and tumors in animals and, in some cases, humans, it has provided a new window through which to study the physiology and pathology of tissue. In-depth microscopy is now achievable with a clinical modality, which can provide a wealth of new information, such as the density and tortuosity of microvessels or even minute modifications in flow patterns. The authors referred to in this review have participated in the early development of super-resolution ultrasound imaging. Its scope has now widened and international endeavors will define the various clinical applications of the technique, targeting the major killers, which affect the smallest blood vessels of the body, such as cancer, diabetes (Ghosh et al. 2019) artheriosclerosis. In the future, we foresee new volumetric ultrasound scanners, which could map the 3-D microvasculature of an entire organ, allowing precise diagnosis with reduced user dependency. College Confidence-in-Concepts.
Human studies cited in this article were approved by an institutional review board (IRB) except for historical studies. IRBs were introduced after the first revision (1975) of the Declaration of Helsinki. Likewise, if animals were studied, there was institutional animal care committee approval.
Declaration of Competing Interest-OC and MT declare being coinventors on a super-resolution patent (PCT FR2011/052810). MT is a co-founder and shareholder of Iconeus company for ultrasound neuroimaging. PAD declares co-inventorship on patents and applications for super-resolution techniques, dual-frequency imaging (US9,532,769), and phase change nanodroplets (US9,427,410) and is a co-founder of SonoVol, Inc., and Triangle Biotechnology, Inc., which have licensed some of these patents. KH and MAO declare being co-inventors on a super-resolution patent (PCT/US2014/036567).