## Abstract

Ultrasound contrast agents (UCAs) have opened up immense diagnostic possibilities by combined use of indicator dilution principles and dynamic contrast-enhanced ultrasound (DCE-US) imaging. UCAs are microbubbles encapsulated in a biocompatible shell. With a rheology comparable to that of red blood cells, UCAs provide an intravascular indicator for functional imaging of the (micro)vasculature by quantitative DCE-US. Several models of the UCA intravascular kinetics have been proposed to provide functional quantitative maps, aiding diagnosis of different pathological conditions. This article is a comprehensive review of the available methods for quantitative DCE-US imaging based on temporal, spatial and spatiotemporal analysis of the UCA kinetics. The recent introduction of novel UCAs that are targeted to specific vascular receptors has advanced DCE-US to a molecular imaging modality. In parallel, new kinetic models of increased complexity have been developed. The extraction of multiple quantitative maps, reflecting complementary variables of the underlying physiological processes, requires an integrative approach to their interpretation. A probabilistic framework based on emerging machine-learning methods represents nowadays the ultimate approach, improving the diagnostic accuracy of DCE-US imaging by optimal combination of the extracted complementary information. The current value and future perspective of all these advances are critically discussed.

## Key Words

## Introduction

The value of dynamic contrast-enhanced ultrasound (DCE-US) for physiological quantification is today well established in many clinical applications. DCE-US relies on the intravascular injection of ultrasound contrast agents (UCAs) consisting of gas-filled microbubbles encapsulated in a biocompatible shell (

Correas et al., 2001

; Cosgrove, 2006

; Quaia, 2007

; Mischi et al., 2018

). These act as echo-enhancers and generate a strong acoustic response when insonified around their resonance frequency (Correas et al., 2001

; Sboros, 2008

; Mischi et al., 2018

). This acoustic response, highly non-linear, is exploited by contrast-specific ultrasound imaging modes, which, by suppressing linear tissue echoes, considerably enhance the visualization of echoes originated from microbubbles (Frinking et al., 2000

; Bouakaz et al., 2002

). The echogenicity of UCAs is characterized by the backscatter coefficient, which can be modeled as a function of the total scattering cross-section of the microbubbles and their concentration (Bouakaz et al., 1998

; Frinking and de Jong, 1998

). The latter paves the way for quantitative imaging by application of the indicator dilution theory.To achieve interpretable quantification, the rheology of UCAs has been investigated and characterized (

Ismail et al., 1996

; Fisher et al., 2002

; Chatterjee and Sarkar, 2003

; Leen et al., 2012

). Indeed, with a diameter of 1–10 μm, the rheology of microbubbles is similar to that of red blood cells (Ismail et al., 1996

; Lindner et al., 1998

; Gorce et al., 2000

; Fisher et al., 2002

; Quaia, 2011

; Leen et al., 2012

; Turco et al., 2016

); comparable distributions and velocities have been observed down to the capillary level (Leen et al., 2012

).In the past few years, the introduction of targeted UCAs has advanced DCE-US to a molecular imaging modality (

Abou-Elkacem et al., 2015

; Smeenge et al., 2017

; Willmann et al., 2017

). These novel agents are obtained from conventional UCAs by decorating the microbubble shell with targeting ligands to bind to specific receptors expressed on the vascular endothelium; this produces selective enhancement in areas where these receptors are overexpressed (Abou-Elkacem et al., 2015

). More recently, the development of nanobubbles able to cross the vascular endothelium is opening up new possibilities for assessment of vascular permeability, a marker of several pathological conditions (Maeda et al., 2000

; Nagy et al., 2008

), and for molecular imaging of targets outside the blood-pool compartment (Perera et al., 2015

).After UCA injection, quantification is performed by analysis of time–intensity curves (TICs) extracted from DCE-US loops, which describe the temporal evolution of the signal enhancement (acoustic intensity) caused by UCA transport through the investigated organ or tissue (Fig. 1). Indicator dilution theory and pharmacokinetic modeling can be used for analysis of TICs extracted from DCE-US, enabling the estimation of quantitative parameters related to, for example, blood volume, flow, perfusion, extravasation and molecular binding (

Strouthos et al., 2010

; Turco et al., 2016

; Mischi et al., 2018

).The use of indicators for physiological quantification has, in fact, a long history. The first studies from Haller (1761), on the transpulmonary transit time (

Fox, 1962

), and Hering (1827), on the total circulating blood volume (Tigerstedt, 1921

), were extended by Stewart (1897) to measure the cardiac output (CO), and corrected by Hamilton (1931) to take recirculation of the indicator into account (Stewart, 1921a

; Fox, 1962

). In its early stage, indicator dilution theory was used mainly for cardiovascular quantification and required invasive sampling of the indicator (usually dye) concentration by central catheterization (Stewart, 1921b

; - Stewart G.N.

Researches on the circulation time and on the influences which affect it: V. The circulation time of the spleen, kidney, intestine, heart (coronary circulation) and retina, with some further observations on the time of the lesser circulation.

*Am J Physiol Legacy Content.*1921; 58: 278-295

Fox, 1962

)In parallel, the field of pharmacokinetic modeling was born in the early 1990s (

Wagner, 1981

), with the introduction of the first compartmental model by Widmark and Tandberg, 1924

. A few years later, Teorell, 1937

presented the first multicompartment model to describe drug transport and clearance through the kidneys, which is considered today the foundation of modern pharmacokinetic modeling (Himmelstein and Lutz, 1979

; Wagner, 1981

). Since then, pharmacokinetic modeling has found widespread applications for assessment of drug transport and deposition in “non-accessible” tissues, based on multiple measurements performed in “accessible” compartments by, for example, blood sampling (Himmelstein and Lutz, 1979

; Wagner, 1981

).The advent of imaging, and especially contrast-enhanced imaging, has widened the spectrum of possibilities for minimally invasive physiological quantification. In DCE-US, combining the injection of UCAs with contrast-specific scanning sequences has enabled

*in vivo*sampling of indicators (contrast agents) at unprecedented temporal resolution also in sites that were previously inaccessible, with the addition of the spatial dimension determined by the ultrasound field of view. Although indicator dilution and pharmacokinetics studies were previously limited to whole-body or system-level analysis, novel imaging techniques are now permitting spatiotemporal analysis of contrast kinetics down to suborgan and molecular levels (Mischi et al., 2012

).In the past decades, several methods for quantitative analysis of DCE-US loops have been developed and established, evidencing clinical value in many areas, especially in cardiology and oncology. Quantification of myocardial perfusion was the first application (

Wei et al., 1998

), followed by the assessment of the angiogenic microvasculature associated with cancer and inflammatory processes. Angiogenesis, that is, the formation of a chaotic network of blood vessels, is known to be involved in several pathological processes. Through leveraging of their intravascular character, UCAs have been successfully employed for non-invasive assessment of abnormal angiogenic vasculature in many pathological conditions, including unstable atherosclerotic plaque, diabetes, hyperlipidemia, placental insufficiency, inflammation and cancer (Lindner, 2004

; Arthuis et al., 2013

; Roberts et al., 2016

; Turco et al., 2016

; Herbst et al., 2017

; Mischi et al., 2018

). Besides diagnostic applications, DCE-US has shown promise for therapy monitoring and follow-up (Lamuraglia et al., 2010

; Ryu et al., 2013

) and for mediating drug delivery (Lindner, 2004

; Sennoga et al., 2017

).Starting from hepatic applications, different guidelines and recommendations for clinical use of DCE-US have been issued over the years, evidencing the broad range of established and emerging applications (

Claudon et al., 2013

; - Claudon M.
- Dietrich C.F.
- Choi B.I.
- Cosgrove D.O.
- Kudo M.
- Nolsøe C.P.
- Piscaglia F.
- Wilson S.R.
- Barr R.G.
- Chammas M.C.
- Chaubal N.G.
- Chen M.H.
- Clevert D.A.
- Correas J.M.
- Ding H.
- Forsberg F.
- Fowlkes J.B.
- Gibson R.N.
- Goldberg B.B.
- Lassau N.
- Leen E.L.
- Mattrey R.F.
- Moriyasu F.
- Solbiati L.
- Weskott H.P.
- Xu H.X.

Guidelines and good clinical practice recommendations for contrast enhanced ultrasound (CEUS) in the liver—Update 2012: A WFUMB–EFSUMB initiative in cooperation with representatives of AFSUMB, AIUM, ASUM, FLAUS and ICUS.

*Ultraschall Med.*2013; 34: 11-29

Cantisani et al., 2015

; Nolsøe and Lorentzen, 2016

; Dietrich et al., 2018

; - Dietrich C.F.
- Averkiou M.
- Nielsen M.B.
- Barr R.G.
- Burns P.N.
- Calliada F.
- Cantisani V.
- Choi B.
- Chammas M.C.
- Clevert D.A.
- Claudon M.
- Correas J.M.
- Cui X.W.
- Cosgrove D.
- D'Onofrio M.
- Dong Y.
- Eisenbrey J.
- Fontanilla T.
- Gilja O.H.
- Ignee A.
- Jenssen C.
- Kono Y.
- Kudo M.
- Lassau N.
- Lyshchik A.
- Franca Meloni M.
- Moriyasu F.
- Nolsøe C.
- Piscaglia F.
- Radzina M.
- Saftoiu A.
- Sidhu P.S.
- Sporea I.
- Schreiber-Dietrich D.
- Sirlin C.B.
- Stanczak M.
- Weskott H.P.
- Wilson S.R.
- Willmann J.K.
- Kim T.K.
- Jang H.J.
- Vezeridis A.
- Westerway S.

How to perform contrast-enhanced ultrasound (CEUS).

*Ultrasound Int Open.*2018; 4: E2-E15

Sidhu et al., 2018

). In addition to promoting the standardization of UCA injection and acquisition protocols, and of qualitative interpretation of DCE-US exams, important steps have been made in quantitative analysis, with the publication of consensus papers providing guidelines for quantification of tumor perfusion by DCE-US (- Sidhu P.S.
- Cantisani V.
- Dietrich C.F.
- Gilja O.H.
- Saftoiu A.
- Bartels E.
- Bertolotto M.
- Calliada F.
- Clevert D.A.
- Cosgrove D.
- Deganello A.
- D'Onofrio M.
- Drudi F.M.
- Freeman S.
- Harvey C.
- Jenssen C.
- Jung E.M.
- Klauser A.S.
- Lassau N.
- Meloni M.F.
- Leen E.
- Nicolau C.
- Nolsoe C.
- Piscaglia F.
- Prada F.
- Prosch H.
- Radzina M.
- Savelli L.
- Weskott H.P.
- Wijkstra H.

The EFSUMB guidelines and recommendations for the clinical practice of contrast-enhanced ultrasound (CEUS) in non-hepatic applications: Update 2017 (short version).

*Ultraschall Med.*2018; 39: 154-180

Dietrich et al., 2012

; Leen et al., 2012

). Moreover, important steps to combine quantitative parameters extracted from DCE-US and other ultrasound modalities in a multiparametric fashion are also being made (Acharya et al., 2011

; - Acharya U.R.
- Faust O.
- Sree S.V.
- Molinari F.
- Garberoglio R.
- Suri J.S.

Cost-effective and non-invasive automated benign & malignant thyroid lesion classification in 3D contrast-enhanced ultrasound using combination of wavelets and textures: A class of ThyroScan algorithms.

*Technol Cancer Res Treat.*2011; 10: 371-380

Brock et al., 2013

; Postema et al., 2015

; Mannaerts et al., 2018

, Mannaerts et al., 2019

; - Mannaerts C.K.
- Wildeboer R.R.
- Remmers S.
- van Kollenburg R.A.A.
- Kajtazovic A.
- Hagemann J.
- Postema A.W.
- van Sloun R.J.G.
- Roobol M.
- Tilki D.
- Mischi M.
- Wijkstra H.
- Salomon G.

Multiparametric ultrasound for prostate cancer detection and localization: Correlation of B-mode, shearwave elastography and contrast-enhanced ultrasound with radical prostatectomy specimens.

*J Urol.*2019; 202: 1166-1173

Theek et al., 2018

).In this article, we provide a comprehensive review of quantification methods that have been developed and used in DCE-US. The aim is to provide the reader with the theoretical basis for proper selection and use of pharmacokinetic models in DCE-US, with special focus on the physiological interpretation of the obtained results. To this end, recommendations are also provided for acquisition and (pre)processing of DCE-US data. After introducing the principles of the indicator dilution theory and discussing the relevance and implications of signal acquisition and calibration, we provide a review of the quantification methods for temporal, spatial, spatio-temporal and molecular analysis of UCA kinetics. We then discuss the emerging role of machine learning toward multiparametric DCE-US. Finally, in addition to our recommendations for quantitative DCE-US, upcoming possibilities and future perspectives are also discussed.

