Bayesian Quantification of Contrast-Enhanced Ultrasound Images With Adaptive Inclusion of an Irreversible Component

Contrast Enhanced Ultrasound (CEUS) is a sensitive imaging technique to assess tissue vascularity and it can be particularly useful in early detection and grading of arthritis. In a recent study we have shown that a Gamma-variate can accurately quantify synovial perfusion and it is flexible enough to describe many heterogeneous patterns. However, in some cases the heterogeneity of the kinetics can be such that even the Gamma model does not properly describe the curve, with a high number of outliers. In this work we apply to CEUS data the single compartment recirculation model (SCR) which takes explicitly into account the trapping of the microbubbles contrast agent by adding to the single Gamma-variate model its integral. The SCR model, originally proposed for dynamic-susceptibility magnetic resonance imaging, is solved here at pixel level within a Bayesian framework using Variational Bayes (VB). We also include the automatic relevant determination (ARD) algorithm to automatically infer the model complexity (SCR vs. Gamma model) from the data. We demonstrate that the inclusion of trapping best describes the CEUS patterns in 50% of the pixels, with the other 50% best fitted by a single Gamma. Such results highlight the necessity of the use ARD, to automatically exclude the irreversible component where not supported by the data. VB with ARD returns precise estimates in the majority of the kinetics (88% of total percentage of pixels) in a limited computational time (on average, 3.6 min per subject). Moreover, the impact of the additional trapping component has been evaluated for the differentiation of rheumatoid and non-rheumatoid patients, by means of a support vector machine classifier with backward feature selection. The results show that the trapping parameter is always present in the selected feature set, and improves the classification.


