Super-Resolution Ultrasound Bubble Tracking for Preclinical and Clinical Multiparametric Tumor Characterization

Super-resolution imaging methods promote tissue characterization beyond the spatial resolution limits of the devices and bridge the gap between histopathological analysis and non-invasive imaging. Here, we introduce Ultrasound Bubble Tracking (UBT) as an easily applicable and robust new tool to morphologically and functionally characterize fine vascular networks in tumors at super-resolution. In tumor-bearing mice and for the first time in patients, we demonstrate that within less than one minute scan time UBT can be realized using conventional preclinical and clinical ultrasound devices. In this context, next to highly detailed images of tumor microvascularization and the reliable quantification of relative blood volume and perfusion, UBT provides access to multiple new functional and morphological parameters that showed superior performance in discriminating tumors with different vascular phenotypes. Furthermore, our initial clinical results indicate that UBT is a highly translational technology with strong potential for the multiparametric characterization of tumors and the assessment of therapy response.


Introduction
Ultrasound (US) is among the most frequently used diagnostic modalities in clinical routine, and its spatial and temporal resolution as well as tissue contrast have been steadily improved. The application of gas-filled microbubbles (MB) as US contrast agents further enhances the diagnostic accuracy of US by adding morphological and functional information about the tissue vascularization 1 . This is particularly relevant in oncology, since the vascular structure of tumors contains essential information for their differential diagnosis 2-4 , prognostication 5 , and for the prediction and monitoring of therapy responses [6][7][8] . In particular, some vascular features have already been shown to be capable of identifying patients not responding to antiangiogenic therapy 9 , who, then, can be reoriented towards alternative approaches 10 .
Different qualitative and quantitative techniques have been developed to extract the information about tumor vasculature contained in contrast-enhanced US (CEUS) scans. However, in state-ofthe-art CEUS imaging, e.g. using Maximum Intensity Over Time (MIOT) 11 or replenishment kinetics analysis 12 , voxels are usually much larger than the majority of tumor blood vessels, whose diameters are in the range of 5-80 µm 13 . This limitation in the spatial resolution makes it difficult to gain a comprehensive overview of the vascular architecture and its heterogeneity. In addition, since the probability is high that every voxel contains at least one blood vessel, the tumor vascularization tends to be overestimated whenever the relative blood volume (rBV) is determined based on the area that exhibits MB signals 14 . Voxel-wise analyses are further complicated by high background noise, which can make the assessment of functional vascular parameters difficult and unreliable at the single voxel level 15 .
To overcome these issues, several postprocessing algorithms for CEUS image analysis have recently been proposed to reveal and quantify vascular features at super-resolution, which means at a resolution beyond the resolution limits of the device 16,17 . Here, individual MB are localized, and a line with the thickness of a MB is drawn connecting the most closely localized MB in two subsequent frames. This line represents the track of a MB and thus, the course of a (micro) vessel. The approach was successfully applied to characterize MB flow tracks in brain 16 and ear vessels 17 . However, in case of ambiguous assignment possibilities, this approach could lead to underestimation of flow velocities and might be particularly difficult to apply to more complex tumor vascular networks. Therefore, Errico and colleagues 16 used an experimental imaging system with a very high frame rate (500 frames per second) to avoid ambiguous assignments.
However, comparable frame rates are not realized in the majority of clinical US systems so far, which makes clinical translation of this method difficult. Therefore, we present here an alternative super-resolution CEUS approach called "Ultrasound Bubble Tracking (UBT)", which is an advanced tracking technique that is adapted to clinical settings. With UBT, within less than a minute and using a conventional US devices operating at standard frame rates, super-resolution images and novel parameters could be extracted, which enabled an accurate discrimination of tumors with different vascular phenotypes. Furthermore, the preliminary clinical data that are presented in this manuscript show that rapid translation of UBT is realistic and that this technology may improve the diagnostic potential of CEUS in future clinical practice.