## Principles of indicator dilution theory

DCE-US quantification requires the intravenous administration of UCAs, followed by acquisition and analysis of the UCA time evolution in a certain region of interest (ROI). The latter is properly positioned on the investigated tissue or organ. UCA administration can follow two main protocols: continuous infusion or bolus injection. The first requires the employment of an infusion pump, and the second is usually performed manually. Although continuous infusion is aimed at establishing a steady plateau concentration in the region of interest, bolus injection is aimed at mimicking a mathematical (Dirac) impulse function, which represents a fundamental assumption of most kinetic models. The latter can be approximated by executing a fast bolus injection, but the associated pressure gradients can cause bubble collapse and affect the accuracy and reliability of the measurement.

Common UCAs, such as SonoVue (Bracco Imaging S.p.A., Colleretto Giacosa, Italy) and Definity (Lantheus, North Billerica, MA, USA), are prepared before the injection by reconstitution of a lyophilized powder or by emulsification activation using a mechanical shaking device, respectively. When multiple injections are performed, the agent should be reactivated by gentle shaking prior to the injection (

Sugimoto et al., 2012

). For optimal infusion, dedicated pumps have been developed for continuous mixing of the agent (Lohmaier et al., 2004

).Indicator dilution theory is based on a fundamental concept: if the concentration of an indicator that is uniformly dispersed in an unknown volume

where the dose

*V*is determined, and the amount of the indicator (injected dose) is known, the unknown volume can then be determined. This requires instantaneous and uniform mixing of the agent. After the injection at time*t*= 0 of a UCA dose,*m*, into a dilution system with constant carrier (volumetric) flow,*F*, and UCA concentration,*C*(*t*), application of the mass conservation law leads to$F=\frac{m}{\underset{0}{\overset{\infty}{\int}}C\left(t\right)dt},$

(1)

where the dose

*m*can be interpreted as the number of injected microbubbles, and the concentration*C*(*t*) as the number of microbubbles per unit volume (Mischi et al., 2003

). Introduced by the pioneers of indicator dilution theory, eqn (1) is referred to as the Stewart–Hamilton equation (Stewart, 1921a

; Hamilton et al., 1928

).As the circulatory system is a closed system, recirculation of the indicator occurs. Hamilton extended the application of eqn (1) to the presence of recirculation, extracting the first pass of the indicator by an exponential fit of the early downslope (

Zierler, 2000

; Stewart, 1921a

; Hamilton et al., 1928

). This avoided the underestimation of the flow *F*in eqn (1) caused by overestimation of the integral of*C*(*t*). Most importantly, the exponential model corresponds to a mono-compartment model (Himmelstein and Lutz, 1979

; Wagner, 1981

), represented by a mixing chamber of volume *V*with flow*F*going from the inlet to the outlet of the chamber, as illustrated in Figure 2.With

When a dose

with τ

*C*(*t*) being the indicator concentration in the compartment, and*C*_{in}(*t*) being the concentration entering the compartment via the inlet, the mass conservation law can be formulated as$\frac{dC\left(t\right)}{dt}V=\left({C}_{\text{in}}\left(t\right)-C\left(t\right)\right)F.$

(2)

When a dose

*m*=*VC*_{0}is administered in the compartment at time*t*= 0 as a Dirac function,*C*_{in}(*t*) =*C*_{0}*u*_{0}(*t*), the compartmental concentration is given as$C\left(t\right)={C}_{0}{e}^{-\frac{t}{\tau}},$

(3)

with τ

*= V/F*the time constant of the compartment. Equation (3) corresponds to the exponential decay introduced byHamilton et al., 1928

.This simple model has found several applications. By interpreting the left ventricle (LV) as a single compartment, Holt et al. proposed a method for assessment of the LV ejection fraction (EF) based on the wash-out curve of an indicator (Evans blue dye) injected directly in the LV cavity through central catheterization (

Holt, 1956

). With *Cn*and*Cn*_{+1}being the concentration at two subsequent systoles, EF can be derived as$\text{EF}=1-\frac{{C}_{\mathrm{n}+1}}{{C}_{n}}=1-{e}^{-\frac{\Delta t}{\tau}}.$

(4)

The last term is obtained by modeling the LV wash-out as in eqn (3), with Δ

*t*being the average cardiac period during wash-out (Rovai et al., 1987

; Mischi et al., 2005

). This improves the reliability of the estimate by including multiple cardiac cycles.The mono-compartment model in eqn (2) also allows for assessment of myocardial perfusion. Assuming the myocardium to be represented by a mono-compartment dilution model,

Wei et al., 1998