I. INTRODUCTION
C ONTRAST-ENHANCED Ultrasound (CEUS) is an imaging technique for the visualization and assessment of tissue vascularization and perfusion, using microbubbles as contrast agent. In particular, CEUS has been widely studied for the early detection and grading of arthritis diseases, since it allows a non-invasive assessment of the synovial neovascularization and of the variations in local perfusion [1]- [7]. CEUS data are generally quantified at the region of interest (ROI) level, i.e. analyzing the time intensity curve (TIC) obtained by averaging all pixel TICs within a specific userdefined region. In a previous work, we have tackled the CEUS quantification problem at the pixel level and we have demonstrated that pixel-wise quantification allows such an effective characterization of different perfusion patterns that they can be then used to discriminate between arthritis subtypes [8]. Moreover, we showed that a physiologically motivated perfusion curve, such as the Gamma-variate model, can accurately quantify synovial perfusion and it is more flexible than a mono-exponential or logarithmic model (generally employed for CEUS kinetic description [9]) to describe many heterogeneous patterns [8]. In fact, even though the mono-exponential model was originally presented the model of choice for the region-based curve [1], [10], [11], it does not represent the appropriate model for the analysis at the pixel level and in this case the Gamma-variate model is to be preferred [9], [12].
However, there are two issues that were not fully addressed in our previous work: first, in some cases the heterogeneity of the CEUS kinetics can be such that even a Gammavariate model does not properly describe the time courses. One possible source of this heterogeneity, generally ignored in CEUS data modeling, can be an apparent trapping of the microbubbles in newly formed vasculature, when characterized by high tortuosity [13], [14]. This apparent trapping can also be explained by the limited temporal duration of the CEUS scan, that prevents the identification of a very slow washout from the data. This problem has already been considered in other contrast-based imaging techniques (e.g. in dynamic susceptibility contrast magnetic resonance imaging, DSC-MRI), where the original Gamma-variate was augmented with its integral to explicitly model the recirculation of the contrast agent (single compartment recirculation, SCR) [15], [16]. However, nothing similar has ever been investigated for CEUS images quantification.
The second issue is the choice of the pixel-wise estimator of the perfusion kinetic model to perform parametric mapping. The use of a non-linear estimator, the standard method for ROI-based quantification, leads to a high number of unreliable estimates (between 40% and 50% of pixels) when applied at pixel-level, thus losing an important part of the physiological information in the synovia [8]. A valid alternative for pixelwise quantification is represented by Bayesian methods, which incorporate prior information on the tissue kinetic and are widely used in other imaging techniques but seldom applied to CEUS data.
In the current work, our aim is two-fold: first, to improve the CEUS data quantification by including an irreversible component in the model, and, second, to present a more robust and efficient estimation algorithm for the derivation of the pixel-wise estimates of perfusion information.
In particular, we applied the SCR model of [15], [16] to CEUS data, but we solved it at pixel level within a Bayesian framework using a Variational Bayesian estimator [17], already applied to MR [18]. Bayesian modeling allows to formally incorporate prior knowledge (in form of probability distributions) into the estimation process. In this work, we set up a user-independent setting of the priors, that were generated automatically from the image data, in an empirical fashion [19]. First, the model under study was applied to a set of synovial sub-region time courses solved with a nonlinear estimator (that represent the gold standard for region-wise analysis). These estimates were then used to determine the prior information for each sub-region, thus obtaining tailored priors for each region.
The SCR model here proposed is more complex than the single Gamma-variate model previously proposed in [8]. In order to avoid increasing the model complexity where it is not necessary, we implemented the Variational Bayesian method including the automatic relevance determination (ARD) algorithm [18], [20], [21] to exclude the trapping where not supported by the data.
We used clinical CEUS data from patients with arthritis disease to validate the use of the SCR model and assess the VB estimator performance in clinical settings. We further employed the perfusion estimates obtained to train a support vector machine (SVM) classifier to distinguish different types of arthritis. In particular, we focused in the analysis of rheumatoid (RA) arthritis versus other types of less aggressive arthritis forms. Rheumatoid arthritis represents the worst outcome among the different types of arthritis, causing premature mortality, disability and compromised quality of life [8]. In a classical view of rheumatic disease, rheumatologists do not need automatic methods to diagnose rheumatoid arthritis. However, the validation and application of precise and quantitative imaging biomarkers is critical for the early assessment of the patient's response to treatment, and for the differentiation of rheumatic disease at their early stage, where the clinical findings are most similar.
A precise quantification of synovial perfusion and synovial angiogenesis is thus becoming of increasingly importance in the field of vascular-targeted therapeutics for the treatment of arthritis, that is now in its pre-clinical phase [22]. Moreover, from a translational point of view, the coupling of a precise description of perfusion with biomarker discovery methods could help in understanding the biological basis of rheumatic arthritis, highlighting the basis of the pathogenesis and progression of RA.

II. MATERIALS AND METHODS A. Dataset
Data from previously reported studies [8] of 115 subjects affected by arthritis were used in the current study. The criteria for subject inclusion and the procedure for CEUS exams are described in detail in [8].
All patients showed signs of inflammatory finger joint involvement: the joint with the highest disease activity was chosen for CEUS examination. All patients gave their informed consent to the examination, to the intravenous administration of the contrast agent and to the participation of the study that was approved by the local institutional ethical committee. Each subject underwent a 2-min CEUS study with a 7-MHz transducer US device (MyLab25, EsaOte, pixel size 0.6 mm, dynamic range between 45 and 63 dB, gain between 0 and 20 dB) equipped with Contrast tuned Imaging (CnTI; Esaote), after a bolus injection of 4.8 ml of microbubbles contrast agent. Gray-scale US (anatomical B-mode image) was acquired before the bolus to define the boundaries of the synovial tissue (semi-automatically outlined by the radiologists as in [23]).
CEUS images were motion corrected and coregistered to the anatomical image as described in [8]. We also linearized the log-compressed data as in [24]- [26] and normalized each pixel curve by the maximum value in the image data, in order to have the data in the range [0; 1].
From the synovial boundaries, we derived an image mask that was then modified to include only pixels that showed a significant enhancement, considered as the difference between baseline value and peak value larger than a fixed threshold, set to 0.2 (20% of the data maximum value) based on previous studies [27]. We excluded sixteen subjects with no active pixels from the analysis. Our final dataset comprised 99 subjects (53 RA, 21 dpPSA, 9 RA-like PSA, 8 uSPA and 8 CTD).
The kinetic analysis was performed on the specific region of interest of the synovia.