Ultrasound Bubble Tracking for structural imaging of tumor vasculature
The UBT method reliably captured the movement of MB in tumors and could be successfully applied to all contrast-enhanced scans using a commercial US device operating at frame rates of approximately 50 frames/s and only required measurement times of 40 s (for more information about the algorithm, see Fig. 1a. as well as the detailed description in the "Materials and Methods" section). A spatial resolution of approximately 5 µm was achieved, and the vascular architecture of tumors was visualized in fine detail (Fig. 1b). Functional parameters were calculated for single vessels and combined with textural features, which so far could not been obtained with standard CEUS imaging (Fig. 1c).
At super-resolution, differences in the vascular texture in different tumor models could be clearly depicted. For instance, in line with histological staining, the highly angiogenic A431 tumors showed a fine network of very small vessels, homogeneously distributed throughout the entire tumor tissue (Fig. 1b, 2). In contrast, A549 tumors, which are known to be less vascularized and characterized by a more mature vascular system, displayed a higher vascular hierarchy in superresolution UBT images, with larger vessels at the periphery, that branch into smaller vessels towards the tumor center. MLS tumors with their heterogeneous vascular pattern were most difficult to classify. In the super-resolution UBT images, as in histology, they were characterized by highly and poorly vascularized regions, and by more or less dense and branched vascular areas (Fig. 2).
As an additional feature of the UBT approach, velocities of MB and their directions of movement can be calculated for individual vessels and displayed in parametric maps (Fig. 1b). In these parametric direction maps information about arterial and venous supply as well as 6 branching and connections between vessels is provided. While the velocity profiles of the tumors were similar, the analysis of flow directions showed differences. The parametric direction maps of A431 tumors indicated that blood flow in the central fine vascular network was predominantly directed towards the tumor core, while in the periphery the blood flow directions were chaotic ( Fig. 1b). Thus, a balance between feeding and draining vessels was found in the periphery, while draining vessels were less apparent in the center. This lack of venous drainage is known to be a typical characteristic of highly angiogenic tumors 18 . In contrary, the parametric maps of A549 and MLS tumors, which have a more mature vascularization ( Supplementary Fig. 1), showed a balanced mixture of feeding and draining vessels (Fig. 2c), which is also reflected by their higher flow direction entropy values (Fig. 3b).

Characterization of vascular tumor phenotypes using Ultrasound Bubble Tracking
While some UBT parameters can also be obtained by state-of-the-art CEUS postprocessing techniques, others represent new parameter classes that so far have been difficult to assess (Fig.   1c and Supplementary Table 1). Parameters determined by UBT include the relative blood volume (rBV), the mean, variance, maximum and median values of MB velocities, distances to the closest vessel, and distances to vessels with low and high velocities, as well as the flow direction entropy as a measure for the organization of the vessel networks.
In good agreement with their histological vascular phenotypes, A431 tumors had the highest rBV values, followed by MLS and A549 tumors (Fig. 3a). However, due to the high heterogeneity of rBV within the respective tumor models, our group sizes were not sufficiently large to generate significant differences.
As expected from the visual inspection of the parametric maps (Fig 2c), the flow direction entropy values, describing the order of blood direction profiles, were lower in A431 than in MLS and A549 tumors. While this parameter could not unambiguously separate all groups, significant differences were found between A431 and MLS tumors (Fig. 3b).
Surprisingly, parameters related to MB velocity were very similar across the tumor models, indicating that, despite their different angiogenic phenotypes, these tumors tend to preserve a similar flow pattern (Fig. 3c).
The textural parameters of the tumor vascularization showed a significantly higher discriminatory potential ( Fig. 3d-f). In this context, super-resolution images obtained by UBT enabled us to determine the distances to the closest vessel and to evaluate their mean, variance, maximum and median values. Strikingly, the first three of the above parameters had the power to discriminate all three tumor groups. However, among all distance parameters, the maximum of distances to the closest vessel was one of the best performing ones, which precisely discriminated all three tumor models (A431 vs. MLS, p< 0.01; A431 vs. A549, p< 0.001 and A549 vs. MLS, p<0.01) (Fig. 3d).
Two new parameter classes were introduced, which combine textural and functional information, i.e. 1) distances to vessels with low velocities and 2) distances to vessels with high velocities.
While based on the parameters associated with distances to vessels with high velocities only one or two out of three possible combinations revealed significant differences, mean and maximum values of distances to vessels with low velocities differed significantly between all tumor models ( Fig. 3e,f). Thus, the latter parameters were considered in the further analysis.