proposed the assessment of myocardial perfusion by the so-called *destruction–replenishment technique*; under a UCA infusion protocol, UCA replenishment is imaged after UCA destruction by a high-pressure ultrasound burst. The obtained replenishment curve (Fig. 3) can then be interpreted as the step response of a mono-compartment dilution system to assess perfusion as*C*_{0}/*τ*, with*C*_{0}being the plateau UCA concentration prior to the destruction burst.A key asset provided by the use of indicators consists of the estimation of volumes. After the pioneering work of

Here,

Stewart, 1921a

on the transpulmonary circulation time and blood volume, Meier and Zierler, 1954

provided an analytical formulation of the problem based on the concept of mean transit time (MTT). Following a fast bolus injection at time *t*= 0 through the inlet of a dilution system, given as a Dirac function,*C*_{in}(*t*)=*C*_{0}*u*_{0}(*t*), the MTT between injection and detection sites can be estimated by the first statistical moment of the detected dilution curve,*C*(*t*). After normalization by the dilution curve integral, referred to as area under the curve (AUC), the MTT can be derived as$\text{MTT}=\frac{\underset{0}{\overset{\infty}{\int}}tC\left(t\right)dt}{\underset{0}{\overset{\infty}{\int}}C\left(t\right)dt}=\frac{\underset{0}{\overset{\infty}{\int}}tC\left(t\right)dt}{\text{AUC}}=\underset{0}{\overset{\infty}{\int}}th\left(t\right)dt.$

(5)

Here,

*h*(*t*) represents the fraction of injected indicator leaving the system per unit time. It corresponds to the probability density function (PDF) of UCA transit times from the injection to the detection site (Mischi et al., 2004

; Kuenen et al., 2013a

). Indeed, a statistical interpretation of the normalized indicator dilution curve is typically provided by a PDF of transit times.Once the MTT is assessed, the volume between injection and detection sites can be derived as the product between flow and MTT, provided that the indicator is injected and detected at the only inlet and outlet of the dilution system. Therefore, the MTT (inverse of the transit rate) can also be interpreted as the ratio between distribution volume and flow. This condition applies to transpulmonary measurements of the indicator concentration on the right (inlet) and the left (outlet) side of the heart. Multiplication of the MTT by the CO provides an estimate of the pulmonary blood volume (PBV) or the intrathoracic blood volume (ITBV, including the blood volume in the cardiac chambers), depending on the position of the injection and detection sites (

Sakka et al., 2000

).Based on eqns (3) and (5), the MTT of a mono-compartment model corresponds to its time constant τ.

Newman et al., 1951

interpreted the time constant of the exponential decay of a measured indicator dilution curve as representing the largest volume in a cascade of compartments. On the basis of this result, the PBV can be separated from the ITBV even when the measurement includes all cardiac chambers. This assumption was used by Sakka et al., 2000

to assess the transpulmonary circulation by thermodilution. The use of DCE-US allows this assumption to be overcome by positioning of the regions of interest for measuring the dilution curves at the input and output of the transpulmonary dilution system. Multiplication of the CO by the difference in MTT between the input and output dilution curves provides a measure of PBV (Herold et al., 2015

). *In vitro*measurements have proven the DCE-US method superior to thermodilution (Herold et al., 2013

).In general, any normalized dilution curve modeling the output of a dilution system assuming a Dirac input can also be interpreted as the impulse response of the system, assumed to be linear and time invariant. The output curve is then interpreted as the result of a convolution between any input function and the impulse response of the system. This leads to several relevant applications; when different indicator dilution curves represent the input and output of a dilution system, this can be analyzed irrespective of the actual indicator injection site by a

*deconvolution*operator (Mischi et al., 2007

).The impulse response of a mono-compartment model is given as

and the concentration

$h\left(t\right)=\frac{1}{\tau}{e}^{-\phantom{\rule{0.25em}{0ex}}{\textstyle \frac{t}{\tau}}},$

(6)

and the concentration

*C*(*t*) in the compartment can be derived as the convolution $C\left(t\right)=h\left(t\right)*{C}_{\text{in}}\left(t\right)$. By use of a deconvolution operator in combination with DCE-US, the invasive procedure for EF measurement based on eqn (4) can be replaced by a simple peripheral intravenous injection. Based on the dilution curves measured in the left atrium and LV, the LV impulse response can be estimated by a deconvolution approach, leading to a minimally invasive estimation of the EF (Mischi et al., 2005

, Mischi et al., 2007

).A similar approach can be used for estimation of pulmonary circulation time and PBV without central catheterization. In fact,

*h*(*t*) in eqn (5) corresponds to the impulse response of the transpulmonary dilution system. Dilution curves measured in the right and left sides of the heart (before and after the lungs) can be processed by deconvolution to estimate the impulse response representing the transpulmonary circulation, which can then be analyzed to extract relevant physiological parameters (Mischi et al., 2007

, Mischi et al., 2009

; Saporito et al., 2016

).- Saporito S.
- Herold I.H.
- Houthuizen P.
- van Den Bosch H.C.
- Den Boer J.A.
- Korsten H.H.
- van Assen H.C.
- Mischi M.

Model-based characterization of the transpulmonary circulation by dynamic contrast-enhanced magnetic resonance imaging in heart failure and healthy volunteers.

*Invest Radiol.*2016; 51: 720-727

## Acquisition and calibration

Application of the indicator dilution theory is based on the indicator concentration,

*C*(*t*). Therefore, DCE-US quantification, whether based on bolus kinetics or replenishment kinetics following bubble destruction, requires a direct relationship between echo signal amplitude and contrast-agent microbubble concentration in a way that is independent of system or user settings. For contrast-agent microbubbles, the instantaneous echo power is a quantity proportional to the local concentration (Arditi et al., 2006

; Lampaskis and Averkiou, 2010

; Payen et al., 2013

), as long as their backscattered ultrasound energy originates from an ensemble of randomly spaced and statistically uncorrelated microbubbles (de Jong and Hoff, 1993

; Uhlendorf and Hoffmann, 1994

). This relation is described by the backscatter coefficient, accounting for the scattering cross-section of all the microbubbles (Bouakaz et al., 1998

; Frinking and de Jong, 1998

). Unfortunately, quantification of the acoustic backscatter is complicated in clinical practice by the effect of microbubble size distribution, attenuation, dynamic range and implemented post-processing (Rovai et al., 1987

; Mischi et al., 2003

).Any method or software tool dedicated to quantification should ideally operate on signals proportional to instantaneous echo power, which can be obtained by signal linearization. When raw echo signals are provided, that is, radiofrequency or demodulated signals before scan conversion, linearization can simply be obtained by squaring these raw signals (

Rognin et al., 2008

; Payen et al., 2013

). In most ultrasound systems, however, raw signals are not available, and clinical users only have access to logarithmically compressed video data. A practical approach here consists then of collecting scan-converted video data, log-compressed and palettized as gray-scale or color-coded 8-bit data, in the form of DICOM files (Rognin et al., 2008

; Payen et al., 2013

). Proper linearization consists then of reversing the effects of log compression, and possibly non-linear palette rendering, followed by squaring. As mentioned before, this is necessary to obtain reliable estimates of relative echo power in the images, representing local microbubble concentration.In both cases of data input, that is, either raw or video data, linearization re-establishes proportionality between signal amplitude, or echo power, and UCA concentration (

Arditi et al., 2006

; Lampaskis and Averkiou, 2010

; Payen et al., 2013

). As a result, data obtained with different ultrasound systems, between departments within a hospital or in multicentric studies can be compared, though only in a relative sense. This is nicely illustrated by Payen et al., 2013

, who reported direct proportionality between echo power and local contrast agent concentration over more than three orders of magnitude in concentration for both types of data (Fig. 4). Note, however, the decibel scale on the ordinate axis indicating the relative nature of the different data and that overlapping of the curves (i.e., system/method independence) is obtained only after appropriate scaling.Thus, processing log-compressed video data without appropriate linearization induces quantification results, which are dependent on user and system settings. However, even in cases where log compression or color palette rendering is known, issues such as signal saturation, truncation and quantization may hamper acquisition of echo-power signals with appropriate proportionality to concentration (

with

Rognin et al., 2008

; Gauthier et al., 2011

; Payen et al., 2013

). Rognin et al., 2008

summarized that perfusion quantification based on linearized log-compressed video data can be equivalent to that based on raw data only under two conditions: (i) the gain settings should be adjusted in such a way to prevent data truncation or saturation in images, and (ii) a sufficiently high dynamic range of log compression (approx >45 dB) should be used. Similar conclusions were drawn by Gauthier et al., 2011

. In usual practice, both conditions can be fulfilled, and they also correspond to settings that provide good visualization of contrast perfusion. Proper data linearization permits the direct application of the indicator dilution theory based on the measured TIC, *I*(*t*), which is then related to the UCA concentration as$I\left(t\right)=aC\left(t\right)+b,$

(7)

with

*a*being a proportionality constant related to the total scattering cross-section of the microbubble population (Bouakaz et al., 1998

; Frinking and de Jong, 1998

), and *b*a baseline constant accounting for background intensity*.*This condition enables quantitative analysis of the hemodynamics based on the spatial distribution and temporal evolution of the relative bubble concentration. However, quantification based on the absolute relation between injected microbubble dose and measured concentration, as described in eqn (1), is not feasible because of the unknown values of*a*and*b*in eqn (7). In particular, it is difficult to account for the effect of attenuation, which is especially evident at high microbubble concentrations. Moreover, in many cases ultrasound systems apply non-linear and sometimes non-reversible signal processing or image processing, such as variations of log compression with depth, frequency or angular compounding and space-dependent equalization after scan conversion. In such cases, the linearization endeavor becomes futile, and results are bound to be inaccurate.## Temporal analysis

After a UCA bolus injection, the TIC describes the wash-in and wash-out of the microbubbles in the investigated area (Fig. 1). As a first approach, empirical semiquantitative parameters can be extracted from time and amplitude characteristics of the measured TIC. Although indirectly, these semiquantitative parameters are related to physiological quantities such as the injected contrast dose, blood flow and distribution volume. The most commonly used parameters are summarized in Table 1. In practice, however, the measured TIC is often corrupted by several sources of noise, such as acoustic reverberation, poor contrast mixing, shadowing and speckle noise (

Tang et al., 2011

); moreover, recirculation of the contrast agent often hides the wash-out. This makes direct calculation of semiquantitative parameters from the measured TIC unreliable. Fitting the measured TIC to a suitable mathematical model permits filtering out the noise and extracting the first pass of the contrast bolus. Moreover, when physiological models are used, the model parameters can be directly related to physiological quantities of interest, such as the MTT.Table 1Semiquantitative parameters typically extracted from a bolus injection investigation

Symbol | Description | Related to |
---|---|---|

AUC | Area under the curve | Blood flow, contrast dose |

TA | Appearance time | Blood flow |

PE | Peak enhancement | Blood volume, contrast dose |

TTP | Time-to-peak | Blood flow/volume |

WIT | Wash-in time | Blood flow/volume |

WOT | Wash-out time | Blood flow/volume |

WIR | Wash-in rate | Blood flow/volume |

WOR | Wash-out rate | Blood flow/volume |

Over the years, several models have been developed to describe the dilution of intravascular contrast agents, such as conventional UCAs (

Strouthos et al., 2010

; Harabis et al., 2013

; Turco et al., 2016

; Mischi et al., 2018

). In the following overview, we categorize them into compartmental and non-compartmental models.The first example of a compartmental model is the mono-compartment model in eqn (2). To describe more realistically the combination of atrium and ventricle in the heart, the mono-compartmental model can be extended to a two-compartment model by adding a compartment at the output (

With initial condition

with τ being the time constant of the compartments, and MTT =

where Γ(

von Reth and Bogaard, 1983

). From this, the *n*-compartmental model is easily obtained by a cascade of compartments, for which the input of the*i*th compartment is equal to the output of the (*i*– 1)th compartment (Fig. 5a) (von Reth and Bogaard, 1983

; Krejcie et al., 1996

). Assuming *n*chambers of equal volume*V*, the partial differential equation for the*i*th compartment reads$V\frac{d{C}_{i}\left(t\right)}{dt}=-F\left({C}_{i}\left(t\right)-{C}_{i-1}\left(t\right)\right).$

(8)

With initial condition

*C*(*t*= 0) =*C*_{0}=*m*/V in the first compartment and*C*(*t*= 0) = 0 elsewhere, the solution of (8) is given by an Erlang PDF multiplied by the AUC as$C\left(t\right)=\text{AUC}\frac{{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{t}{\tau}}{\left(\frac{t}{\tau}\right)}^{n-1}}{\tau (n-1)!},$

(9)

with τ being the time constant of the compartments, and MTT =

*nτ*, with*n*a positive integer (*n*∈ ℤ^{+}). Besides the compartmental approach, the same solution was obtained byMischi et al., 2008

by using a binomial distribution to describe the position of an indicator particle in a tube composed of *n*connected sections. Admitting a non-integer number of compartments (*n*∈ ℝ^{+}), eqn (9) is generalized by the Gamma distribution as (Mischi et al., 2008

)$C\left(t\right)=\text{AUC}\frac{{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{t}{\tau}}{\left(\frac{t}{\tau}\right)}^{n-1}}{\tau \phantom{\rule{0.16em}{0ex}}\Gamma \left(n\right)},$

(10)

where Γ(

*n*) is the Gamma operator, and MTT =*nτ*.Among the non-compartmental models, the log-normal distribution is one of the models most commonly used to fit TICs measured by DCE-US (

Wise, 1966

; Band et al., 1997

; Krix et al., 2004

; - Krix M.
- Plathow C.
- Kiessling F.
- Plathow C.
- Kiessling F.
- Herth F.
- Karcher A.
- Essiq M.
- Schmitteckert H.
- Kauczor H.J.
- Delorme S.

Quantification of perfusion of liver tissue and metastases using a multivessel model for replenishment kinetics of ultrasound contrast agents.

*Ultrasound Med Biol.*2004; 30: 1355-1363

Strouthos et al., 2010

; Harabis et al., 2013

; Turco et al., 2016

; Mischi et al., 2018

). In statistics, it is obtained as the product of a large number of independent random variables. Taking the time *t*as the independent variable, the mean μ and standard deviation σ of its logarithm are the controlling parameters of the log-normal distribution as (Chen et al., 2017

)$C\left(t\right)=\frac{\text{AUC}}{\sqrt{2\pi}\sigma t}{e}^{{\frac{(ln(t)-\mu )}{2{\sigma}^{2}}}^{2}}.$

(11)

To describe the transport of an indicator bolus from a large artery to the microvascular bed,

where μ and σ are the mean and standard deviation of the Gaussian input, and τ

Bassingthwaighte et al., 1966

proposed the use of the Lagged-normal function. Mathematically, the family of Lagged-normal functions is the result of the convolution of a Gaussian density function with *n*exponentials;*n*also defines the order of the Lagged-normal function. From a compartmental modeling viewpoint, this can be also seen as a cascade of*n*compartments whose input is a Gaussian PDF instead of a Dirac impulse (Fig. 5). This results in$C\left(t\right)=\frac{\text{AUC}}{\sqrt{2\pi {\sigma}^{2}}}{e}^{\frac{{}^{-\phantom{\rule{0.25em}{0ex}}{(t-\mu )}^{2}}}{2{\sigma}^{2}}}*\frac{1}{{\tau}_{1}}{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{1}{{\tau}_{1}}t}*\frac{1}{{\tau}_{2}}{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{1}{{\tau}_{2}}t}*...*\frac{1}{{\tau}_{n}}{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{1}{{\tau}_{n}}t}*,$