B. Gamma-Variate and Irreversible Component Model
The model for the vector of CEUS measurements c pi xel (t) (i.e. the tracer concentration in a pixel over time) is: where G( p, t) is the model with the parameter vector p (which represents all the individual rate constants of the model) and e(t) is the measurement error. For the sake of clarity, the dependency on the time t has been omitted in the following passages. We assumed the measurement error to be additive, independent and identically distributed with a normal density function with zero mean and φ −1 unknown variance, i.e. e ∼ N(0, φ −1 I), as defined in [18]. We propose two different models G( p, t) for CEUS data quantification, the Gamma-variate and the SCR model.
The Gamma-variate model equation is: where a 0 is a scaling factor, α and β determine the bolus shape (raise and washout of the dye from the vascular bed) and t 0 is the contrast arrival time. The amplitude a 0 of the Gamma-variate function is normalized for c gamma max , the maximum value of c gamma (t) (i.e. the peak value), in order to have a 0 varying between 0 and 1. From the Gammavariate estimated microparametersp gamma = [â 0 ,α,t 0 ,β] it is possible to derive a set of additional macroparameters, such as the peak value G max , the time of peak t max , the raise time t raise and the washout time t wash (computed as the time needed to raise the intensity from the baseline value to half maximum and from the peak value to half maximum respectively), the mean transit time MT T , the blood flow index B F I and the blood volume index BV I , for a total of 11 parameters for each curve. The equations for the calculation of the macroparameters are reported in Appendix A.
The classical Gamma-variate model will consider negligible any apparent trapping (that can be seen as an irreversible component) of the contrast agent, compared to the general picture of the system dynamics.
An alternative model to describe CEUS kinetics is the single recirculation component (SCR), previously proposed for DSC-MRI perfusion imaging. This model can be used to describe the irreversible component in CEUS: where a 1 is the scaling factor (normalized by the area under the curve AU C g of c gamma (t)) to range between 0 and 1). The integral term describes the trapping term (i.e. the fraction of contrast trapped in the vasculature). From the set of estimated parametersp SCR = [â 0 ,α,t 0 ,β,â 1 ] for the SCR model, it is possible to derive the same macroparameters obtained for the Gamma-variate model, i.e. G max , t max , t raise t wash , MT T , B F I and BV I , for a total of 12 parameters for each curve (Appendix A).

C. Variational Bayes Theory
Both the Gamma-variate and SCR model have been identified separately at the pixel level using a Variational Bayesian (VB) estimator, as proposed in [18].
In the Bayesian framework, one is interested in the posterior distribution of the parameters given the data and the model (P( | y, w)) (or P( | y) written neglecting the dependence on the chosen model w), where is the vector of the parameters of a chosen model w from a set of measured data y. In real applications, the numerical integrations needed for the direct computation of the posterior are usually intractable. Briefly, VB uses an analytical approximation Q( ) to describe the true posterior, therefore shifting the problem to the minimization of the distance between the approximation Q( ) and the actual posterior distribution of the parameters, usually measured via the Kullback-Leibler divergence [28], that is equal to [29]: where F is the Free Energy which is defined as: Since log P( y) does not depend on and the Kullback-Leibler divergence is always non-negative, the minimization of the K-L divergence corresponds to maximizing the (negative) variational Free Energy F, that is also the lower bound to the log model evidence.
In this work we applied a mean field approximation of Q( ), as previously done in [18]: it consists in collecting the vector of parameters into two separate groups, one for the model parameters ( p) and one for the single noise parameter (φ). Each parameter has its own approximate posterior distribution (Q p ( p| y) and Q φ (φ| y)), assumed independent (i.e. Q( ) = Q p ( p| y) · Q φ (φ| y)), with φ −1 the unknown variance of the measurement error. We defined the prior distributions as a multivariate Gaussian distribution for the vector of model parameters p and a Gamma distribution for the variance scaling factor φ.
Furthermore, we used prior conjugated with the likelihood, i.e. with the same parametric form of the posterior.
where m and −1 are the mean vector and the covariance matrix of the multivariate normal distribution (with m 0 and −1 0 as corresponding prior mean and covariance) and s and c are the shape and scale parameter of the Gamma distribution (with s 0 and c 0 as corresponding prior values).
This simplifies the computation of the factorized posteriors, as the VB update becomes a process of updating the posterior hyper-parameters. The interested reader is referred to the original publications for a detailed derivation of the method and a full description of the parameter update and Free Energy equations [18].