Confusion and correlation matrices of novel parameters obtained by UBT
Parameters that could distinguish all three tumor models by statistical evaluation with one-way ANOVA and Bonferroni post-hoc tests were used for further analyses. These parameters were: 1) mean, 2) variance and 3) maximum of distances to the closest vessel as well as 4) mean and 5) maximum of distances to vessels with low velocities (Fig. 4a).
In order to determine the discriminatory power of the parameters at the basis of individual tumors, for each of these parameters the nearest neighbor classifier (NN) was applied in a leaveone-out-cross-validation, and confusion matrices were generated (Fig. 4b). It is clearly indicated that the maximum of distances to the closest vessel and maximum of distances to vessels with low velocities were best suitable for classifying three tumor models. With both parameters, a completely correct classification of all tumors was achieved (100%). With variance of distances to the closest vessel, 83% of the classifications were correct. However, all A431 tumors were classified correctly, and only one MLS tumor was wrongly assigned as an A549 tumor and one A549 tumor as a MLS tumor. Furthermore, with the mean of distances to the closest vessel and mean of distances to vessels with low velocities, 67% and 58% of the correct classification could be achieved, respectively (Fig. 4b).
In order to decide which parameters should be combined to correctly classify the tumors, it is important to investigate their interdependence. The lower the correlation between parameters that have high distinctive power, the higher is the probability that they will provide complementary information. For this purpose, a correlation matrix was generated. The superior distance parameters were all highly correlated and one single out of these parameters was sufficient to distinguish all tumors. Thus, a combination of parameters was not required. Nevertheless, it may become necessary to combine parameters during examinations of animals or patients with more heterogeneous tumors. In this context, a promising parameter is the flow direction entropy, since it can significantly distinguish two tumor groups and shows low correlations (r<0.5) with the distance parameters (Fig. 4c).

Comparison of parameters derived from UBT and reference methods
To assess the robustness and the accuracy of UBT, we firstly compared the level of tumor vascularization (rBV) obtained by UBT to rBV values obtained by three other techniques, i.e.
MIOT postprocessing, ex vivo micro computed tomography (µCT) and immunohistochemical (IHC) analysis of the tumor sections. Although rBV values did not differ significantly across the three tumor models, all modalities showed the same trend, classifying A431 tumors as the most vascularized ones, followed by MLS and A549 tumors (Fig 5a-e). However, at a quantitative scale rBV determined by MIOT revealed higher, µCT comparable and IHC slightly lower values than UBT (Fig. 5e).
In the next step, we compared mean velocities obtained by UBT with mean velocities calculated from replenishment kinetics. We found that both methods did not show differences in perfusion among the tumor models and provided values clearly below 1 mm/s. However, while UBT indicated mean velocities of approximately 0.8 mm/s, the values obtained by replenishment analysis were systematically lower (approximately 0.09 mm/s) (Fig. 5f).
Finally, the quantitative values (mean, variance and max) of distances to the closest vessel obtained by IHC und UBT presented with the same order, i.e. A431 had the smallest, A549 the largest and MLS intermediate values (Fig. 5g).

Clinical proof of concept
In order to demonstrate that UBT can be performed using conventional US devices and clinically approved US contrast agents, three patients were investigated. The first patient was a 53-year old woman with an invasive "no special type (NST)" breast cancer examined after the third cycle of neoadjuvant chemotherapy with epirubicin and cyclophosphamide. The patient was scanned in B-mode using a 12 MHz transducer of the Toshiba Aplio 500 device (Toshiba Medical Systems GmbH, Otawara, Japan). Three milliliters of SonoVue (Bracco, Milan, Italy) were injected slowly over 6 min. Although the data were not acquired with a contrast-specific mode, UBT successfully captured a large vascular trunk with many feeding and draining vessels and thus identified the tumor areas that were still viable (Fig. 6a).
The second patient was a 30-year old woman with a thyroid nodule, scanned with a Zonare ZS3 The patient was repeatedly imaged after the slow injection of 0.5 ml of SonoVue, before the initiation of the chemotherapy, after the first and after the second cycle with a 12 MHz transducer of the Toshiba Aplio 500 device. UBT super-resolution images nicely displayed the tumor vasculature and depicted the change in blood perfusion and flow direction over the course of treatment. At the baseline measurement, the vascularization was mainly located in the central 11 part of the tumor without a dominant blood flow direction. Surprisingly, after the first chemotherapy cycle, we could observe apparent changes in the vascularization pattern. The tumor vascularization appeared much more homogeneous and strongly enhanced at the periphery. Additionally, rBV increased while the tumor size decreased from 25.6 cm 3 to 3.0 cm 3 .
After the second cycle of the chemotherapy the tumor size decreased remarkably to 0.8 cm 3 and also its vascularization was substantially lower (Fig. 7). We hypothesize that tumor cell death induced by the first chemotherapy administration decreased the solid stress and/or interstitial fluid pressure 19 , which then caused vascular decompression and thus, improved tumor perfusion and drug delivery in subsequent chemotherapy cycles.