(12)

where μ and σ are the mean and standard deviation of the Gaussian input, and τ

_{1}, … ,τ*n*represent the time constants of the exponential functions. As described byDavis and Kutner, 1976

, analytical solutions of eqn (12) can be obtained for any *n*by using an error function approximation of the convolution integrals, from which the MTT can be calculated as$\text{MTT}=\mu +{\tau}_{1}+{\tau}_{2}+...+\phantom{\rule{0.16em}{0ex}}{\tau}_{n}.$

(13)

Another class of models that have been extensively used to describe UCA dilution are derived as a solution of the diffusion with drift equation (

Mischi et al., 2003

, Mischi et al., 2018

; Sheppard and Householder, 1951

; Turco et al., 2016

), which is given as$\frac{\partial C(x,t)}{\partial t}=\frac{\partial}{\partial x}\left(D\frac{\partial C(x,t)}{\partial x}\right)-v\frac{\partial C(x,t)}{\partial x}.$

(14)

Equation (14) describes the transport along the longitudinal direction

where μ is the MTT, and λ is the skewness parameter. The latter is related to the Péclet number Pe =

where the dispersion parameter

*x*of an indicator bolus in a tube of infinite length, where the microbubbles are dragged by a carrier fluid at a velocity,*v*, and diffuse with a diffusion coefficient,*D*, given by vascular flow profile and molecular diffusion projected over the longitudinal dimension. The local density random walk (LDRW) model and the first-time passage (FTP) model are both solutions of eqn (14), which mainly differ in the boundary condition at the output: while the LDRW allows for multiple passages of a microbubble across the detection site, the FTP assumes an absorbing barrier at the output so that only one passage of each microbubble is allowed (Sheppard, 1962

; von Reth and Bogaard, 1983

; Bogaard et al., 1986

; Mischi et al., 2004

, Mischi et al., 2007

). For constant *D*, the resulting LDRW and FTP models are given, respectively, as$C\left(t\right)=\text{AUC}\frac{{e}^{\lambda}}{\mu}\sqrt{\frac{\lambda \mu}{2\pi t}}{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{\lambda}{2}\left(\frac{t}{\mu}+\frac{\mu}{t}\right)},$

(15)

$C\left(t\right)=\text{AUC}\phantom{\rule{0.28em}{0ex}}{e}^{\lambda}\sqrt{\frac{\lambda \mu}{2\pi {t}^{3}}}{e}^{-\phantom{\rule{0.25em}{0ex}}\frac{\lambda}{2}\left(\frac{t}{\mu}+\frac{\mu}{t}\right)}.$

(16)

where μ is the MTT, and λ is the skewness parameter. The latter is related to the Péclet number Pe =

*vL*/*D =*2λ, which reflects the ratio of the contributions by convection and diffusion to the microbubble transport. When modeling UCA transport through a microvascular network, the dominant contribution to diffusion is given by the multiple transit times across the multipath intravascular trajectories, referred to as apparent diffusion or dispersion (Taylor, 1953

). Under the hypothesis that dispersion features of the contrast transport better reflect the microvascular architecture than perfusion, Kuenen et al., 2011

proposed the modified LDRW model (mLDRW), by which eqn (14) is solved locally, using as input boundary condition a Gaussian distribution in space of the contrast concentration, prior to entering the detection site. This results in$C\left(t\right)=\text{AUC}\sqrt{\frac{\kappa}{2\pi t}}{e}^{-\frac{\kappa}{2t}{(t-\mu )}^{2}}.$

(17)

where the dispersion parameter

*κ*=*λ*/*μ*=*v*_{l}^{2}/2*D*_{l}represents the ratio between the square of the local velocity,*v*_{l}, and the local dispersion coefficient,*D*_{l}. In a pilot study on prostate cancer patients, quantification of dispersion yielded higher diagnostic accuracy compared with perfusion-related parameters (Smeenge et al., 2011

).Although the analysis is most commonly performed on an ROI, providing an averaged assessment of the investigated area, parametric maps can also be obtained when the analysis is performed on a spatial grid, down to the single-pixel resolution. Besides permitting more accurate localization, parametric maps are particularly interesting in the assessment of tumor angiogenesis and microenvironment, for which heterogeneity has been reported to be a key feature (

Fidler, 2001

). However, this comes at the cost of computational time, as time-consuming iterative non-linear curve-fitting needs to be repeated for all pixels. Moreover, fitting single-pixel TICs is especially challenging because of lower signal-to-noise ratios. Next to standard least-squares approaches, several authors have introduced maximum-likelihood approaches that account for the statistics of arrival times (measured TIC) as well as noise, which has non-stationary characteristics depending on the UCA concentration levels (Barrois et al., 2013

; Kuenen et al., 2013a

; Bar-Zion et al., 2015

). In Figure 6 is an example of a parametric map of local dispersion obtained in prostate tissue.As an alternative to the bolus injection protocol, the contrast agent can be administered by a constant infusion, and the destruction–replenishment technique can then be used to assess perfusion. By this technique, a “negative” bolus is created by a series of high-pressure ultrasound pulses, referred to as a flash, which destroy all the microbubbles in the insonated area. The consequent inflow of microbubbles is sampled repeatedly at low ultrasound pressure so as to obtain a “replenishment” curve (Figure 3), which can be fitted by suitable models to extract perfusion parameters. Table 2 provides an overview of semi-quantitative parameters typically extracted from a destruction-replenishment investigation.

Table 2Semiquantitative parameters typically extracted from a destruction–replenishment investigation

Symbol | Description | Related to |
---|---|---|

A | Replenishment plateau | Blood fractional volume |

τ | Replenishment time constant | Blood flow |

β | Replenishment rate | Blood flow |

The first and simplest model of destruction–replenishment kinetics was proposed by

where

Wei et al., 1998

for assessment of myocardial perfusion. It consists of the step response of a mono-compartment dilution system (Fig. 2) given as$C\left(t\right)=A\left(1-{e}^{-\frac{t}{\tau}}\right),$

(18)

where

*A = aC0·*(see eqn [7]) represents the replenishment plateau intensity and is proportional to the fractional blood volume (blood volume per unit tissue volume), while τ is the replenishment time constant and is related to the blood flow,*F*, by τ =*V*/*F*, with*V*the vascular volume. Indeed, the model in eqn (18) corresponds to the response of a mono-compartmental system to an input of the form*Au*_{1}(*t*), with*u*_{1}(*t*) the unitary step function. Perfusion is derived as the ratio*A*/τ, which is proportional to blood flow per unit tissue volume. In fact, based on the combination of eqns (5) and (6), τ corresponds to the MTT of the insonated volume (Metoki et al., 2006

). In clinical practice, the replenishment rate β = 1/τ is more commonly adopted, and perfusion is derived as the product *A*β. Quantification of replenishment kinetics by the Wei model showed promise for differentiation between focal nodular hyperplasia and other hypervascularized liver focal lesions (Huang-Wei et al., 2006

) and between liver metastases and normal liver parenchyma (Krix et al., 2004

).- Krix M.
- Plathow C.
- Kiessling F.
- Plathow C.
- Kiessling F.
- Herth F.
- Karcher A.
- Essiq M.
- Schmitteckert H.
- Kauczor H.J.
- Delorme S.

Quantification of perfusion of liver tissue and metastases using a multivessel model for replenishment kinetics of ultrasound contrast agents.

*Ultrasound Med Biol.*2004; 30: 1355-1363

During imaging of microbubble replenishment, bubble destruction can occur and affect the measurement, depending on the adopted pressure and frame rate (

Shi et al., 2001

). Moreover, the establishment of a stable plateau is often complicated. To overcome these problems, Krix et al., 2003b

proposed the combination of a bolus injection with intermittent imaging, by which repeated flashes are transmitted during a bolus passage, and only one measurement is taken for each replenishment at varying delay. After compensating for the baseline variation by TIC fitting, a replenishment curve can be reconstructed that is not affected by the adopted imaging sequence.Extension to a multicompartment model was proposed by

where

Lucidarme et al., 2003

to take into account the refilling evolution in different subvolumes of the imaging and destruction areas. Considering a cascade of *n*subvolumes, describing the full path covered by the microbubbles prior to entering the observed area, the concentration in the*n*th compartment is given by$\frac{d{C}_{n}\left(t\right)}{dt}=\frac{F}{V}{C}_{n-1}\left(t\right)-\frac{F}{V}{C}_{n}\left(t\right)-\lambda {C}_{n}\left(t\right),$

(19)

where

*C*_{n – 1}(*t*) is the microbubble concentration in the subvolume*n*– 1,*F*is the blood flow,*V*is the blood volume and λ is the rate of microbubble destruction. This model can be solved analytically by assuming*Cn*(*t*= 0) for all*n*≠ 0 and is equal to*C*_{0}for*n*= 0, where*C*_{0}represents the concentration of microbubbles entering the first subvolume within the ultrasound field (Lucidarme et al., 2003

).A model describing a more complex vascular architecture, accounting for multiple vessels at different angles and with different refilling dynamics, was proposed by

with

Krix et al., 2003a

. The model separates between a linear part, which ends when the vessels perpendicular to the imaging plane are completely refilled (*t*_{max}), and a non-linear part, for the replenishment of the remaining vessels, which starts at time*t*_{min}. The resulting piecewise model is described by$\{\begin{array}{l}C\left(t\right)=M\sum _{i}{g}_{i}{v}_{i}t=M{v}_{\text{mean}}t\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}t<{t}_{\text{max}}\\ C\left(t\right)=max\left[1-\sum _{i}{g}_{i}\frac{{d}^{2}}{3}{v}^{2}{t}^{2}\right]\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.28em}{0ex}}t>{t}_{\text{min}}\end{array},$

(20)

with

*gi*being the weight of each blood velocity*vi, M*a constant and*d*the effective width of the ultrasound beam. Promising applications of this model include perfusion assessment in liver metastases and skeletal muscles (Krix et al., 2004

, - Krix M.
- Plathow C.
- Kiessling F.
- Plathow C.
- Kiessling F.
- Herth F.
- Karcher A.
- Essiq M.
- Schmitteckert H.
- Kauczor H.J.
- Delorme S.

Quantification of perfusion of liver tissue and metastases using a multivessel model for replenishment kinetics of ultrasound contrast agents.

*Ultrasound Med Biol.*2004; 30: 1355-1363

Krix et al., 2005

).Besides the complexity of the vascular tree, the shape of the refill curve is also influenced by the point-spread function of the imaging system and by the effective shape of the imaging and bubble destruction zones (

where

Potdevin et al., 2004

). To take the characteristics of the ultrasound beam into account, more complex models were proposed by Arditi et al., 2006