D. Model Parameter Prior Setting
We developed a hierarchical approach for the empirical data-driven prior definition [19]. We used the estimates obtained from the model quantification at the region level with nonlinear least squares (NLLS) to inform the priors at the pixel level. The model used for the region-based analysis was the same subsequently used at the pixel level (Gammavariate or SCR). The regions of interest were defined within the synovia by a functional clustering (k-means partitioning method [30], Euclidean distance, 6 clusters) to detect the main kinetic behaviors of CEUS data as previously presented in [31]. The clustering was applied on the parametrized CEUS kinetics as previously presented in [32] and the parameters used to characterize the pixel concentration curves were the area under the curve, the time to peak and the slope of the initial and the last part of the curve (i.e. the first and the last 10 points). All these parameters were derived from the raw pixel time activity curves through linear interpolation. The resulting centroids were then fitted with the model under study solved with NLLS.
In particular, NLLS regional estimates were used as the multivariate Gaussian prior mean of the model parameters for all the pixels composing the selected region. As regard the prior variance 0 , we set it based on the region-wise NLLS estimates, multiplied by a constant factor λ, which gives a measure of the estimates variability across the synovia as 0 = 1/(λm 0 ) 2 . We set λ to 50%. This value was set based on previous studies [33].
We also incorporated the Automatic Relevance Determination (ARD) [20], [21] within the algorithm to allow for the automated reduction of model complexity (from SCR to Gamma-variate model). This was achieved by using a shrinkage prior for the parameter to be subjected to ARD, that in our case was the amplitude a 1 of the irreversible component. The prior distribution for a 1 was Gaussian with zero mean and unknown variance (set initially very large, i.e. 10 6 ) with a Gamma hyper-prior on the variance of a 1 to be estimated from the data. The ARD variance becomes very small if the data does not support the extra complexity introduced by this parameter, so that the a 1 is set to zero as the model collapses to the simple Gamma-variate.
In short, the VB estimation with the inclusion of ARD implies the quantification of CEUS data with a Gamma-variate model at the region level to determine the prior distributions for the other parameters (a 0 , α 0 , t 0 , β 0 ). Then the SCR model is solved at the pixel level, with the additional shrinkage prior that allows for the identification of a simpler model.
A summary of the models included in the analysis and their relationships along with the quantification methods is presented in Supplementary Figure 1. First, Bayesian model selection (BMS) was employed to compare the performance of Gamma-variate and SCR models to quantify CEUS kinetics. We used the Free Energy obtained from the VB algorithm as approximation of the log-model evidence (lower bound) [34]. We assumed equal prior probabilities for the models and we computed the log of the Bayes Factor, approximated as the difference of the Free Energies calculated for the two models in competition (SCR and Gamma-variate model). From the Bayes Factor it is possible to determine the strength of evidence (from weak to very strong) that one particular model generated the data for the population, following the scale of interpretation proposed by [35] (Table 1).
We performed the analysis on all pixels together regardless the subject they belong to, but excluding pixels for which at least one model failed to give physiological or reliable estimates (failure or outliers). Failures were defined as nonphysiological (null estimates) or unreliable values (coefficient of variation >200%) or pixels where the VB algorithm reached the maximum number of iterations without convergence. We also tested the randomness of the residuals, hypothesized to follow a standard normal distribution, with a Kolmogorov-Smirnov test.
Computational time was also taken into account in the model comparison. Both models were implemented with no parallelism using MATLAB Release 2015b, The MathWorks, Inc., Natick, Massachusetts, USA. The analysis was carried out on a quad core Intel Xeon E5-1620 Processor (3.50 GHz).
After assessing the presence and performance of the two models separately, we included ARD in the VB algorithm. We present the corresponding results in terms of failure percentage (as defined above), computational time, and model complexity ratio, i.e. the percentage of pixels where ARD identified a Gamma-variate or a SCR model.