Discussion
In this study, we evaluated UBT for the structural and functional imaging of vascular features in murine and human tumors at super-resolution. We show that UBT opens new avenues for textural and functional tumor analysis at the individual vessel level and it provides novel classes of vascular parameters. In this context, reliable and robust quantification was achieved and different (vascular) phenotypes of tumors could be accurately discriminated. We postulate that this comprehensive and quantitative vascular characterization can be clinically highly valuable since the level of vascularization, microvessel density, and the distribution of vessels are often highly correlated with tumor invasiveness, aggressiveness, metastatic potential and prognosis of the disease 20,21 , as already shown in different types of tumors e.g. brain tumors and melanoma 22,23 .
When comparing UBT to reference techniques, we found that rBV values of tumors obtained by MIOT, µCT and IHC showed the same trend across the tumor models. At a quantitative scale rBV values obtained by UBT and µCT were very similar. However, as expected, MIOT overestimated the rBV since this technique counts every US voxel showing a positive MB signal as vessel, even if the vascular fraction within the respective voxel is small. In contrary, the somewhat lower rBV values obtained by IHC are explained by the fact that we preserved samples in formaldehyde, which is known to lead to tissue shrinkage 24,25 .
Subsequently, mean velocities obtained by UBT and replenishment kinetics analysis were compared. In this context, it should be noted that in 25% of the cases, replenishment curves could not be fitted due to high noise levels in the US images, while all measurements were reliably postprocessed with UBT, which demonstrates its higher robustness. Nevertheless, both methods indicated that velocity values of the three tumor models were very similar. However, the values obtained by replenishment kinetics analysis were systematically lower than in UBT.
Considering velocity values of mouse tumors reported in literature (1.1-1.5 mm/s) from multiphoton laser scanning microscopy 26 , the quantitative numbers provided by UBT appear more realistic. In line with this, in previous publications the authors already reported that quantitative values obtained by replenishment kinetics analysis may not always be absolutely accurate [26][27][28] . This may be explained by the fact that within a region of interest blood flow in the image plane cannot be detected and therefore remains unconsidered 28 . Furthermore, the majority of replenishment analyses do not consider the influence of the beam elevation characteristics on the replenishment curve shape, which may also make the velocity values less accurate 27 .
There is substantial need for further improvement of the UBT technique. In this context, our first patient measurements indicated several issues that, as long as unresolved, stand in the way of a broad clinical implementation. All measurements suffered from tissue motions that need to be compensated. While this is considerably easy for in-plane motion, out-of-plane movements cannot be corrected in 2D measurements, except when removing the non-matching slices, which leads to a loss of valuable data. Furthermore, the injection speed and concentration of MB need to be optimized. In case of the thyroid nodule, for example, the injection rate was too high, and individual MB could hardly be distinguished. Therefore, we used the CEUS sequences acquired in the early phase of the injection, which reduced the number of exploitable image frames.
Consequently, some vessel trees might not have been completely reconstructed. Thus, CEUS scans for UBT should be performed under slow MB injection and will require much lower MB doses than for conventional methods, which in turn, however, may help to decrease potential side effects and concerns that have been raised, e.g. for using CEUS in patients with instable cardiopulmonary conditions and pulmonary hypertension 29 . Furthermore, our data suggest that in clinical settings, contrast-specific modes should be applied to acquire US data for UBT, since the low signal-to-noise ratio of the B-mode, as used for the first case (breast cancer), made the MB detection more difficult than in the contrast mode scan used for the second (thyroid nodule) and third patient (neoadjuvant breast cancer therapy). Additionally, we observed relatively high UBT may not only be applied in oncology but may also be relevant for other indications, e.g. for characterizing inflamed tissues (e.g. inflammatory bowel disease), organ fibrosis (e.g. in liver and kidney), and immunological disorders (e.g. rheumatological disorders and organ rejection after transplantation). Furthermore, it may be applied to monitor vascular remodeling in the cardiovascular field (e.g. revascularization in ischemic tissues) and to assess antiangiogenic treatment effects e.g. in retinopathies. Based on the presented data, we are confident that, after further refinements, UBT will experience a rapid translation into clinical practice.