and Hudson et al., 2009

. The former proposed a replenishment function averaged over the probability distribution *P*(τ) of indicator transit times in the vascular tree, which is described by$C\left(t\right)=A\underset{0}{\overset{\infty}{\int}}P\left(\tau \right)\text{perf}\left[1.94\frac{Kd}{2\tau}\left(t-\tau \right)\right]d\tau ,$

(21)

where

*A*is an estimate of the blood volume;*K*is a transmit–receive diffraction parameter, which depends on beam elevation and ultrasound wavelength;*d*is the width of the destruction slice, τ is the time needed to refill half of the bubble destruction slice, and perf(*q*)*=*0.5⋅(1+erf(*q*))*,*with erf(⋅) the error function. As suggested by the authors, any model describing the distribution of indicator transit times (*e.g.,*bolus kinetic models) can be used for*P*(τ).Hudson et al., 2009

proposed a more general representation of the replenishment kinetics by averaging over the investigated volume a flow function *F*(

*x, y, z, t*) weighted by a beam function

*B*(

*x, y, z*) as

$C\left(t\right)=\underset{V}{\int}B\left(x,y,z\right)F\left(x,y,z,t\right)dV.$

(22)

A possible description for the flow function is derived by a bifurcating model of the vascular tree (

Qian and Bassingthwaighte, 2000

), which results in a log-normally distributed velocity profile (Hudson et al., 2009

). The beam function is typically modeled as a sinc function, parametrized by the beam profile in the elevation plane and the shape of the boundary between the destruction and imaging zones, which is band-pass filtered to truncate the main lobe (Hudson et al., 2009

).A framework combining the bolus and replenishment methods was proposed by

Here

Jirik et al., 2013

. Although performing both the bolus injection and constant infusion protocols is not standard in clinical practice, this framework can potentially provide absolute values of blood flow and volume from linearized ultrasound data, without need for contrast calibration. Defining the arterial input function, AIF(*t*), as the contrast concentration in the local arterial input to the ROI, the contrast concentrations in the ROI following a bolus injection and a constant infusion are modeled as, respectively,${C}_{\text{bolus}}\left(t\right)=F\text{AIF}\left(t\right)*R\left(t\right),$

(23)

${C}_{\text{repl}}\left(t\right)=F{C}_{0}{u}_{1}\left(t\right)*R\left(t\right).$

(24)

Here

*F*is the blood flow per unit tissue volume (mL*min/mL),*C*_{0}is the plateau concentration before the destruction pulse and*R*(*t*) is the tissue residue function, which relates to the impulse response of the dilution system as$R\left(t\right)=1-\underset{0}{\overset{t}{\int}}h\left(\xi \right)d\xi .$

(25)

Simultaneous fitting of eqns (23) and (24) permits estimation of the MTT and absolute blood volume

*V*. However, a key assumption is that the same proportionality constant between the measured linearized acoustic intensity and the contrast agent concentration can be assumed for AIF(*t*) and*C*(*t*). Often this is not realistic because the AIF shows higher concentrations that are close to the upper limit of the dynamic range.## Spatial analysis

Temporal features of a measured TIC relate to perfusion and UCA dispersion, which reflect the underlying vascular architecture and distribution volume. In abnormal vascular tissue, for example, due to angiogenesis, local features of the microvascular network such as density, tortuosity, vessel size and multipath trajectories influence local perfusion and UCA dispersion kinetics. Therefore, direct analysis of the spatial distribution of UCAs may reveal additional information compared with that revealed when only the temporal aspects are considered. For instance, vascular heterogeneity is a key feature of many solid tumors (

Fidler, 2001

), which may also influence the response to therapy (Knieling et al., 2013

).Although not as established as temporal analysis, a few studies have shown the potential of spatial analysis of DCE-US, with applications in cancer diagnostics and atherosclerotic plaque characterization (

Molinari et al., 2010

; Zhang et al., 2014

; Saidov et al., 2016

). This is achieved by either direct analysis of the spatial distribution of UCAs at peak enhancement or by analysis of parametric maps, obtained by maximum intensity projection or other quantitative temporal analysis.Three-dimensional characterization of thyroid nodules by DCE-US was proposed