F. Classification
In order to test the ability of the estimated parameters to distinguish different perfusion patterns, we trained a supervised classifier on the estimated parameter statistics (features). We trained two different classifiers in order to assess the discriminant power of the perfusion parameters obtained with Gamma and SCR model. In particular, we considered the 11 synovial perfusion parameters obtained with Gamma model solved with VB and the 12 synovial parameters obtained with SCR model solved with VB including ARD, plus the percentage of pixels where ARD detected a slow component, and the ratio a 1 /a 0 , as measure of the relative importance of the slow component on the CEUS kinetic. We also included for both classifiers the percentage of active pixels in the synovia. For each parameter distribution we considered several statistical descriptors (mean, standard deviation, 25th percentile and 75th percentile), to explore the importance of both the mean perfusion pattern in the synovia and its variability.
We sought to discriminate RA versus non-RA patients, in order to assess whether the identification of the apparent trapping can help distinguish this aggressive arthritis type. Subjects with less than 100 pixels were discarded from the analysis (12 non-RA and 2 RA patients), in order to maintain robust descriptors of the parameters for the classification. The dataset for the classification therefore comprised 54 RA and 40 non-RA patients.
Additionally, we built in the same way a classifier to discriminate RA and RA-like PSA from all other types of arthritis. In this case the dataset comprised 68 patients in the RA and RA-like PSA class, and 40 in the other arthritis class.
Since we have an unbalanced representation of the classes, the classifiers were validated using an hold-out scheme, by randomly selecting 70% of dataset as training set, but enforcing a balanced representation of the two classes in the training set (500 random samplings). Sensitivity, specificity and balanced accuracy (to take into account the unbalanced number of representatives for the two classes in the testing set) were used as performance indexes.
A backward feature selection was applied to discard noninformative features, and identify those more relevant for the classification.

A. Comparison of Gamma-Variate and SCR Model
Gamma-variate model solved with VB took 0.06 ± 0.02 seconds per pixel, for an average of 1.8 minutes analysis per subject. The median failure rate for Gammavariate model across subjects was 9% ± 14%, among which 6% of outliers were due to pixels where β 0 was estimated unrealistically high (>1000), indicating a very slow washout of the kinetics. The computational time required by VB to solve a more complex model such SCR did not increase significantly (elapsed time, 0.07 ± 0.03 seconds per pixel, for an average of 1.8 minutes per subject). The median failure rate was slightly higher for SCR model (16% ± 21%) mostly due to pixels were the algorithm did not converge to a reliable solution.
When we compared Gamma-variate and SCR model performance in the group analysis (i.e. considering all the pixels together), there was positive evidence for SCR to be the model of choice in 51% of pixels (Fig. 1 A). The Gammavariate model probability was important as well, with positive evidence in 44% of pixels. This result held also at the single subject level, where all possible combinations were found: in some patients, Gamma-variate model was more probable than SCR model, in others it was far less probable than SCR model, in other cases it was equally probable, depending on the subject considered ( Fig. 1 B-D).
These results supported the implementation of ARD, to automatically infer the model complexity from the data.