Study design
The objective of this study was to establish UBT as a CEUS postprocessing method for distinguishing tumors with different vascular phenotypes at super-resolution, as well as to provide proof of principle for applying UBT on clinical data. The study consisted of three major parts. In the first part, the ability of UBT to distinguish tumor models with different vascular phenotypes was assessed 32 . For this purpose, A431, MLS and A549 tumor xenografts were induced in female immunodeficient CD1-nude mice (n=4 mice per tumor model). CEUS imaging was performed, and various parameters were extracted that describe morphological and functional vascular characteristics. Statistical tests were applied to investigate the diagnostic potential of these parameters for discriminating three tumor models and to find their ideal combination. In the second part of the study, we evaluated the robustness of UBT and the accuracy of the outcome parameters using available reference techniques. For that purpose, CEUS cine loops were analyzed with MIOT and replenishment kinetics to assess rBV values and mean velocities, respectively. For further validation of the rBV values, high resolution µCT scans of Microfil perfused tumors and histological analyses of tumor sections were evaluated (n=4 mice per tumor model). In the third part of the study, we applied our technology on CEUS data from three patients, in order to demonstrate its translational potential.

Study approval
All animal experiments were approved by the local and governmental committee on animal care.
The patients' examinations with CEUS were performed as part of the clinical routine diagnostic procedure. All patients gave their informed consent to retrospectively extract the relevant US images from the stored routine data and use them for the image postprocessing.