where IP is the planar angle between two vectors connecting three nodes, and TP is the angle defining the torsion of the curve in a node. From eqn (26), the inflection count metric, ICM, is calculated by dividing DM for the number of flexes in the path. Validation of this method on 20 thyroid nodules revealed that spatial features extracted by DCE-US aid the distinction between benign and malignant nodules (

Molinari et al., 2010

as a combination of morphological operations to reconstruct a skeleton of the vasculature and computation of quantitative features from the obtained skeleton (Fig. 7). By their approach, the skeleton is mapped as a series of vascular trees *Ti*, each containing a sequence*m*of nodes*pj*, that is*Ti*= [*p*_{1},*p*_{2}, …,*p*_{m}]. From this, features such as number of trees, density of the vascular structure, number of branching nodes and average vessel radius are extracted to analyze the vascular network. Tortuosity is characterized in 2-D, by calculation of the distance metrics, DM, and in 3-D by the sum of angle metrics, SOAM, as$\text{DM}=\frac{|{p}_{m}-{p}_{i}|}{{\sum}_{i=1}^{m-1}|{p}_{i}-{p}_{i+1}|},$

(26)

$\text{SOAM}=\frac{{\sum}_{k=1}^{m-3}\sqrt{\mathrm{I}{{\mathrm{P}}_{k}}^{2}+\mathrm{T}{{\mathrm{P}}_{k}}^{2}}}{{\sum}_{k=1}^{m-1}|{p}_{k}-{p}_{k+1}|},$

(27)

where IP is the planar angle between two vectors connecting three nodes, and TP is the angle defining the torsion of the curve in a node. From eqn (26), the inflection count metric, ICM, is calculated by dividing DM for the number of flexes in the path. Validation of this method on 20 thyroid nodules revealed that spatial features extracted by DCE-US aid the distinction between benign and malignant nodules (

Molinari et al., 2010

).Similarly,

where

Zhang et al., 2014

proposed quantitative analysis of the UCA spatial distribution for characterization of atherosclerotic plaques. By the proposed method, TICs are first pre-filtered to separate the perfusion TIC from the cardiac cycle curve; the intraplaque neovascularization (IPN) is then segmented by a framework comprising temporal averaging, interactive delineation and binarization by thresholding. From the obtained binary image, which contains *n*_{p}“white” pixels representing the IPN, quantitative features can be derived such as the center deviation degree (CDD) the radial deviation degree (RDD) and the dispersion degree, DD, as$\text{CDD}=\frac{1}{{n}_{p}}\sum _{i=1}^{{n}_{p}}\left(\frac{{d}_{i}}{{D}_{i}}\right),\phantom{\rule{0.33em}{0ex}}\phantom{\rule{-0.5em}{0ex}}\text{RDD}=\sqrt{\frac{1}{{n}_{p}-1}\sum _{i=1}^{{n}_{p}}{\left(\frac{{d}_{i}}{{D}_{i}}-\text{CDD}\right)}^{2}},\phantom{\rule{0.28em}{0ex}}\phantom{\rule{-0.5em}{0ex}}\text{DD}=\frac{1}{R{n}_{p}}\sum _{i=1}^{{n}_{p}}\left({r}_{i}\right),$

(28)

where

*d*_{i}is the distance between the center of the plaque and the*i*th white pixel,*Di*is the distance from the center of the plaque and the point of intersection between the plaque boundary and the extension line from the*i*th pixel; thus, CDD represents the average normalized distance from the IPN to the plaque, and RDD its standard deviation;*ri*is the distance between the center of the IPN and*i*th pixel, and*R*is the mean radius of the plaque. A combination of these parameters with the ratio between IPN and plaque areas (*A*_{r}), extracted from standard B-mode imaging, was used to obtain the combined area ratio as*A*_{r,com}=*A*_{r}·DD/CDD, which achieved a sensitivity and specificity of 0.79 and 1.00, respectively, in the distinction between low-grade and high-grade plaques.Based on the well-known concept that the vascular tree can be described by a fractal bifurcating tree, fractal analysis has been proposed for spatial characterization of DCE-US. In fact, previous studies on microscopic analysis of tissue samples have indicated that the fractal dimension, FD, is higher in chaotic vascular networks, such as those typical of malignant tissue, than in more regular networks, typical of benign tissue (). On this basis,

where

Saidov et al., 2016

proposed a method to extract the FD from perfusion maps obtained from DCE-US. To calculate the FD, a sample of size *L*is first divided into subsamples of smaller sizes (*e.g., L*/4,*L*/16,*L*/64); then, the relative dispersion, RD, defined as the ratio between the standard deviation and the mean value of the extracted parameter, can be modeled at different scales*m*by a power law of the FD as$\frac{\text{RD}\left(m\right)}{\text{RD}\left({m}_{\text{ref}}\right)}={\left(\frac{m}{{m}_{\text{ref}}}\right)}^{1-\text{FD}},$

(29)

where

*m*_{ref}is an arbitrary sample size chosen as reference. From eqn (29), FD can be extracted by curve fitting. In a feasibility study in mice, good agreement was shown between the immunohistological analysis of tumor tissue, reflecting geometric properties of the microvascular network, and the FD estimated on parametric maps of wash-in rate and peak enhancement (Saidov et al., 2016

).## Spatiotemporal analysis

The transport of UCA particles across the vascular bed is inherently a 3-D spatiotemporal process. Thus, exploitation of both spatial and temporal characteristics of DCE-US loops, possibly in 3-D, may provide a more complete description of the UCA transport process and a better characterization of perfusion, contrast dispersion and vascular architecture.

An early example of spatiotemporal analysis estimates the fractal dimension (FD) over time for better characterization of the vascular anatomy and perfusion in ischemic tissue (

where

Charalampidis et al., 2006

). Differently from eqn (29), FD is calculated by a modified variation method, repeated over all frames so as to obtain a temporal FD sequence, FD(*t*), of length*L*_{FD}. This permits the extraction of spatiotemporal features by averaging over different time windows as${F}_{\text{FD}}^{\left(i\right)}=\frac{1}{T}\sum _{t={t}_{i}-T+1}^{{t}_{i}}\text{FD}\left(t\right),$

(30)

where

*T*is chosen according to the oscillation period empirically observed in FD(*t*), and the time points*t*_{i}are taken at different percentages of the full time window of FD(*t*)_{.}The idea of analyzing TICs in a local spatiotemporal area within a DCE-US loop was first proposed by

Angelelli et al., 2011

and used for the development of an interactive software to aid analysis and characterization of liver lesions. Defining a spatiotemporal pixel neighborhood, *N*(**x**,*t*), at location**x**and time*t*, they extracted several enhancement metrics in the neighborhood, including statistical measures such as the average, first and third quartiles of enhancement and percentage of enhancement, calculated as percentage of pixels in*N*(**x**,t) with enhancement higher than an empirically chosen threshold. From these spatially averaged enhancement metrics, more robust perfusion parameters can be extracted.While in

Angelelli et al., 2011

similarity analysis was used as an intermediate step to aid manual lesion segmentation, Mischi et al. (Kuenen et al., 2013b

; Mischi et al., 2012

) proposed a fully automated method to extract similarity maps from DCE-US loops and found that spatial similarity is related to dispersion of UCAs through the multipath trajectories of the microvascular network (Mischi et al., 2012

; Kuenen et al., 2013a

, Kuenen et al., 2013b

). Besides linear correlation (Kuenen et al., 2013b

), they proposed to measure similarity between neighboring pixels by their spectral coherence, defined as the correlation coefficient of the TIC amplitude spectra in a given frequency bandwidth. Focusing on the amplitude spectrum overcomes any influence of the UCA appearance time, which is instead captured by the phase spectrum.As illustrated in Figure 8, a ring-shaped kernel is used to obtain parametric maps of temporal correlation and spectral coherence (

Mischi et al., 2012

; Kuenen et al., 2013b

), by comparing the TIC in each pixel with neighboring TICs within a ring of 1-mm inner radius and 2.5-mm outer radius. While the kernel shape is chosen to overcome any dependence on the direction of blood perfusion, the kernel size is based on the well-known concept of “angiogenic switch”: tumors require angiogenesis to grow beyond the size limit of 2–3 mm in diameter, and thus, a minimum spatial resolution of 3 mm is necessary to visualize early angiogenesis. Prior to the analysis, the images are regularized by Wiener deconvolution to obtain a homogeneous and isotropic resolution (speckle size) across the image (Kuenen et al., 2013b

). This is essential for correct application of similarity measures. Validation of the method for prostate cancer localization revealed that the proposed spatiotemporal features yielded better performance compared with time-amplitude features, related to blood perfusion and UCA dispersion (Mischi et al., 2012

; Kuenen et al., 2013a

, Kuenen et al., 2013b

). Moreover, an analytical link was also established between spatiotemporal dispersion measures by similarity analysis and temporal measures of the dispersion coefficient by eqn (17) (Kuenen et al., 2013b

).To explore possible non-linear relationships between neighboring pixels, a more general statistical measure of similarity, the mutual information, was also proposed (

Schalk et al., 2017

). This can be obtained by modeling the center TIC as a random variable *A*obtained by sampling*M*independent identically distributed data points, with*M*the number of frames, and the surrounding TICs in the kernel as sampled from a random variable*B*^{N}, with*N*the number of surrounding pixels in the kernel (Schalk et al., 2017

). Similarity analysis can also be extended to 3-D by substituting the ring kernel by a spherical kernel (Schalk et al., 2015

, Schalk et al., 2018

).The proposed estimators of UCA dispersion-based spatiotemporal similarity, as well as those based on temporal TIC analysis by eqn (17), cannot distinguish between the contributions of convection and dispersion to the transport process, as the estimated dispersion parameter, κ refers to the ratio between convection (

*v*^{2}) and dispersion (*D*). More recently,van Sloun et al., 2017a

proposed a local system identification approach using the spatiotemporal impulse response (Green's function) of the 1-D convective dispersion equation in eqn (14), which is given as$g(x,t|v,D)=\frac{1}{\sqrt{4\pi Dt}}exp\left(-\frac{{(x-vt)}^{2}}{4Dt}\right).$

(31)

Adopting again a ring-shaped kernel with the same rationale as described above, the center TIC is considered the input of the system, while the TICs in the kernel are all possible outputs. Based on the known distance between the ring center and each pixel in the kernel, and including causal contributions only, input–output analysis by model-based system identification enables separate estimation of

*v*and*D*, from which the Péclet number can also be estimated (van Sloun et al., 2017a

).In addition to magnitude, directionality of blood velocity can also be assessed by estimation of the local propagation of the UCA bolus (

where is a vector describing time delays between all the TICs in the circular kernel,

van Sloun et al., 2017b

). Assuming the convective component to be dominant in the UCA transport, the average velocity vector $\overrightarrow{v}=\left[{v}_{x,}{v}_{y}\right]$ in a circular kernel of radius *R*is described by the relationship${\overrightarrow{v}}^{T}\overrightarrow{\tau}=\mathbf{L}$

(32)

where is a vector describing time delays between all the TICs in the circular kernel,

**L**is a matrix representing the inter-pixel distance, and*T*represent the transpose operation. At each pixel, the time delay is first estimated by maximization of the cross-correlation between each pair of TICs in the kernel, and then least-square minimization of eqn (34) is applied for estimation of $\overrightarrow{v}$.Despite the limitations imposed by a 2-D approach, where only projections along the imaging plane can be visualized, statistical analysis of the obtained velocity vector field allows characterization of the heterogeneous and chaotic nature of angiogenic vasculature. In particular, the Shannon entropy

*H*(*V*), reflecting the degree of disorder in the velocity field pattern, and the conditional entropy*H*(*V*_{1}|*V*_{2}), reflecting the predicting value of a velocity vector (represented by random variable*V*_{1}) with respect to its surroundings (V_{2}), have been adopted to reflect vascular heterogeneity (van Sloun et al., 2017b

).The use of eqn (34) can be extended to the analysis of 3-D DCE-US data by replacing the circular kernel by a spherical one, leading to the estimation of all three components of the local velocity vector, that is, $\overrightarrow{v}=\left[{v}_{x},{v}_{y},{v}_{z}\right]$ (

describing the trajectory of a particle starting at

van Sloun et al., 2018

). From the velocity vector field, the streamlines describing the flow trajectories can be extracted to infer the underlying vascular architecture (van Sloun et al., 2018

). This can be achieved by a method referred to as DCE-US tractography, where the streamlines of the velocity field are obtained by solving the ordinary differential equation${\partial}_{t}x\left(t\right)=v\left[x\left(t\right)\right]\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}\text{with}\phantom{\rule{0.33em}{0ex}}x\left(t=0\right)={x}_{0}$

(33)

describing the trajectory of a particle starting at

*x*_{0}and moving within the vector field given by*v*[*x*(*t*)]. From the obtained tractographs, features such as the ICM and spatial density of trajectories (*i.e.,*number of trajectories per unit volume) can be extracted to assess vessel tortuosity and density. In a proof-of-concept validation in prostate cancer patients, 3-D parametric maps of the proposed tractographic features showed good agreement with the presence of malignancy as indicated by histopathology (van Sloun et al., 2018

). By extension of the analysis to 3-D convective-dispersion velocity fields, Wildeboer et al., 2018b

showed probabilistic-tractography representations of velocity and dispersion in the prostate. This required solving the convective-dispersion equation in 3-D, based on a 3-D DCE-US acquisition. Figure 9 shows an example of 3-D UCA dispersion and tractographic analysis of a prostate.The first spatiotemporal model describing UCA transport kinetics as a 3-D convective-dispersion process was proposed and applied to prostate cancer localization (

Wildeboer et al., 2018a

). This model is still based on the interpretation of UCA transport as a convective dispersion process, similar to eqn (14), but in 3-D. By assuming *D*and*v*to be locally constant and assuming*D*_{xy}=*D*_{yx},*D*_{xz}=*D*_{zx}and*D*_{yz}=*D*_{zy}, eqn (14) expands as$\begin{array}{c}{\partial}_{t}C\left(t\right)={D}_{xx}{\partial}_{xx}C\left(t\right)+2{D}_{xy}{\partial}_{xy}C\left(t\right)+2{D}_{xz}{\partial}_{xz}C\left(t\right)+{D}_{yy}{\partial}_{yy}C\left(t\right)\hfill \\ \phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}+2{D}_{yz}{\partial}_{yz}C\left(t\right)+{D}_{zz}{\partial}_{zz}C\left(t\right)-{v}_{x}{\partial}_{x}C\left(t\right)-{v}_{y}{\partial}_{y}C\left(t\right)-{v}_{z}{\partial}_{z}C\left(t\right).\hfill \end{array}$

(34)

The concentration gradients are computed by convolving the acquired data with Gaussian derivatives in space and time. In this way, a linear model is obtained and solved by least-squares minimization. From the obtained estimates of the dispersion tensor

**D**and velocity vector $\overrightarrow{v}$, the apparent dispersion coefficient,*D*_{ADC}, and the velocity magnitude,*v*, are calculated as${D}_{ADC}=\left|\frac{1}{3}\left({D}_{xx}+{D}_{yy}+{D}_{zz}\right)\right|,\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}v=\sqrt{{v}_{x}^{2}+{v}_{y}^{2}+{v}_{z}^{2}}.$

(35)

In a preliminary validation in prostate cancer patients, improved diagnostic performance was demonstrated compared with 3-D similarity analysis (

Wildeboer et al., 2018a

). Rather than using Gaussian high-pass filters, an alternative approach has also been proposed that makes use of a finite-element formulation of the convective-dispersion equation, solved by use of spatiotemporal boundary conditions measured in 3-D (Wildeboer et al., 2018b

). The use of the Green's function in eqn (33) has also been proposed to estimate convection and dispersion in 3-D (Wildeboer et al., 2019

). However, to deal with low signal-to-noise ratios and spatiotemporal sampling, spatiotemporal filtering by singular value decomposition analysis is required prior to the analysis (Wildeboer et al., 2019

).An alternative method for quantification of UCA velocity describes the spatiotemporal changes in ultrasound gray-level intensity with the fluid-dynamic model (

by which the temporal variations in intensity,

where

de Senneville et al., 2018

)${\partial}_{t}I+\overrightarrow{v}\nabla I=0,$

(36)

by which the temporal variations in intensity,

*I*, are taken into account by the transient term, while the spatial variations are represented by the transport term $\overrightarrow{v}\nabla I$. Similar tovan Sloun et al., 2017a

, van Sloun et al., 2017b

), the UCA transport is here considered as purely convective; in fact, if a linear relationship between the contrast concentration and the gray-level intensity can be assumed, eqn (36) is equivalent to a 2-D convective-dispersion equation where the dispersion coefficient is set to zero. As the proposed transport equation is underdetermined, estimation of the vector velocity field is achieved through numerical resolution by the Euler–Lagrange method. To summarize the pixelwise information in an ROI, Γ, over a defined time window, Δ*T*, the overall microbubble transport amplitude is obtained from the velocity vector field as$\gamma =\frac{{\sum}_{t={t}_{0}}^{{t}_{0}+\Delta T}{\sum}_{\overrightarrow{r}\in \Gamma}I(\overrightarrow{r},t){\parallel V(\overrightarrow{r},t)\parallel}_{2}}{{\sum}_{t={t}_{0}}^{{t}_{0}+\Delta T}{\sum}_{\overrightarrow{r}\in \Gamma}I(\overrightarrow{r},t)},$