B. On-Line Model Selection With ARD
VB including ARD took 0.08 ± 0.03 seconds per pixel (median ± median absolute deviation, MAD), for an average of 2.2 minutes analysis per subject. The median failure rate across subjects for VB with ARD was 17% ± 17%.
When we applied VB including ARD to describe CEUS kinetics, the algorithm selected the SCR model in 47% of pixels, while the Gamma-variate was chosen in the other 53%, confirming the results of the Bayesian model selection.
When we analyzed separately the different disease types in a group analysis, the algorithm identified Gamma-variate and SCR model almost in the same percentage of pixels for RA, dpPSA and CTD. Interestingly, Gamma-variate was identified by ARD as the model of choice in 86% of pixels for uSPA and SCR was the best model in 61% of pixels in RA-like patients. Despite the slow component was found thoroughly in all patients, its amplitude changed significantly across different types of arthritis. Mean values of a 1 and its variability were statistically higher in RA patients (168% increase in mean, p = 0.002 and 140% increase in standard deviation, p = 0.001 respectively).
An example of the utility of the application of VB with ARD is presented in Fig. 2, where CEUS kinetics are best described by the SCR model at the region level (panel A), according to all parsimony criteria considered (BMS, AIC, BIC) [36], [37]. Within the same cluster, there is a variety of kinetics, some best described by Gamma-variate model (an example is presented in panel B), others by SCR model (panel C).
The information on model complexity returned by ARD can also be used to visually confirm the kinetic heterogeneity of CEUS data in the synovia, as can be clearly seen by the cluster analysis results to derive the regions of interest (Fig. 3 A). By considering the parametric map of ARD results, i.e. the map representing for each pixel whether the best model was the Gamma-variate or SCR one, we found that in general the Gamma-variate model is identified at the boundaries of the synovia (Fig. 3 B). Moreover, it is interesting to note that the model complexity does not overlap with the cluster boundaries, thus it is not dependent on the prior definition. Another important aspect is related to the amplitude of the trapping component, a 1 : the spatial distribution of this parameter reflects an important heterogeneity in the synovia (Fig. 3 C).

C. Classification
The SVM classifier trained to detect RA patients versus non-RA patients with synovial perfusion parameters obtained with VB and ARD (SV M AR D ) achieved a balanced accuracy of 87%, with a sensitivity of 88% and a specificity of 86%, with 24 features selected. The reduced set of 24 features included the percentage of pixels where ARD selected SCR On the contrary, the SVM classifier trained with perfusion estimates obtained from Gamma model (SV M Gamma ) achieved a balanced accuracy of 81%, a specificity of 94% but a sensitivity of only 68%, with 12 features selected. The reduced set of 12 features included the mean and variability of the wash-out time and the blood flow index. When setting the operating point of the SV M AR D classifer to 81% accuracy to match that of SV M Gamma , only 7 perfusion features where necessary, and those selected included the mean of a 1 and the standard deviation and 75 th percentile of a 1 /a 0 .
When evaluating the SVM classifier to distinguish between RA and RA-like PSA from the other forms of arthritis, the one making use of the perfusion parameters obtained with the SCR model and ARD selection achieved 100% balanced accuracy, versus 83% of the SVM trained on the parameters obtained by the Gamma model.