Xenograft tumor models
The mice were housed in groups of four per cage under specific pathogen-free conditions with a 12h light and dark cycle in a temperature-and humidity-controlled environment (according to

Contrast-enhanced ultrasound imaging
Hard-shell polybutylcyanoacrylate (PBCA) MB were used as US contrast agent to assess the potential of UBT for tumor characterization. PBCA-MB were freshly synthetized as described before 33 . For animal experiments, the PBCA-MB suspension was diluted in sterile sodium chloride to a concentration of 2x108 MB/ml. Each mouse was injected with a 50 µL bolus containing 1x10 7 PBCA-MB over approximately 3 seconds, followed by a 20 µl saline flush, into the lateral tail vein. The correct placement of the probe on the tumors was controlled prior to the measurements using real-time B-mode imaging. For the measurements, image series were recorded during the destruction-replenishment sequence at a frame rate of 50 frames/s. The images were acquired in digital raw radio frequency (RF) mode to receive uncompressed IQ data. The gain was set to 22 dB and the transmit power to 2% to minimize MB destruction. For each tumor, the number of processed frames was limited to 2000, which was equivalent to 40 s measurement time.

Image processing
After the US measurements, the tumor border was outlined manually and the further processing steps were carried out inside the segmented area.
In a first step, a rigid motion estimation and compensation were carried out. For this, the B-mode images were interpolated on a finer grid (4-fold) to increase the accuracy of the motion compensation. The motion profiles over time typically exhibited periods of small movements disturbed by spikes of large movements due to breathing. These frames of large movements were excluded because they typically included also out-of-plane movements.
To detect single echoes of individual MB, the B-Mode images were separated into static background and moving foreground images. The rolling background was calculated applying a temporal rank filter of rank 3 over ± 10 frames around the actual frame. The foreground was computed by subtracting the background from the original frames. After applying an adaptive threshold to the foreground images, the MB were localized by calculating the intensity weighted centroid for each MB. Choosing a too low threshold leads to false detections due to noise in the image, choosing a too high threshold leads to missed detections. Therefore, for each measurement set, the threshold was adapted to result in no detections immediately after the destruction event in the destruction replenishment sequence.

Tracking of microbubbles
To reconstruct the vasculature, the detected MB must be tracked over several frames. As discussed in 16 , usually this is either achieved by a nearest-neighbor tracking for very high frame rates in the kHz range, which need non-standard ultrafast US scanners, or by very low MB concentrations and long observation times up to several minutes.
Here, we solved the tracking challenge that occurs when using clinically recommended MB concentrations and standard US systems, by using a novel, more robust UBT method. This was necessary, since tracking by using the nearest MB in the next frame to continue a track is prone which also led to a higher a posteriori probability .
The tracking algorithm yielded not only the trajectories but also the velocities and directions of the moving MB.

Definition and extraction of parameters
For further evaluations, images of MB tracks, velocities, and flow direction were reconstructed with a pixel size of 5 × 5 µm 2 for each tumor. Even though we were able to achieve a higher resolution with the 40 MHz transducer, we used these settings due to the MB's size (3µm), the limited resolution in the graphical illustrations and to restrict the data amount.
Microbubble track images were generated using Bresenham's line algorithm to connect the MB positions of the estimated tracks. The MB track image was a binary image indicating the pixels which were passed by the MB. In the flow velocity map and the flow direction map, the corresponding quantity was assigned to the pixels along the track.
For the evaluation of the parameters, each tumor was divided into two regions: a rim of 0.5 mm thickness and the core which was the remaining area when excluding the rim from the whole tumor area. For the characterization of the tumor vasculature, we used only the core region to exclude the large feeding vessels in the rim.
From the MB track image the rBV was derived as the ratio of the area covered by the tracks to the respective total area, which was expected to be proportional to the rBV.
Additionally, from the MB track image a track distance map was generated applying the Euclidean distance transform (bwdist function, Matlab, MathWorks, Natick, MA, USA). For each pixel, the track distance map provided the shortest distance to the next vessel. For each track distance map, mean, variance, maximum, and median of the distances to the closest vessel were calculated for the respective areas. Small distances are characteristic for a fine meshwork of vessels. The larger the maximum distance, the larger are the non-perfused areas.
From the flow velocity map the statistics of the velocities were derived. Again, mean, variance, maximum and median were calculated.
Since we were interested in the structure of the vasculature, we divided the vessels into two groups of high and low flow velocities, respectively, and calculated the distance parameters separately for the resulting two groups. We defined the mean value of the median velocities of the tumors as the threshold between high and low flow velocities (0.7 mm/s) and calculated the mean, variance, maximum, and median of the distances to vessels with low and high velocities.
We were also interested in parameterizing the directions of MB flow. To characterize a locally

Statistical analysis
Data are presented as mean ± standard deviation (s.d). The one-way ANOVA and Bonferroni post-hoc test were applied to evaluate differences between groups considering a p-value of <0.05 to be significant. All analyses were performed using GraphPad Prism 5.0 (GraphPad Software, San Diego, CA).

Confusion and correlation matrices
For all parameters extracted by UBT, one-way ANOVA and Bonferroni post-hoc tests were applied. Parameters that distinguish all three tumor models were used for the further analyses.
Then, for each parameter the nearest neighbor classifier (NN) was applied in an exhaustive leave-one-out-cross-validation. The results are presented in confusion matrices, which plot the actual classes versus the classes predicted by the classification. The numbers in the diagonal elements of the matrix represent the correct classifications; the remaining numbers indicate the false assignments. The correct classification is expressed in percentage.
Additionally, a correlation matrix of all parameters obtained by UBT was generated to depict the pairwise dependencies among them measured by the Pearson's correlation coefficient (r).
Confusion and correlation matrices were generated by Matlab, MathWorks, Natick, MA, USA.

Reference methods
We

Data availability
The data that support the findings of this study are available from the corresponding authors on request. yellow dot: the information is available but its assessment is less robust, less accurate or not quantitative; red dot: the parameter cannot be obtained with the respective method).   Therefore, a correlation matrix (Pearson's correlation coefficient (r)) of all UBT parameters was generated to indicate the parameters providing complementary information. The highly discriminating distance parameters strongly correlated and thus, their combination may not be advantageous. However, the parameter flow direction entropy that distinguished two tumor models showed a low correlation with the distance parameters and could be selected as a potential candidate for a multi-parameter readout. While these parameters obtained by UBT analysis were significantly different across tumor models, by IHC analysis we could observe significant difference in A431 tumors compared to A549 (*=p<0.05, **=p<0.01) (data are presented as mean ± s.d. one-way ANOVA with Bonferroni post-hoc analysis).