(37)

where

*t*_{0}is the bolus arrival time, $\overrightarrow{r}\phantom{\rule{0.33em}{0ex}}\phantom{\rule{0.33em}{0ex}}$represents the spatial location and the symbol ∥ · ∥_{2}represents the L2 norm. In a pre-clinical study on assessment of placenta perfusion, microbubble transport fields were successfully visualized, showing agreement with expected perfusion patterns in normal versus ligature placentas; moreover, the proposed parameter γ outperformed perfusion-related parameters obtained from standard temporal analysis for classification of placental insufficiency (de Senneville et al., 2018

).## Molecular analysis

The recent introduction of targeted ultrasound contrast agents (tUCAs), which can specifically bind to molecules present in the vessel wall or in the bloodstream, has permitted quantification of pathophysiological processes down to the molecular level (

Deshpande et al., 2010

; Pochon et al., 2010

; Abou-Elkacem et al., 2015

; Smeenge et al., 2017

; Willmann et al., 2017

). Novel tUCAs are obtained by decoration of the microbubble shell with targeting ligands, designed to bind specifically to receptors expressed on endothelial cells. To date, the only clinical-grade tUCA targets the vascular endothelial growth receptor factor 2 (VEGFR2), which is overexpressed on the endothelial wall of many solid tumors, including breast, ovarian, prostate, colorectal and pancreatic cancer (Pochon et al., 2010

; Abou-Elkacem et al., 2015

; Smeenge et al., 2017

; Willmann et al., 2017

).Analysis of TICs derived from molecular DCE-US (

*i.e.,*DCE-US with injection of tUCAs) is complicated by the fact that two populations of microbubbles with different kinetics need to be distinguished, namely, the free microbubbles in the circulation, whose transport can be assumed equivalent to that of conventional UCAs, and the bound microbubbles, which have attached to the vessel wall.Under the assumption that the concentration of bound microbubbles is proportional to the concentration of the target molecule, quantification in molecular DCE-US has been based mainly on the assessment of microbubble binding. As a first approach, the late enhancement, that is, the acoustic intensity in the late phase of the DCE-US loop, has been used for semiquantitative assessment of binding (

Abou-Elkacem et al., 2015

; Smeenge et al., 2017

; Willmann et al., 2017

). This requires long acquisition times of several minutes to ensure that all microbubbles in the circulation have washed out from the investigated area; moreover, the reproducibility is limited by the dependence on scanner and operator settings and on the choice of the analysis time point.To isolate the signal coming from bound microbubbles, an approach similar to the destruction–replenishment method has been proposed (

Frinking et al., 2012

; Abou-Elkacem et al., 2015

). By this approach, as illustrated in Figure 10, a high-pressure ultrasound burst is applied in the late phase to destroy all the microbubbles in the investigated area. Assuming that the subsequent replenishment is owing only to circulating microbubbles, the difference in the acoustic signal before and after the application of the destruction burst (*i.e.,*the differential targeted enhancement) is used to assess microbubble binding (Frinking et al., 2012

; Abou-Elkacem et al., 2015

). Table 3 provides an overview of semi-quantitative parameters typically extracted from molecular DCE-US.Table 3Semiquantitative parameters extracted from molecular dynamic contrast-enhanced ultrasound

Symbol | Description | Related to |
---|---|---|

LE | Late enhancement | Blood volume, contrast dose, contrast binding |

dTE | Differential targeted enhancement | Blood volume, contrast dose, contrast binding |

More quantitative approaches are based on mathematical modeling of microbubble transport and binding. By compartmental modeling, quantification has been achieved by dividing the intravascular space into two “non-physical” compartments: a free compartment accounting for the concentration of microbubbles in the circulation,

*C*_{f}(*t*), and a binding compartment, accounting for the concentration of bound microbubbles,*C*_{b}(*t*).The first example of a compartmental model for molecular DCE-US was proposed by Sirsi et al. (

where

Sirsi et al., 2012

) and applied by Chen et al., 2012

for assessment of microbubble binding in kidney vasculature. The proposed model is described by the equations${C}_{t}\left(t\right)={C}_{f}\left(t\right)+{C}_{b}\left(t\right),$

(38)

${C}_{f}\left(t\right)=\frac{{C}_{0}{k}_{1}}{\left({k}_{2}+{k}_{3}\right)-{k}_{1}}\left[{e}^{-{k}_{1}t}-{e}^{-\left({k}_{2}+{k}_{3}\right)t}\right],$

(39)

$\begin{array}{c}{C}_{b}\left(t\right)={C}_{0}{k}_{1}{k}_{3}\frac{{e}^{\left({k}_{4}-{k}_{1}\right)t}-1}{{{k}_{1}}^{2}-{k}_{1}{k}_{2}-{k}_{1}{k}_{3}-{k}_{1}{k}_{4}+{k}_{2}{k}_{4}+{k}_{3}{k}_{4}}\hfill \\ \phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.28em}{0ex}}\phantom{\rule{0.16em}{0ex}}+\phantom{\rule{0.33em}{0ex}}{C}_{0}{k}_{1}{k}_{3}\frac{{e}^{\left({k}_{4}-{k}_{3}-{k}_{2}\right)t}-1}{{{k}_{2}}^{2}+{{k}_{3}}^{2}-{k}_{1}{k}_{2}-{k}_{1}{k}_{3}+{k}_{1}{k}_{4}+2{k}_{2}{k}_{4}-{k}_{3}{k}_{4}},\hfill \end{array}$

(40)

where

*C*_{0}is the dose of contrast agent input to the free compartment,*k1*describes the contrast influx in the tissue vasculature and*k2*and*k4*describe the efflux of microbubbles from the free and bound compartment, respectively. The transport from the circulating and bound compartments is assumed unidirectional (*i.e.,*there is no unbinding) and described by*k*_{3}. Parameter estimation is performed by fitting simultaneously eqns (38)–(40) to TICs and time-fluctuation curves, which are derived by a dedicated speckle-decorrelation procedure (Chen et al., 2012

; Sirsi et al., 2012

). A relative measure of microbubble binding is given by the adhesion ratio AR = *k*_{2}*k*_{3}+*k*_{2}(Sirsi et al., 2012

).Another compartmental approach was proposed by

where

Turco et al., 2017

and applied for angiogenesis imaging in prostate cancer and therapy monitoring in colorectal cancer (Turco et al., 2019

). Focusing on the first pass of a tUCA bolus, the microbubble transport in the circulation compartment is modeled by the mLDRW model (see eqn [17]), and the transfer from the circulation to the bound compartment by the unidirectional binding rate constant *K*_{b}. The resulting “first-pass binding model” provides an analytical solution for*C*_{t}(*t*) as${C}_{t}\left(t\right)={v}_{f}{C}_{\text{mLDRW}}\left(t\right)+{K}_{b}{u}_{1}\left(t\right)*{C}_{\text{mLDRW}}\left(t\right),$

(41)

where

*v*_{f}is the fractional volume of free microbubbles and*C*_{mLDRW}(*t*) is the concentration of microbubbles in the circulation as modeled in eqn (17).Non-compartmental approaches typically include the empirical combination of intravascular contrast dilution models, such as those described under Temporal Analysis, with a mathematical function, which takes into account the contribution of bound microbubbles to the total acoustic intensity. The first example of such a model was developed to describe the retention of microbubbles in the myocardium, which can occur in the case of cardiac pathologies such as ischemia and coronary stenosis (

where

Lindner et al., 1998

). The total concentration of microbubbles is modeled as the sum of a first-order Gamma variate function and its integral, weighted by the retention fraction *f*as${C}_{t}\left(t\right)=f\underset{0}{\overset{\infty}{\int}}{C}_{\text{GV}}\left(t\right)dt+\left(1-f\right){C}_{\text{GV}}\left(t\right),$

(42)

where

*C*_{GV}(*t*) =*At*e^{–}*αt*, with*A*and α empirical parameters. This model has also been extended to describe a double-bolus experiment, in which a constant infusion of a conventional UCA is combined with the bolus injection of a tUCA. The combined model allows for separate assessment of blood flow and volume and of microbubble binding by the retention fraction,*f*(Carr et al., 2011

). This is especially relevant to study biologic processes such as angiogenesis, which influence blood perfusion. Figure 11 illustrates an application of this model for assessment of the angiogenic response to ischemic injury in wild-type and diabetic mice, obtained from an Institutional Animal Care and Use Committee-approved pre-clinical study. Similarly, the sum of a first-order Lagged-normal function (see eqn [12]) and its integral was proposed to model microbubble persistence in the myocardium for assessment of myocardial inflammation (Fisher et al., 2002

).In another empirical approach, the LDRW model in eqn (15) is used to describe the kinetics of microbubbles in the circulation, while microbubble binding is described by a ramp function, leading to (

where β and Δ

Sugimoto et al., 2012

)${C}_{t}\left(t\right)={C}_{\text{LDRW}}\left(t\right)+\beta \left(t-\Delta t\right)$

(43)

where β and Δ

*t*are empirical parameters describing the slope and time-delay of the ramp function, respectively.## Multiparametric analysis by machine learning

As multiple parameters obtained through quantification of DCE-US recordings have been found to correlate with specific pathologies, multiparametric image analysis has come into view. Multiparametric analysis entails the combination of multiple parameters to improve sensitivity and specificity for a certain diagnosis. A few groups have applied this in a cognitive manner, asking readers to weigh different DCE-US-derived parametric maps during the decision-making process (

Postema et al., 2015

). Machine learning, on the other hand, offers a data-driven approach to optimally combine parameters and exploit their complementarity for a given task.To perform classification, supervised machine-learning algorithms have to be trained, that is, the underlying weights of the algorithm have to be tuned based on a set of labeled training data. Moreover, the design of the algorithm itself is often dependent on a few hyperparameters that determine its structure. It is important to note that the performance of such algorithms can be assessed only in an independent test data set not previously seen by the algorithm. Examples of widely used machine-learning algorithms are

*k*-nearest neighbor (*k*-NN) classifiers (which classify new data according to the label of the*k*most comparable training instances [Altman, 1992

]), support vector machines (SVMs, which classify new data according to their position with respect to a multiparametric hyperplane separating benign and malignant data [Cortes and Vapnik, 1995

]), random forests (which classify new data based on consensus by an ensemble of decision trees trained on different partitions of the training data [Ho, 1995

]) and artificial neural networks [ANNs], which classify new data according to a series of fully connected non-linear functions [Hornik, 1991

]). Deep neural networks (*i.e.,*deep learning) form a family of ANNs that is probably the most powerful approach to machine learning. They consist of multiple layers and possibly tens of thousands of weights, but generally require a large number of training data (Altman, 1992

).Encouraging results have been reported on multiparametric combination of quantitative ultrasound features for disease classification in various organs. Several authors have studied the combination of DCE-US-derived (semi)quantitative features for liver lesions classification, using either SVMs (

Căleanu et al., 2014

; - Căleanu C.D.
- Simion G.
- David C.
- et al.

A study over the importance of arterial phase temporal parameters in focal liver lesions CEUS based diagnosis.

in: 2014 11th International Symposium on Electronics and Telecommunications (ISETC), Timisoara, Romania, November 14–15, 2015, Piscataway, NJ IEEE,
2014

Kondo et al., 2017

) or ANNs (Shiraishi et al., 2008

; Sugimoto et al., 2009

), usually reaching an appreciable accuracy in distinguishing hepatocellular carcinoma, hemangioma and benign lesions. In the prostate, Wildeboer et al., 2017