A. Gamma-Variate and SCR Model
Gamma-variate model has been shown to be a more appropriate model to describe the CEUS perfusion curves compared to the simple mono-exponential or logarithmic model [9], [38]  and it has been applied both at the region and pixel level [4], [8], [39], [40]. In fact, even though the mono-exponential can accurately describe CEUS kinetics at the region level, it lacks the flexibility to proper follow the pixel-wise CEUS time courses (Supplementary Figure 2 A), especially when the data present an important wash-out ( Supplementary  Figure 2 B). However, the heterogeneity of CEUS kinetics can be such that even Gamma-variate model is unable to fully describe the behaviour of the tracer during the scan time length. For this reason, we propose the use of the SCR model to describe CEUS data in order to take into account the irreversible component in the data, generally considered negligible. The SCR model consists in a single Gamma-variate plus an additional component (the integral of the Gamma), originally proposed to describe the recirculation of the contrast within the vasculature. The integral component of the SCR model could also describe the microbubbles behavior in the neovascularized synovial area in arthritis disease. In fact newly formed vessels are characterized by high tortuosity [13], [14] and the bulk of proliferated vascular network, by increasing the resistance to synovial vessel flow, might determine the prolongation of the refilling time [5]. A reduction of blood volume fraction of synovial capillaries has been already observed in rheumatoid arthritis [41], [42]. Our results are in agreement with these previous findings: in RA patients, the mean and variability of the amplitude of the slow component detected in the data is significantly higher than in patients affected by other, less aggressive, forms of arthritis.
Among the potential factors which may affect the observed synovial refilling time, there are the vessel intrinsic factors, as tortuosity of the microvasculature, reduced vessel diameter, or capillary obliteration [43]. In this case, the passage of the microbubbles in the vessels might slow down so that it is not possible to identify the washout within the experimental time of the CEUS exam, thus appearing as a trapping.
Another possible explanation could be the increase in the intra-articular pressure due to the increase in synovial fluid in chronic arthritis [44]. However, this increase would imply a more homogeneous trapping effect in the synovia, which is in conflict with the spatial gradient observed in the amplitude of the integral component.
When applied to CEUS data, SCR resulted the most probable model in 51% of pixels, highlighting the necessity to consider the trapping component in at least half of the data. However, the single Gamma-variate was the model of choice with positive evidence in the other 44% of pixels (there was a 5% of pixels with weak evidence for both models and nothing can be said about the optimal model in that subset). In principle, one should fit both models for each pixel and then choose the best one based on some parsimony criterion. More efficiently, since the two models are nested, we employed the automatic relevance determination, originally presented in [20], to automatically infer the model complexity directly from the data. The results obtained applying VB including ARD confirmed the initial results: the algorithm selected the SCR model in 47% of the total pixels, while the simpler Gamma-variate was sufficient to describe the CEUS data in the other 53% of pixels. The inclusion of the ARD in the Variational Bayesian method did not increase the computational time requested to complete a single subject analysis (2.2 minutes on average). However, the model on-line selection comes at a price, namely a slight increased percentage of outliers (between subject mean: 17% ± 17%). Interestingly, the model complexity follows a specific trend within the synovial area, being spatially dependent. In most subjects, the irreversible component is not necessary at the boundaries of the synovia, where the single Gamma-variate model is selected (Fig. 3).
When we used the perfusion estimates obtained with VB including ARD to train a SVM classifier, we were able to differentiate between RA versus other forms of arthritis with a sensitivity of 88%, a specificity of 86% and an overall balanced accuracy of 87%. Although this level of accuracy and specificity might not be useful in terms of straightforward clinical diagnosis, it is important from a translational point of view. Being able to differentiate more aggressive types of arthritis on the base of perfusion patterns could help not only to unveil the etiopathology of arthritis but also to better depict treatment response. For example, the use of a more complex model allows to find significant differences between patients only on the base of the slow component, that also plays a predominant role in the feature selection. The SCR model has never been proposed for the analysis of CEUS data, and in this sense any new information can be extremely important for the study of the pathology.
It is interesting to note that these results are not affected by the orientation of blood vessels relative to the imaging plane. The locality of the analysis ensures that the timeintensity curve is related to the appearance of a quantity of microbubbles in a imaged region (pixel), with very little influence from the direction (whether the vessel is longitudinal or cross-sectional to the plane) or speed of the blood flow.
CEUS has great potential for many other clinical applications (with data quantified with SCR and solved with VB), where it is important to assess the microvasculature structure. The most prominent examples are the study of tumours (e.g. to assess tumour angiogenesis or for lesion detection in brain, liver or kidney), assessment of myocardial perfusion and vascular pathologies (for example, aneurysms or dissections). In fact, CEUS great advantage with respect to other imaging modalities is the ability to image small capillaries (mean diameter 8-9 μm, [45]): the only possible limitation is the ultrasonic penetration (and echo capture) within the tissue to be studied during the exam (e.g. in arthritis, the presence of the bone that limits the visibility of the whole articulation).

B. Definition of the Bayesian Priors
The main issue in the implementation of a Bayesian estimator is the definition of an appropriate set of priors. The prior distribution of the model parameters can be set to literature values, however CEUS data heterogeneity prevents the use of strict prior information for the totality of the pixels in the synovia. Moreover, the presence of different types of arthritis, each characterized by a different morphogenesis and vascular patterns, precludes the definition of general priors.
A valid alternative is to retrieve robust and reliable information directly from the data, in a user-independent top-down approach: following a hierarchical scheme the prior distributions to use at the pixel level are obtained for each subject from the region-wise model estimates solved with NLLS. In our case, regions of interest are derived by functional parametric clustering, with k-means algorithm and 6 clusters (which were reduced to 2 clusters when the synovia comprised less than 100 pixels). The number of clusters was set empirically to detect the different main components in the data [31]. In general, the number of clusters determines in a way the quality of the prior distributions: the use of too many clusters will result in detailed segmentation of the synovia and therefore noisy region-wise CEUS kinetics that could be lead to difficulties in the numerical estimation of the curves. On the other hand, the use of too few clusters will result in a coarse segmentation of the synovia leading to less informative priors (due to mixed kinetic components in the same ROI).
In this way we can determine tailored priors for each region of interest, with the prior variance set based on the prior mean, with a 50% of coefficient of variation. This value was set based on previous studies [33].
It is important to note that despite the definition of a common prior for all the pixels in the region, variation in the parameter value at the pixel level is permitted in the inference procedure.

C. Variational Bayesian Performance
VB is a Bayesian method that has been shown to outperform conventional NLLS by providing accurate and robust estimates, with a substantial lower percentage of outliers, when applied to both magnetic resonance [18].
To further confirm its validity also to quantify CEUS data, we compared the results obtained with VB with those obtained solving the model at the pixel level with NLLS, generally considered as the reference method for CEUS quantification at the region and pixel level. We analyzed separately the Gammavariate and the SCR model. VB and NLLS estimates were compared in the intersection of pixels where both methods gave reliable estimates (as detailed in the methods section). We compared outliers percentage, correlation (as Pearson's correlation coefficient R, slope and intercept of the regression line) and mean relative difference (MRD) between VB and NLLS estimates for both micro-and macroparameters of interest.
The main impact of the method was the percentage of outliers: for Gamma-variate model NLLS did not converge in 22% ± 24% of the pixels while for SCR model NLLS outliers raised to 28% ± 25%.
The results obtained with NLLS and VB (where both methods converged to reliable solutions) were all highly correlated with small relative difference ( Table 2, Fig. 4). Gamma-variate model showed excellent agreement for macroparameters and good agreement for microparameters (as expected). The same held for SCR model, albeit the variability in the mean relative difference was greater, due to the highest model complexity (Fig. 4).
VB outperformed NLLS also in terms of computational time: for a Gamma-variate model NLLS took an average of 0.08 ± 0.02 seconds per pixel, while for SCR NLLS took 0.10 ± 0.03 seconds per pixel.   VB and NLLS performances in terms of computational time and outliers percentage are summarized in Table 3.

V. CONCLUSIONS
The inclusion of an irreversible component better describes CEUS data in arthritis disease in half of the pixels. The heterogeneity of the data however makes it fundamental the implementation of a method that automatically infer the model complexity from the data. Variational Bayes with automatic relevance determination is applied for the first time to CEUS data and it provides robust and accurate model estimates with low percentage of outliers, recovering physiological information in the majority of the pixels. Computational time required for a single subject analysis is less than 3 minutes and therefore compatible with clinical practice.

APPENDIX A DERIVATION OF MACROPARAMETERS
For both Gamma-variate and SCR model is possible to determine a set of 7 macroparameters of interest, namely the peak value y max , the time of peak t max , the raise time t raise and the washout time t wash , the mean transit time MT T , the blood flow index B F I and the blood volume index BV I . We report here the equations for the calculation of these parameters.
since a 0 has been normalized by the value of peak. The parameter t raise and t wash cannot be calculated in closed form, since often there is no washout in the Gammavariate model: where L is the Lambert W function. Interestingly, the formulas for the calculation of the macroparameters for SCR model do not change, since these are calculated on the Gamma-variate course and the trapping component is excluded from their calculation.