found that the combination of (in particular) perfusion- and dispersion-related DCE-US parameters by means of a Gaussian mixture model (GMM) improved the accuracy of prostate cancer classification. GMMs model the probability density distribution of benign and malignant instances over multiple parameters by a mixture of multidimensional Gaussians, and classify new pixels by assessing the likelihood of their being drawn from the benign or malignant distribution. Factor analysis and subsequent multiple linear regression were used by Zhang et al., 2015

to first condense 10 quantitative parameters in five independent predictors that were collectively analyzed to correlate with the microvascular density in carotid plaques. Finally, even for glioblastoma, the use of model-based features extracted from intra-operative DCE-US scans of the brain as input to an SVM classifier was reported to be a feasible approach to highlight tumor borders (Ritschel et al., 2015

).Aside from the parameters extracted through (pharmaco)kinetic modeling, radiomic features are increasingly used in the classification process. The field of radiomics comprises a range of spatial, statistical and structural features based on the original imaging or one of the extracted parameters (

Avanzo et al., 2017

). The underlying idea of radiomics is that when properly combined, meaningful information on texture, function and morphology can be captured that would normally remain invisible to the naked eye. Exploiting DCE-US discrete wavelet transform and texture features, for example, Acharya et al., 2011

found that a - Acharya U.R.
- Faust O.
- Sree S.V.
- Molinari F.
- Garberoglio R.
- Suri J.S.

Cost-effective and non-invasive automated benign & malignant thyroid lesion classification in 3D contrast-enhanced ultrasound using combination of wavelets and textures: A class of ThyroScan algorithms.

*Technol Cancer Res Treat.*2011; 10: 371-380

*k*-NN classifier had high performance in distinguishing benign and malignant thyroid nodules. Likewise,Theek et al., 2018

studied more than 200 DCE-US radiomics for the classification of tumor phenotypes in xenograft models.With the introduction of more powerful deep-learning approaches, feature extraction is today more and more “learned” during the training process, alleviating the need for explicit feature selection. Examples are convolutional encoder–decoder and auto-encoder networks that are designed to extract a hierarchy of higher-level features useful for the classification task.

Wu et al., 2014

applied deep learning for TIC classification, outperforming classic machine-learning approaches to the classification of focal liver lesions. Also for the prostate, despite using mice models, Feng et al., 2018

reported that a 3-D convolutional neural network might be able to predict prostate cancer using targeted DCE-US recordings with reasonable accuracy.The outcome of machine-learning technology is either (multi-)class classification or a related multiparametric “score” that reflects likelihood of classification. Examples of these scores are the Thyroid Malignancy Index proposed by

Acharya et al., 2011

or the classification confidence in GMMs (- Acharya U.R.
- Faust O.
- Sree S.V.
- Molinari F.
- Garberoglio R.
- Suri J.S.

Cost-effective and non-invasive automated benign & malignant thyroid lesion classification in 3D contrast-enhanced ultrasound using combination of wavelets and textures: A class of ThyroScan algorithms.

*Technol Cancer Res Treat.*2011; 10: 371-380

Wildeboer et al., 2017

). Mostly, however, machine-learning approaches are developed as computer-aided diagnostic tools appending conventional DCE-US imaging.## Discussion and recommendations

Over the years, several methods have been developed for the quantification of hemodynamic parameters based on the combination of DCE-US and pharmacokinetic modeling. Global parameters representing the central circulation, as well as local parameters reflecting tissue perfusion, have been proposed and employed for a large spectrum of diagnostic applications, ranging from myocardial perfusion defects to cancer angiogenesis.

The main challenge in the application of indicator dilution theory to DCE-US relates to calibration, that is, the establishment of the relationship between ultrasound intensity and UCA concentration. In an ideal world, ultrasound system manufacturers would provide data complying with a quantification-compatibility standard. With such a technical standard, clinicians, researchers and software solution developers could be provided with an operational mode and data paths in ultrasound systems suitable for performing perfusion quantification in a robust and reproducible manner. Such a standardized operational mode may represent a valuable step toward making DCE-US a truly functional, widely used and quantifiable imaging modality.

Currently, also because of the inherent limitations related to the inhomogeneous pressure field, absolute calibration is not feasible, and quantitative DCE-US imaging relies on linearized data. As a result, use of the mass conservation law and Stewart–Hamilton equation (eqn [1]) is unfeasible, and researchers focused on the assessment of temporal dynamics. The use of system identification methods, in both the temporal (impulse response) and spatiotemporal (Green's function) domains, has paved the way for important developments in quantitative DCE-US, enabling the assessment of local hemodynamic parameters.

Temporal analysis by the infusion protocol would in principle permit quantification in multiple planes with a single bolus injection and standard 2-D acquisitions; however, it is more challenging to implement in practice, as it requires the use of an injection pump, a dedicated imaging sequence and the application of a high-pressure burst; moreover, a steady-state plateau can be rather difficult to reach (

Dietrich et al., 2012

). In comparison, the bolus injection method offers a much simpler protocol and has been adopted more widely. Several models have been proposed to describe the bolus kinetics of conventional and targeted UCAs. In general, models that are based on the physics and physiology underlying the transport and distribution of UCA particles should be preferred, as they enable the estimation of interpretable parameters, which can often be directly related to the pathophysiologic processes of interest. Spatiotemporal analysis, which exploits information on UCA kinetics in both the spatial and temporal domains, may further overcome the limitations of the analysis in one domain only (be it temporal or spatial), by which the information present in the other dimension is disregarded or often reduced by averaging.As DCE-US recording can take up to several minutes, several sources of intrinsic (

*e.g.,*breathing) and extrinsic (*e.g.,*probe tilting) motion unavoidably corrupt DCE-US data. As a result, motion-compensation strategies should be employed prior to quantitative analysis. In 2-D DCE-US, compensation of motion perpendicular to the imaging plane is rather challenging and often limited to identification and exclusion of out-of-plane frames (Schäfer et al., 2015

); on the other hand, several methods are available to compensate for in-plane motion in 2-D DCE-US, which are typically based on image registration (Akkus et al., 2012

; Bouhlel et al., 2014

; Carvalho et al., 2014

; Schäfer et al., 2015

). In addition to enabling in-plane motion compensation, registration of DCE-US data should also be performed to ensure accurate clinical validation, for which matching of quantitative analysis with a ground truth is often required (Wildeboer et al., 2018c

). Moreover, image registration represents an important step for accurate training and validation of machine learning models, whereby data sets from different ultrasound acquisitions and modalities are combined (Mannaerts et al., 2019

; - Mannaerts C.K.
- Wildeboer R.R.
- Remmers S.
- van Kollenburg R.A.A.
- Kajtazovic A.
- Hagemann J.
- Postema A.W.
- van Sloun R.J.G.
- Roobol M.
- Tilki D.
- Mischi M.
- Wijkstra H.
- Salomon G.

Multiparametric ultrasound for prostate cancer detection and localization: Correlation of B-mode, shearwave elastography and contrast-enhanced ultrasound with radical prostatectomy specimens.

*J Urol.*2019; 202: 1166-1173

Wildeboer et al., 2019

). Although lower frame rates may limit tracking of fast kinetics, 3-D DCE-US provides important opportunities to overcome the current limitations imposed by out-of-plane motion and speckle decorrelation, enabling more accurate registration and motion compensation; thus, 3-D acquisitions should be the preferred option when available.Both for 2-D and especially for 3-D DCE-US, the signal-to-noise ratio may be poor, and pre-processing of the images is required prior to the analysis. Developments in singular-value-decomposition filtering are showing promise in this respect, enabling the enhancement of those components that are related to the UCA kinetics while suppressing clutter (

Wildeboer et al., 2019

). This method has also proven valuable in reducing motion artifacts, which can be present especially in the course of long recordings (Desailly et al., 2017

). In addition to the employment of standard motion compensation based on speckle tracking, improved signal quality was also obtained by employment of TIC fitting quality as cost function for optimizing motion compensation (Barrois et al., 2015

). Spatial and spatiotemporal analysis also require dedicated strategies to account for depth-dependent and anisotropic resolution. This can be achieved by the employment of spatial regularization or adaptive kernel sizes (Schalk et al., 2015

; Wildeboer et al., 2019

).## Conclusions and future perspectives

When proper acquisition and calibration are performed, advanced pharmacokinetic modeling has established its value for assessment of relevant hemodynamic parameters by quantitative DCE-US.

Today, ultrasound systems allowing 3-D DCE-US are increasingly available for a variety of 3-D ultrasound probes (

Huang and Zeng, 2017

). From the perspective of quantification, 3-D acquisitions have the advantage of capturing the volumetric contrast behavior, alleviating the need for multiple acquisitions to cover a 3-D field-of-view and more readily allowing motion compensation. Moreover, given the 3-D nature of the underlying dilution systems, 3-D DCE-US may lead to more accurate hemodynamic assessment, accounting for full boundary conditions and 3-D flow directions. On the other hand, the acquisitions are computationally more demanding and have a lower spatiotemporal resolution (Schalk et al., 2015

). Schalk et al., 2015

reported that relevant bolus kinetics in the prostate is still captured at frame rates of ∼0.25 Hz, but this might differ among applications. Indeed, in both the prostate (Schalk et al., 2015

) and the liver (Cao et al., 2019

), 3-D quantitative features were in good agreement with the corresponding values in 2-D, although the estimation of 3-D features in deep regions seems somewhat limited. In general, these results already evidence the feasibility of quantitative 3-D DCE-US. Moreover, ultrafast ultrasound imaging methods that enable DCE-US imaging at thousands of frames per second are rapidly being developed (Montaldo et al., 2009

; Couture et al., 2012

; Viti et al., 2016

), providing a path to overcome the limitations in spatiotemporal sampling; however, the implementation of a 3-D ultrafast imaging system poses important hardware requirements. Recent developments in ultrafast imaging, in combination with sparse sampling, are making important advances toward efficient implementation of ultrafast 3-D DCE-US imaging (Roux et al., 2016

).To infer the underlying (micro)vascular architecture, usually spatial analysis is applied to parametric maps extracted from DCE-US loops. However, ultrasound localization microscopy is also being developed that permits ultrasound imaging at a virtual resolution down to the microbubble size (