Review ULTRASOUND CONTRAST AGENTMODELING : A REVIEW

* Physics of Fluids Group, MESA+ Institute for Nanotechnology, Technical Medical (TechMed) Center, University of Twente, Enschede, the Netherlands; yDepartment of Engineering Science, Institute of Biomedical Engineering, University of Oxford, Oxford, UK; and zCentre National de la Recherche Scientifique (CNRS), Laboratoire Interdisciplinaire de Physique (LIPhy), Universit e Grenoble Alpes, Grenoble, France


INTRODUCTION
While in the 17th century Robert Hooke had already predicted that we would be able to image the human body with sound (Tyndall 1875), it was not until the 1940s that the first ultrasound scan of soft tissue was realized (Howry 1952). Today, ultrasound is one of the most widely used medical imaging techniques because of its safety and relatively low expense compared with magnetic resonance imaging and X-ray computed tomography. As with the other modalities, ultrasound allows for functional imaging, including Doppler (Evans et al. 1989;Christopher et al. 1996), vector flow imaging (Jensen et al. 2016) and neuroimaging (Deffieux et al. 2018). Moreover, ultrasound machines are small and portable and can be used at the bedside. Finally, ultrasound provides real-time imaging, offering typically 50 frames per second for conventional line-by-line scanning and up to 10,000 frames/s for ultrafast plane wave imaging (Tanter and Fink 2014).
The best-known application of medical ultrasound is fetal scanning. The amniotic fluid around the fetus contains very few scatterers in comparison to the fetal tissue, and as a result the image of the fetus stands out clearly against a black background. Such high contrast can also be found in echocardiography, where the blood in the left ventricular cavity is a poor ultrasound scatterer as compared with the cardiac muscle tissue. This high contrast then aids in the delineation of the endocardial border and the diagnosis of regional wall motion abnormalities (Mulvagh et al. 2008;Senior et al. 2009). However, the poor backscatter from the blood pool is problematic for perfusion imaging and blood flow quantification.
Ultrasound contrast agents (UCAs) can dramatically boost the native contrast of the blood pool. A UCA consists of a suspension of stabilized gas microbubbles (Quaia 2005) that is injected intravenously. As the microbubbles are smaller than red blood cells, they can traverse even the smallest capillaries. Despite their diminutive size, their scattering cross section can be up to nine orders of magnitude larger than that of red blood cells, and this highly efficient scattering enables the quantification of the perfusion of the myocardium and other organs, including liver, brain and kidney.
The potential of microbubbles as UCAs was in fact discovered by accident in 1968, during an intravenous injection of a saline solution while performing a cardiac ultrasound examination (Gramiak and Shah 1968). The microbubbles that were formed in solution by hand agitation were found to scatter the administered ultrasound highly efficiently. Given the poor stability of the airfilled bubbles in the electrolyte solution, considerable effort was put into finding more stable alternatives. This led to the first class of agents-including Quantison (Frinking and de Jong 1998), Albunex (Feinstein et al. 1984) and Optison (Skyba et al. 1996;Cohen et al. 1998)-containing coated microbubbles that were formed by sonication of human serum albumin. With their relatively rigid shells and corresponding decreased scattering response, a second generation of UCAs was developed that was inspired by the more flexible nature of lipid cell membranes. These novel agents-including SonoVue (Schneider 1999), Definity (Unger et al. 1994) and Sonazoid (Sontum 2008)-are composed of a suspension of purified phospholipid-mixtureÀencapsulated gas microbubbles, with diameters ranging from 1 to 10 mm, that are filled with a high-molecular-weight gas, typically perfluorocarbons or sulfur hexafluoride. Examples of UCA microbubbles are given in Figure 1aÀd.
The coating of the microbubbles decreases (or potentially eliminates) the interfacial tension, which stabilizes the bubble against Laplace pressureÀdriven dissolution (Epstein and Plesset 1950). The stability of the microbubbles is further enhanced by both the low solubility of the gas and the resistive barrier that the encapsulation forms against gas diffusion through the monolayer membrane (Borden and Longo 2002). Long-chain polyethylene glycol spacers are added to provide stealth capability against detection by the body's immune system and therefore improve the circulation lifetime in vivo (Klibanov 2006). Furthermore, they prevent bubble coalescence through repulsive intermolecular forces (Segers et al. 2017;Borden and Song 2018). Cumulatively, this leads to an overall bubble lifetime of 5À10 mins after intravenous injection (Kaul 2000;Lindner 2019). Perfluorocarbon gas-filled microbubbles can be stable in suspension for several days to a few weeks when kept in a vial under a protective gas headspace. Freeze-dried bubbles can be stored for several years (European Medicines Agency 2015).
The physical mechanism by which microbubbles enhance image contrast is twofold. First, acoustic waves are scattered by the bubbles owing to the large difference in acoustic impedance between the gas core and the surrounding liquid. Second, and more important, the large compressibility of the gas core leads to volumetric oscillations in response to ultrasound pressure waves. Thus, the oscillating bubble acts as a source of sound with an effective scattering cross section up to several orders of magnitude larger than that of a rigid sphere of the same size and acoustic impedance (Hilgenfeldt et al. 1998). Furthermore, when the bubble s are driven at resonance, their oscillations can be highly non-linear and therefore they can produce scattered ultrasound fields containing a range of frequencies, including whole and fractional harmonics of the driving frequency. In addition, and particularly at lower driving pressures (<100 kPa), further non-linearities are induced by the non-linear viscoelastic properties of the stabilizing shell. To date, contrast-specific imaging modes use the non-linear echoes to detect the contrast agent response while suppressing the linear tissue echoes (Hope Simpson et al. 1999;Frinking et al. 2000;Rafter et al. 2004;Eckersley et al. 2005;Averkiou et al. 2008). By fortunate coincidence, the range of microbubble sizes appropriate for clinical use corresponds to resonance frequencies on the order of 1À10 MHz, the typical frequency range used in diagnostic ultrasound. It is this combination of scattering efficiency and non-linearity that makes microbubbles such excellent UCAs.
While the use of coated bubbles as a perfusion contrast agent is clinically established, the interest in oscillating bubbles as a theranostic agent is also rapidly increasing for the following reasons. First, transfer of momentum to the surrounding fluid by the non-linearly oscillating bubble can set up a steady flow, or microstreaming, that can both affect the transport of therapeutic material and impose shear stresses on biological membranes, thereby altering their permeability and thus mediating therapeutic effects (Ward et al. 2000;Marmottant and Hilgenfeldt 2003;Pereno et al. 2018). Second, absorption of the acoustic emissions from the bubbles can produce localized heating of the surroundings (Prosperetti 1977b;Coussios and Roy 2008) that can directly produce bio-effects or be used, e.g. to release drugs from thermally sensitive carriers. Heating is increased compared with direct ultrasound exposure on account of the higher-frequency components in the bubble emissions that are more readily absorbed by tissue (Hilgenfeldt et al. 2000). Third, the increased temperature and pressure in the core of the bubble during compression can stimulate sonochemistry (Suslick and Flannigan 2008), leading to the production of chemical species including oxygen radicals, although this effect is reduced for the polyatomic gases typically used in contrast agents. Finally, there may be complex interactions between oscillating bubbles and biological structures that promote therapeutic effects. These include microjetting and radiation force phenomena.
In the mid-1990s, a new field emerged that focused on the complex interaction between coated microbubbles and ultrasound. Advances were made through experimental observations, first using acoustic spectroscopy (de Jong et al. 1992) and later using ultra high-speed imaging (Morgan et al. 2000;Chin et al. 2003). Along with the experimental work, great progress has been made in model development. The majority of the available coated-bubble models have in common a starting point of a Rayleigh PlessetÀtype equation (Lauterborn 1976;Matula 1999). The field of UCA modeling has greatly benefited from the field of sonoluminescence, where the dynamics of free gas bubbles in an unbounded fluid have been fully described (Brenner et al. 2002;Lohse 2018). The effect of the viscoelastic shell on bubble dynamics is evident from acoustical and optical characterization. The pressure contributions resulting from the viscoelastic shell are incorporated into the RayleighÀPlesset (RP) equation as additional interfacial pressure terms (de Jong et al. 1994;Church 1995;Hoff et al. 2000). The shell elasticity increases the stiffness of the system, and an increase in the resonance frequency is expected. Similarly, the shell viscosity leads to energy dissipation which results in increased damping and hence a small reduction in resonance frequency and decreased amplitude of oscillation.
The non-linear properties of the bubble shell are especially important at small oscillation amplitudes. The  Segers et al. 2016b. Copyright 2016 American Chemical Society); (b) perfluoropropane bubbles with a protein shell (Optison), seen here against a background of red blood cells (reproduced from Wilson and Burns 2010); (c) lipid-coated microbubbles of perfluoropropane gas (Definity) as seen in dark-field microscopy (reproduced from Wilson and Burns 2010); (d) scanning electron micrograph of a Quantison microsphere, fractured owing to vacuum (reproduced with permission from Frinking and de Jong 1998); (e) vapor bubble release from a laser-heated single oil-filled microcapsule captured under ultrahigh-speed imaging (reproduced with permission from Lajoinie et al. 2014); (f) scanning electron micrograph of monodisperse microcapsules (image courtesy of Nanomi, Oldenzaal, the Netherlands); (g) monodisperse microbubbles with a radius of 2.7 mm produced in a flow-focusing device (reprinted with permission from Segers et al. 2016b. Copyright 2016 American Chemical Society); (h) numerical simulation of the buckling instability of the elastic shell of a polymeric hard-shelled microcapsule (reproduced with permission from Marmottant et al. 2011); (i) buckling instability upon dissolution for two different surfactants (Survanta and dipalmitoylphosphatidylcholine; reprinted with permission from Thomas and Borden 2017. Copyright 2017 American Chemical Society); (j) buckling observed during compression of a lipid-coated microbubble (reproduced with permission from Sijl et al. 2011). densely packed lipid monolayer may buckle during compression, while the bubble coating may rupture upon expansion, forming lipid domains on the otherwise free interface and thereby exposing the gas directly to the surrounding liquid. Such changes in the interfacial properties of the coating leads to surprising contrast response characteristics, including thresholding behavior , compression-only behavior Sijl et al. 2011) and subharmonics ). In addition, gas loss through diffusion, or acoustic deflation, upon ultrasound driving may lead to a higher concentration of lipids at the bubble interface, changing its dynamics. Alternatively, the bubble may expel excess lipids (Luan et al. 2014), leading to higher surface tension and promoting gas loss (Viti et al. 2012). Altogether, the interfacial rheology is an important factor in fully describing the origin and details of the non-linear echoes generated by bubbles. Ultimate control over the shell properties and its physical behavior is elementary for the design of optimized agents for contrast imaging, displaying improved backscatter, less attenuation and superior non-linear properties. Such agents can then also act as vector flow tracers for echo particle image velocimetry (Leow et al. 2015;Poelma 2016;Engelhard et al. 2018) and enable non-invasive pressure measurements (Shi et al. 1999;Frinking et al. 2010). In addition to the dramatic effect of the stabilizing shell, bubble dynamics are highly dependent on the biological environment and the presence of neighboring bubbles. A compliant blood vessel wall will reflect part of the pressure field emitted by the bubble in a non-trivial way, resulting in a complex interaction with the bubbles (Garbin et al. 2007;Doinikov et al. 2009b;Hay et al. 2012;Helfield et al. 2014). Similarly, the sound field generated by neighboring bubbles interacts with that of the bubble, thereby altering its dynamics and inducing acoustic radiation forces (Bjerknes 1906;Crum 1975). These effects are considered higher-order effects and are discussed in more detail later under Modeling higher-order effects, as well as in the companion review on therapeutic applications of microbubbles (Kooiman et al. 2020).
This article reviews the modeling efforts that have been undertaken, starting with bubble fundamentals, bubble shell models and secondary higher-order effects, and finishing with new directions that are explored for improved models to describe bubble dynamics in tissues and those of novel agents such as bubble clusters and bubbles for multimodality imaging. One or more of our cited studies used human subjects, and these were, to the best of our knowledge, granted protocol approval by an ethics committee or institutional review board. Where animals were studied the research protocol was, again to the best of our knowledge, approved by a local institutional care and animal use committee.

Modeling free-bubble dynamics
The free-gas-bubble model. Numerous mathematical formulations have been developed to describe bubble dynamics. Historically these have been motivated primarily by the hydrodynamic effects of bubbles in turbomachinery. The first published use of the term cavitation dates to Thornycroft and Barnaby (1895) and is attributed to the engineer and naval architect Robert Froude. That bubbles might be generated by the motion of a water turbine had in fact been proposed over a century earlier by Euler (1756), but it was not until the late 19th century that the relationship between liquid pressure, bubble formation and, importantly, physical effects on the surroundings became the focus of intense research. Besant (1859) and subsequently Lamb (1879) and Rayleigh (1917) considered an empty spherical cavity in a large volume of inviscid fluid. In this case the bubble undergoes a single, very rapid, collapse, driven by the inertia of the surrounding fluid and generating a shock wave as a result of the high pressure developed at the site of the bubble. Blake (1949), Noltingk and Neppiras (1950) and Neppiras and Noltingk (1951) subsequently considered the more commonly encountered case of a bubble partially or completely filled with gas. In this case the pressure inside the bubble increases as it is compressed, producing a cushioning effect. Finally, an additional term was included to account for viscous dissipation in the surrounding fluid. According to Young (1989) this was originally given by Taylor and Davies (1963), although it is generally attributed to Poritsky (1951). The result was the Rayleigh-Plesset or Ray-leighÀPlessetÀNoltingkÀNeppirasÀPoritsky equation (Lauterborn 1976): This is a pressure balance between the inertial term on the left-hand side, originating from the inertia of the liquid associated with the bubble motion, and the pressure difference on the right-hand side arising from the difference in work done at the bubble wall and done remotely from the bubble (Leighton 2011). Table 1 summarizes all symbols used in this review. In this equation, R is the time-dependent bubble radius (see Fig. 2, left), with initial value R 0 ; _ R and € R are the velocity and acceleration, respectively, of the bubble wall; r L is the density of the surrounding fluid; P B is the pressure in the liquid at the bubble wall; and P 1 is the pressure in the liquid far away from the bubble, defined as the sum of the hydrostatic pressure P 0 and the time-varying driving pressure P A (t) in the liquid.
When the gas in the bubble is an ideal gas, from the pressure balance at the bubble interface P B can be Ultrasound in Medicine & Biology Volume 00, Number 00, 2020 expressed as following the ideal gas law, with s the surface tension at the gasÀliquid interface, k the polytropic constant and m L the viscosity of the surrounding liquid. The vapor pressure p v inside the bubble is typically neglected in the modeling of the dynamics of micrometer-sized bubbles because the vapor pressure of water at atmospheric conditions (2.3 kPa) is low as compared with the hydrostatic and interfacial pressure terms.
There are a number of assumptions implicit in the derivation of eqns (1) and (2): the bubble is assumed to remain spherical at all times, with no heat or mass transfer across the bubble wall during oscillation. The surrounding fluid is assumed to be infinite, incompressible and Newtonian. The gas inside the bubble is assumed to behave ideally. The wavelength of the incident ultrasound field is assumed to be sufficiently large that the bubble experiences a spatially uniform pressure. The effects of gravity and buoyancy are ignored.
Acoustic reradiation. At small amplitudes of oscillation the assumption of an incompressible fluid is justified. However, with increasing driving pressure it becomes necessary to modify eqn (1) to correct for the Fig. 2. Schematic view of a free gas bubble (left) and a lipidcoated microbubble (right). To the right, the coating is indicated by the primary phospholipids and lipids with polyethylene glycol added (not drawn to scale). R 0 = bubble size; P B = liquid pressure at the bubble wall; P 1 = liquid pressure far away from the bubble; m L = liquid viscosity; r L = liquid density; c L = speed of sound in the liquid; s = surface tension of water; k s = shell viscosity; x = shell elasticity; s(R) = sizedependent effective surface tension of the coated bubble. effects of liquid compressibility and energy dissipation associated with generating acoustic emissions. This may become significant for situations in which the velocity of the bubble wall is comparable with the speed of sound in the surrounding liquid ðj _ Rj=c L » 1Þ. Based on the work of Herring (1941), Trilling (1952) rederived eqn (1) to take into account the energy stored in the liquid during bubble oscillation and the damping effect produced by the re-radiation of sound waves. This treatment assumes that the sound velocity c L in the surrounding liquid is constant and is therefore limited in terms of the maximum amplitude of oscillation for which it is valid. Gilmore (1952) improved upon this assumption by employing the KirkwoodÀBethe hypothesis (Bethe 1935;Kirkwood 1938), which states that the wave velocity is equal to the sum of the velocity of the @Eleanor fluid or liquid? and the local speed of sound at the bubble wall.
In 1980, building upon the earlier models, Keller and Miksis (1980) published a more complete treatment for large-amplitude bubble oscillations, which included surface tension and liquid viscosity, both of which were neglected in earlier work: with P B and P 1 as defined before. This expression was subsequently refined by Prosperetti et al. (1988), who removed the assumption of ideal gas behavior and coupled the equation of motion for the bubble to the equations of mass and energy conservation. Note that eqn (3) is no longer valid when j _ Rj=c L » 1. Therefore, a popular form in the context of sonoluminescence and UCA modeling is given by (L€ ofstedt et al. 1995;Barber et al. 1997;Brenner et al. 2002;Marmottant et al. 2005;Paul et al. 2010): For an exhaustive review of the derivation of the different RP-type equations, see Brenner et al. (2002).

Modeling coated-bubble dynamics
Adding the interfacial pressure terms. The gas core of UCA microbubbles is typically coated with a layer of proteins, polymers or phospholipids (Stride and Saffari 2003;Faez et al. 2013). The first class of agents comprised microbubbles with a rigid shell, often referred to as hardshelled agents, which motivated the initial approach to modeling their response to ultrasound. This work was started by de Jong and Hoff, who proposed a model describing the behavior of microbubbles of Albunex contrast agent, where the shell was considered to be perfectly elastic (i.e., obeying Hooke's law with no additional damping [de Jong et al. 1992]). The resulting stiffness pressure term was added in an ad hoc way to the RP equation to describe the bubble dynamics. They later added a semiempirical shell dissipation term to fit the experimental observations (de Jong and Hoff 1993). Church (1995) and Hoff et al. (2000) then followed a similar approach using a more thorough derivation to account for the viscoelastic shell of Albunex. They considered the shell material to be incompressible, but the shear viscosity of the shell was directly included in the derivation of the modified RP equation. Several other models were developed based on linear elasticity (some of them are compared by Sarkar et al. 2005), presenting similar results. Non-linear elasticity considerations were added by Allen and Rashid (2004), who considered a non-linear Hooke's law, or neo-Hookean model, for the shell. However, all these models missed the most important aspects of the behavior of the rigid shells, namely (non-spherical) buckling upon compression and rupture upon large expansion. These aspects were investigated and implemented in a model by Marmottant et al. (2011). Despite its simple approach, this model captures many features observed in experiments (Lensen et al. 2011). The model shows, in particular, that these bubbles mostly behave as a free gas bubble after rupture of the shell, as demonstrated before by Postema et al. (2005a) and Bouakaz et al. (2005). A more extensive review of the hard-shelled models is given by .
To date, the majority of the clinically approved UCAs and research-grade UCAs consist of microbubbles coated with a phospholipid monolayer. In the following, we therefore focus on the modeling of non-linear lipidcoated microbubble dynamics-e.g., the observed driving-pressureÀdependent resonance behavior (Chen et al. 2002;Emmer et al. 2009;Segers et al. 2016a) and asymmetric oscillation amplitudes-with more pronounced compression than expansion, termed compression-only behavior Sijl et al. 2011). Figure 2 (right) gives a schematic of the coated bubble. The viscoelastic lipid coating changes the pressure jump condition across the bubble interface that is used in the RP equation, eqns (1) and (3). The seminal articles by de Jong et al. (1992Jong et al. ( , 1994 and Church (1995) include the presence of a shell through the addition of two independent pressure terms on the right-hand side of eqn (4), one accounting for the elasticity P elas of the shell and one accounting for the viscous dissipation P vis in the shell. These terms can be derived for a shell of finite thickness (Church 1995;Hoff et al. 2000) and for an infinitesimally thin surface layer using material constitutive laws (de Jong et al. 1992;de Jong and Hoff 1993;de Jong et al. 1994;Chatterjee and Sarkar 2003;Marmottant et al. 2005;Sarkar et al. 2005;Tsiglifis and Pelekasis 2008;Doinikov et al. 2009a). The physical principles of the available models for soft-shelled models have been exhaustively reviewed by . For a more extensive review, see Faez et al. (2013).
The viscous pressure term P vis is added to account for viscous energy dissipation in the microbubble shell resulting from the work done to change the distance between the lipid molecules during an acoustic cycle. The most common form of P vis is derived from Newton's linear viscous law for an infinitesimally thin shell (de Jong et al. 1994;Church 1995;Morgan et al. 2000;Chatterjee and Sarkar 2003;Marmottant et al. 2005;Sarkar et al. 2005): with k s the surface dilatational viscosity of the lipid monolayer shell. Typically, it is assumed that k s does not depend on the bubble surface area, surface dilatation, or surface dilatation rate. The validity of eqn (5) will be discussed in the final paragraph of this section. The elastic pressure term P elas accounts for the pressure increase inside the bubble owing to the stiffness of the shell. Shell elasticity x can be described to result from interfacial tension gradients ds induced by (acoustically-driven) bubble surface-area variations dA (Glazman 1983;Mohwald 1995;Marmottant et al. 2005;Sarkar et al. 2005): Assuming a constant and surface-areaÀindependent shell elasticity, eqn (6) can be integrated to find the effective interfacial tension: sðRÞ ¼ xðR 2 =R 2 0 À1Þ, with R 0 the equilibrium bubble radius. With the effective interfacial tension, the elastic pressure term reads Equation (7) is the most common form of the elastic pressure term (de Jong et al. 1994;Church 1995;Chatterjee and Sarkar 2003;Marmottant et al. 2005;Sarkar et al. 2005;Tsiglifis and Pelekasis 2008). The full equation-including the shell pressure terms and eqns (5) and (7)-then reads: where s(R 0 ) represents the surface tension of the bubble at rest. With a complete description of the radial oscillations of coated bubbles, we can now calculate the echo produced by UCA microbubbles. The sound emitted by an oscillating microbubble has two contributions. The first is a passive contribution that simply results from the mismatch in acoustic impedance between the gas core and the surrounding medium. The second one is the active part and results from the volumetric oscillations V(t) of the bubble. Micrometer-sized bubbles driven near resonance (see also the next section) have a size that is typically three orders of magnitude smaller than the incident wavelength. It has been shown that in this case the passive contribution is negligible with respect to the active part (Hilgenfeldt et al. 1998). From mass and momentum conservation it follows that the active scattered pressure wave P s at a distance r from the bubble center equals (Vokurka 1985) Thus, the scattered pressure P s (t) can be directly calculated by inserting into eqn (9) the solution of the RP equation, R(t), and its derivatives.
The numerical solution R(t) of the ordinary differential equation (ODE) given by eqn (8) is displayed in Figure 3 (left) for a bubble with radius 2.0 mm. Three driving pressures were used here-5, 50 and 500 kPacorresponding to bubble oscillation behavior in the linear, non-linear and inertial-collapse regimes, respectively. The scattered pressure P s (t) was calculated by inserting R(t) and its derivatives _ R and € R directly into eqn (9) (see Fig. 3, middle). Figure 3 (right) gives the spectral content of the scattered pressure. It nicely shows the different regimes of bubble oscillation behavior: in the linear regime, the response of the bubble is at the fundamental frequency of 2.0 MHz; for a driving pressure of 50 kPa the response is non-linear, with buildup of second harmonic and superharmonic components and some small amount of broadband emission; and for the highest-pressure regime, despite some harmonic content, the response is primarily broadband owing to the impulsive collapses of the bubble. Figure 4a shows the maximum bubble response, expressed as the normalized difference of maximum and minimum bubble amplitude, as a function of the bubble radius for a driving frequency of 1.0 MHz and a driving pressure of 75 kPa. Bubble shell and liquid properties are the same as in Figure 3. The resonance behavior is clearly visible, with bubbles with a radius of 3.2 mm displaying maximum oscillation amplitudes. Figure 4B shows the fundamental and second harmonic components of the scattered response of the bubbles driven in Figure 4a. It shows first of all that next to the resonance behavior, bubbles of a larger size also contribute to the scattering at the fundamental driving frequency. More importantly, it shows that the non-linear generation of the second harmonic is limited to the resonant bubbles.
Linear bubble-shell modeling: Resonance frequency and damping. The role of the shell in the microbubble dynamics becomes clear through linearization of the RP equation (eqn [8]), by assuming sinusoidal acoustic driving and small amplitudes of oscillation with respect to the equilibrium radius (i.e., R ¼ R 0 1 þ 2 ðtÞ , with e (t) ( 1), and by neglecting higher-order terms to yield the classic equation of a driven harmonic oscillator: with _ e and € e the first-and second-order derivatives of e with respect to time, d tot the total damping of the system, v the angular ultrasound frequency and v 0 the angular eigenfrequency of the oscillator: In these equations, s(R 0 ) represents the surface tension of the bubble at rest. Equation (11) shows that the eigenfrequency of a bubble is inversely proportional to the microbubble size and that the stiffness x of the coating increases the eigenfrequency of the coated bubble compared with the free gas bubble.
The total damping d tot is a summation of the damping contributions that arise from reradiation of sound, viscous dissipation in the surrounding liquid and shell viscous dissipation: The linearization of eqn (8) results in the following expressions for the damping terms at v = v 0 : This expression for radiation damping is in approximate agreement with the exact form d rad ¼ vR 0 =c L (Leighton 1994). A microbubble also experiences thermal damping d th owing to heat diffusion (Prosperetti 1976;Leighton 1994). For micrometer-sized bubbles it is of the same order as the viscous damping for a bubble in water (Devin 1959), and therefore typically the viscosity of the medium is increased phenomenologically by a factor of 2 to account for thermal damping Sijl et al. 2010;Segers et al. 2016a). The damping constants for an uncoated bubble are plotted in Figure 4c. Figure 4d shows the damping constants for a coated bubble (x = 0.5 N/m). Shell damping was calculated using a fixed shell viscosity k s = 1.0 £ 10 À8 kg/s. Note that the damping for a coated bubble is dominated by the shell viscous damping, as it contributes approximately 75% of the total damping for bubbles in the relevant size range of 2À3 mm. The total damping d tot of the system results in a slightly reduced microbubble resonance frequency f res as compared with q . For a typical total damping of 0.25, however, the decrease in resonance frequency is less than 2%. On the other hand, eqn (11) demonstrates that the microbubble shell stiffness can dramatically increase the resonance frequency. Figure 5 compares the resonance curve of a free gas bubble (R 0 = 2.5 mm) with that of lipid-coated bubbles of the same sizes, with respective shell stiffness of 1.0 and 2.5 N/m and shell viscosities of k s = 6 £ 10 À9 and 12 £ 10 À9 kg/s. Note that the linearization performed here provides insightful information only for coated microbubbles driven at very low pressure amplitudes (<5 kPa). The total damping of a coated bubble (x = 0.5 N/m) over the size range relevant for ultrasound contrast agents. The shell damping d shell was calculated using a fixed shell viscosity k s = 1.0 £ 10 À8 kg/s. Note that the total damping of a coated bubble is dominated by shell damping for bubbles in the diameter range 2À3 mm.
Non-linear bubble-shell modeling: Buckling and rupture. For elevated acoustic pressures we will need to describe the non-linear dynamics of the coated bubble, which is then modeled by solving the full RP equation including the shell pressure terms (eqn [8]). The RP equation is non-linear owing to the intrinsic non-linearities in the ODE. However, the shell pressure terms are linear (i.e., the surface tension is linearly dependent on the microbubble surface area) and can adopt negative values or values exceeding the surface tension of the surrounding medium (see the black dotted line in Fig. 6a), and the shell stiffness x is constant throughout and independent of the surface dilatation (see Fig. 6b). Therefore, eqn 8 does not predict the experimentally observed non-linear microbubble dynamics such as compression-only behavior ) and the resonance-frequency shift with increasing acoustic pressure amplitudes for lipid-coated bubbles .
The breakthrough in modeling non-linear bubble dynamics came from Marmottant et al. (2005). Inspired by Langmuir-trough surface pressure-area isotherms of lipid monolayers (Fig. 6c), they limited the surface tension to zero for a compression below a certain buckling radius R b and to the surface tension of the surrounding medium for an expansion above a so-called rupture radius R r : The main idea is that during the compression phase, the lipid monolayer cannot withstand compression beyond the buckling radius, at which the shell buckles, resulting in a tensionless shell ðx ¼ Ads=dA ¼ 0Þ; see Figure 1j, where buckling is captured by ultrahigh-speed imaging. During the expansion phase, the average area per lipid molecule increases, either owing to shell rupture (e.g., along lipid domains ]) or owing to a more uniform increase in the area per lipid molecule, resulting in an increase of the effective surface tension until the surface tension s water of the surrounding medium is reached, again leading to a tensionless state ðx ¼ Ads=dA ¼ 0Þ. In the elastic phase, the surface tension was assumed to increase monotonically with surface area, as in eqn (7). The solid black line in Figure 6a shows the surface-tension curve used in the Marmottant model . Note that purely elastic (linear) oscillations are limited to small amplitudes of oscillation around the initial bubble surface area with radius R 0 . Figure 6b shows the corresponding elasticity x as a function of the surface area.
When the bubble is driven slightly outside the elastic regime, there are three possible oscillation scenarios, depending on the initial surface tension s(R 0 ) of the bubble. At s(R 0 ) < s water /2, the bubble can reach the buckling state during compression while it remains in the elastic state during expansion. During the buckling phase, the shell stiffness vanishes and the bubble therefore experiences less resistance against compression than against expansion. This asymmetry in shell stiffness over an oscillation period induces compression-only behavior, as shown in Figure 6d. Similarly, at s(R 0 ) > s water /2, the bubble experiences less resistance upon expansion than during compression once it is driven beyond the rupture radius. This drives the asymmetric expansion-only behavior (Fig. 6d). At sðR 0 Þ ¼ s water =2, the bubble oscillations are symmetric.
The non-linear interfacial tension and corresponding shell stiffness have another implication for a coated bubble driven outside its elastic regime, namely a downshift of its resonance frequency. The downshift results from the effective averaging of the shell stiffness over elastic and non-elastic regimes. When the acoustic pressure is increased, the effective shell stiffness decreases. At high driving pressures, the contribution of shell stiffness becomes negligible and the resonance frequency approaches that of an uncoated bubble. This sudden increase in microbubble response above an acoustic pressure threshold has also been named thresholding behavior . A detailed explanation of thresholding behavior is given by Overvelde et al. (2010). Next to compression-only and driving-pressur-eÀdependent resonance behavior, the Marmottant model has also proven successful in modeling the generation of subharmonics (Frinking et al. 2010;Sijl et al. 2010). It has been shown that subharmonics are induced by rapid change in shell stiffness with a small change in surface area (Sijl et al. , 2011. Thus, the second-order derivative of the interfacial tension curve with respect to the bubble surface area is the driving factor behind the generation of subharmonics. However, Marmottant et al. (2005) implemented the dilatational interfacial tension and corresponding shell stiffness in a rather unphysical way, as they are discontinuous at the buckling and rupture radii (Fig. 6a, 6b). Therefore, the numerical solutions of eqns (8) and (13) are sensitive to the integration time step of the ODE solver, as this sets d 2 s/dA 2 . To circumvent this problem, Sijl et al. (2010) proposed extending the Marmottant model with two crossover functions for the transition from the regions with zero elasticity to the one with finite elasticity (see the blue curves in Fig. 6a, 6b). With the aim of getting closer to a more physicochemical surface-areaÀdependent surface tension (i.e., with decreasing elasticity for increasing area per lipid molecule) Paul et al. (2010) tested a quadratic elasticity model and an exponential elasticity model (see Fig. 6a, 6b). They once again pointed out the importance of a surface-areaÀdependent elasticity for non-linear bubble dynamics.
It has now become clear that the actual shape of the surface-tension curve and the initial surface tension are important for the correct prediction of non-linear radial bubble dynamics. In the modeling efforts, ad hoc assumptions have been made regarding the rheological properties of the microbubble shell, as measurement of the actual s(R) dependence of a lipid-coated microbubble interface is difficult at megahertz frequencies (Marmottant et al. 2005;Stride 2008). The difficulty lies in the extremely short timescales and small length scales of ultrasound-driven microbubble oscillations, which render unuseful the relatively slow measurements obtained by Langmuir trough (Lozano and Longo 2009) or bubble tensiometer (Lee et al. 2001). Furthermore, UCA bubbles are at least one order of magnitude smaller than the minimum cap radius in a bubble tensiometer (Alvarez et al. 2010), and the curvature of a monomolecular film may play a role in its viscoelastic properties (Kotula and Anna 2016). Recently, the dilatational shell elasticity curve of lipid-coated bubbles has been measured experimentally by controlling the surface area of monodisperse bubbles through the ambient pressure and measuring the resonance behavior from acoustic attenuation measurements (Segers et al. 2018a). The measured xðA=A 0 Þ curve, Figure 6f, has been integrated to find the sðA=A 0 Þ curve, Figure 6e. The measured dilatational interfacial tension curve shows remarkable agreement with the initial assumptions of (Marmottant 2005), explaining the effectiveness of their model in the prediction of non-linear bubble dynamics. Is has also been shown that the x A=A 0 ð Þcurve is independent of the total bubble surface area for a given lipid composition (Segers et al. 2018a); in other words, it is independent of the bubble size. However, it is also shown that the shell stiffness is not constant in the elastic regime (see Fig. 6f), which may result in a change in acoustic properties of a dilute bubble suspension over time through gas loss.

Modeling shell viscosity
In contrast to the bubble-sizeÀindependent shell stiffness, the reported shell viscosity values are strongly dependent on the equilibrium bubble size. Shell viscosity k s has been extensively characterized through a direct fit of eqn (8) to experimental R(t) curves, through optical light scattering (Tu et al. 2009) and through acoustic and high-speed imaging characterization (Morgan et al. 2000;van der Meer et al. 2007;Overvelde et al. 2010;Segers et al. 2018a). The data show an exponential increase of k s with bubble size, which indicates an incomplete implementation of shell viscous dissipation in eqn (5). Shell viscosity has been shown to decrease with dilatation rate ), which indicates shear-thinning behavior. Based on those experimental data, Doinikov et al. (2009a) included shearthinning behavior in the RP equation in an empirical way. Furthermore, Segers et al. (2018a) have measured a dependence of k s on the bubble surface area (lipid-packing density), indicating that not only the dilatation rate is important but also surface dilatation itself. The future challenge is to clarify this matter, including in terms of the role of different shell components on the acoustic bubble response. Attempts have been made to model shell parameters based on molecular interactions (Stride 2008;Borden 2019). This route can be further explored, for example using molecular-dynamics simulations (Boek et al. 2005; van Opheusden and Molenaar 2011). Control over the bubble-size distribution has now been achieved ( Fig. 1g; Segers et al. 2016bSegers et al. , 2017Segers et al. , 2019. Together with monodispersity, control over the acoustic response of coated bubbles would be highly valuable, as this would allow the tailoring of microbubble design for specific clinical applications, e.g.non-invasive bloodpressure measurements, and therapeutic applications [Miller et al. 2018], including sonothrombolysis and drug and gene delivery using bubbles and ultrasound.

Scattering and attenuation
Scattering leads to a loss of energy of the insonifying acoustic wave. The acoustic energy that is lost by the scattering of a single bubble, assuming that the scattering is purely spherical, can be expressed as the ratio of the scattered power to the incident acoustic intensity, also known as the scattering cross section s s : with r the distance from the bubble where the pressure is measured, |P s (r)| the amplitude of the scattered pressure and |P T | the pressure amplitude of the incident acoustic wave (Hilgenfeldt et al. 1998). This expression also holds at the bubble wall, at r = R 0 . The frequency-dependent scattering cross section can be calculated by using the power spectrum of P s (r) in eqn (14) instead of the total scattered power. Note that through P s , the scattering cross section is now a function of the shell parameters and depends on the frequency and the pressure amplitude of the acoustic wave. For small-amplitude bubble oscillations driven in the linear regime at low acoustic pressure amplitudes, the linearized scattering cross section V s can be expressed as follows (de Jong et al. 1992;Leighton 1994): Equation (15) shows that for a bubble driven at resonance ðv ¼ v 0 Þ and with a total damping of 0.25 (which is typical for UCA bubbles), the scattering cross section is two orders of magnitude larger than its geometric cross-sectional area. Next to scattering, energy of the insonifying ultrasound wave is dissipated owing to viscous damping, thermal damping and shell damping. The total energy loss of the acoustic wave owing to the scattering and damping can therefore be expressed as the ratio of the total energy loss per bubble to the intensity of the acoustic wave: the extinction cross section s e . It is found by multiplying the scattering cross section s s and the ratio of the total damping d tot to the radiation damping d rad (Medwin 1977).
The total amount of scattered and absorbed energy of an acoustic wave propagating through a dilute cloud of homogeneously distributed bubbles at concentration n(R) of bubbles with radius R can now be obtained by summation over the scattering and extinction cross sections of all the individual bubbles present in the acoustic beam (de Jong et al. 1992;Raymond et al. 2014). This linearized approach only holds if the pressure wave scattered by a neighboring bubble at a distance d is much smaller than the average insonifying pressure field: | P s (d)|/|P T | ( 1. At a bubble concentration n, the average distance between bubbles is d = n À1/3 , which, using eqn (14), results in the condition n(R) 2/3 V s ( 1 (Commander and Prosperetti 1989;Goertz et al. 2007). The total attenuation a of an ultrasound wave traveling through a bubble suspension is a result of the total energy dissipation resulting from the presence of bubbles. For linear bubble oscillations, it is given by with d the acoustic path length over which the sound wave interacts with the bubbles. Similarly, again for linear bubble oscillations, the backscatter coefficient h can be calculated as follows (Marsh et al. 1999;Gorce et al. 2000): As discussed before, when the acoustic pressure is increased, the effective shell stiffness decreases. At high driving pressures, the contribution of shell stiffness then becomes negligible and the resonance frequency approaches that of an uncoated bubble. The decrease of the resonance frequency with increasing acoustic pressure is shown in Figure 7a for measured attenuation spectra of a monodisperse microbubble suspension (Segers et al. 2016a). The shift in resonance frequency can be extremely useful in contrast imaging schemes such as power modulation (Eckersley et al. 2005). The transmit frequency can be chosen such that it coincides with the microbubble resonance frequency at the higher acoustic pressures, generating a strong echo only when the transmit power exceeds the threshold for the resonance-frequency downshift.
For non-linearly oscillating bubbles, attenuation can be modeled by solving the full non-linear RP equation to obtain the radius time curves of all bubbles in the size distribution. Then the scattered pressure can be calculated, using eqn (9), and subsequently the scattering cross section s s . Attenuation is then calculated though integration over the size distribution (Hilgenfeldt et al. 2000;Segers et al. 2016a). Note that for contrast bubbles, P s is dependent not only on the ultrasound frequency but also on the acoustic driving pressure. As a consequence, the acoustic bubble response will depend not only on the local pressure at the bubble position in the transmit beam (typically extracted from hydrophone measurements in a calibration tank) but also on the attenuation experienced by the transmit beam on its path to the bubble, which will lower the local acoustic driving pressure. Thus, to accurately model attenuation of a sample of contrast bubbles, the full two-way coupled interac-  tion of the non-linearly oscillating bubbles with the complete ultrasound field must be taken into account (e.g., through forward integration) (Segers et al. 2018b).

ARTICLE IN PRESS
For acoustic characterization of UCAs, the effect of a non-uniform acoustic driving pressure can be minimized by using a narrow bubble sample holder such that it forms a sheet of bubbles in the acoustic characterization beam rather than a 3-D volume (Segers et al. 2018b).

BUBBLE NON-SPHERICITY AND ASYMMETRY
While spherical-bubble models provide a basic fundamental understanding of bubble behavior and the influence of the physical parameters of the surrounding environment, the assumption of spherical symmetry is often challenged. Then, modeling the bubble behavior requires the inclusion of the effects of surface instabilities and neighboring membranes or walls. This breaking of symmetry deeply affects the fundamental behavior of a bubble, namely its resonance frequency and quality factor, and enables, or amplifies, higher-order phenomena such as streaming.
Two main mechanisms can destabilize the spherical shape of a bubble during ultrasound-induced oscillations in the free field: RayleighÀTaylorÀtype instability and parametric instability. RayleighÀTaylor instability can occur at the point of maximum bubble compression while the bubble wall is experiencing maximum acceleration. It can only dominate after repeated, strong collapse events of the bubble or during the after-bounces that follow violent bubble collapse occurring at a low driving frequency (e.g., in sonoluminescence) (Hilgenfeldt et al. 1996). In such a configuration, shape destabilization results in pinch-off, jetting or bubble breakup, and the bubble fate becomes unpredictable. Parametric instabilities, on the other hand, arise during long driving pulses. They grow from an infinitesimally small perturbation on the bubble interface and develop with certain mode preferences depending on the bubble size and driving frequency. Such instabilities are often observed for UCAs (see a typical example in Fig. 8a; Raymond et al. 2016), especially when using long driving pulses (e.g., for therapeutic applications).

Non-spherical oscillations: Modeling
One of the first models for spherical instability was developed by Plesset (1954). In essence, the derivation considers a nearly spherical bubble, with a distortion of small amplitude from the spherical shape that can be written in the form of spherical harmonics: with Y n the spherical harmonic mode of order n and a n its amplitude; R(t) is assumed to obey the spherical RP equation (eqn [1]). Using an analogous derivation of the RP equation, and following a first-order Taylor expansion, we find As before, the overdots represent the time derivatives. This equation is a Hill equation, or a generalized Mathieu equation. The eigenfrequency of the mode directly follows from the small oscillation limit: In this potential flow approach, the velocity field is described by the gradient of the velocity potentialrf, which is a function of space and time. This description, however, does not include any form of damping or viscosity, which results in parametric mode behavior with a threshold of _ R ¼ 0. This model was later picked up by Ceschia and Nabergoj (1978) and Nabergoj and Francescutto (1979), who applied a similar principle by using a Lagrangian approach to find the equation of motion of the bubble from the RP equation and the governing equation of the modes. They assumed purely sinusoidal bubble oscillations; the vorticity is neglected and the viscosity is introduced solely by assuming a vanishing tangential stress. Note that this last assumption is debatable for coated microbubbles, owing to the presence of the shell. This derivation, using a known solution for the resulting Mathieu equation, provides a first estimate of the threshold for the growth of each parametric mode (i. e., the pressure leading to bubble instability). Despite these simplifications, the results have demonstrated good agreement with experimental observations (Fig. 8b; Dollet et al. 2008;Versluis et al. 2010).
Around the same time, Prosperetti (1977b) provided a more rigorous theoretical framework to study surface oscillations of a bubble. Vorticity was accounted for by writing the curl as the sum of a poloidal and a toroidal component. Application to a non-driven bubble provides the same mode frequency and decay rate as given by Nabergoj and Francescutto (1979). Based on the reference work of Prosperetti, subsequent studies (Hilgenfeldt et al. 1996;Wu and Roberts 1999) have lifted some assumptions on the vorticity and instead limited the range of the curl to a small boundary layer whose thickness d BL depends on the bubble size and the mode number: Including this short-range vorticity modifies the second term of eqn (19), where 3 _ R=R must be replaced by the damping term B n (t) given by eqn (13) of Hilgenfeldt et al. (1996).
Two major research directions evolved from this founding work. First, the coupling between bubble motion and surface modes is of particular importance for local drug delivery that involves bubble translation and bubble oscillations. This problem has been mainly investigated in an axisymmetric geometry with viscous effects (Benjamin and Ellis 1990;Shaw 2009) or without (Shaw 2006), demonstrating that shape oscillations can result in self-propulsion and bubble translation. This phenomenon is distinct from radiation forces. Some attention has also been given to the coupling between the surface modes themselves (Mei and Zhou 1991;Feng and Leal 1993;Mao et al. 1995;Harkin et al. 2013).
The second direction consists in modeling the impact of the shell on the existence of surface modes. A few recent studies have investigated the effect of an elastic shell on the parametric oscillations of a coated bubble (Liu et al. 2012;Liu and Wang 2016;Liu et al. 2018). These articles consider linear elasticity theory for the shell, with or without the effect of viscosity, and have determined a correction for the eigenfrequency of the surface modes to account for the presence of the elastic shell. However, polymer-coated ultrasound contrast microbubbles typically fail to display such extreme elastic behavior owing to the buckling or breakup of the fairly rigid shell (Fig. 1h; Marmottant et al. 2011). The influence of a non-linear phospholipid shell, which constitutes the vast majority of contrast microbubble shells, remains mostly unexplored.

Bubbles near a wall
The presence of a boundary near the oscillating bubbles is difficult to avoid, both in vitro and in vivo, and is known to influence the response of the microbubble in a non-trivial way, affecting its oscillation amplitude and non-linear response and shifting its resonance frequency. When the bubble is targeted, and therefore adherent to the wall, bubbleÀwall interactions are inevitable and typically desirable. It is therefore of utmost importance to fundamentally understand the effect of a neighboring wall on the bubble response. In vivo, it is key to enabling many novel applications that require a very fine control of bubble dynamics. In vitro, it is essential to extract meaningful and precise mechanical parameters describing the bubble and its shell as well as its interaction with the boundary.
Physically, the presence of a wall has two distinct effects whose magnitudes depend on its mechanical properties: (i) the bubble interacts with its own acoustic reflection returning from the interface, often described as the image bubble, and (ii) the flow around the bubble deviates from the spherical flow derived from the RP equation, owing to the presence of the wall.
The simplest treatment of the boundary was proposed by Lamb (1932) and consists in adding a pressure term to the RP equation representing the pressure reflected at the interface. Since the bubble is considered to be small compared with the acoustic wavelength, phase lag is disregarded. The reflection decreases the resonance frequency by a factor of ffiffiffiffiffiffiffi ffi 2=3 p . Much later,  investigated the problem from a potential flow perspective. They considered a thin flexible wall and added bubble translation as a result of radiation forces. In this picture, the resonance frequency decreases for a rigid wall and increases for an elastic wall, compared with the free-bubble case. This model was later confirmed experimentally by the same group (Doinikov and Bouakaz, 2015). This result is well aligned with that of Oguz and Prosperetti (1990), who used perturbation theory to derive a model describing the extreme case of a bubble oscillating near a free surface. That work shows that, near a free surface, the resonance frequency should increase by about 40%.
With respect to microbubble oscillations near a wall, of special interest are the non-linear echoes of the microbubbles. These non-linearities are essential for harmonic imaging, including pulse inversion and amplitude modulation, designed to boost the contrast-to-tissue ratio. Suslov et al. (2012) have shown that the dominant harmonic frequencies that are emitted depend on the bubble-to-wall distance. They also suggest an increased chance of period-doubling behavior for bubbles in contact with the wall. The most important effect of a boundary on the bubble behavior, however, lies in the loss of spherical symmetry of the bubble oscillations. Such nonspherical oscillations are systematically observed experimentally (e.g., Fig. 8c), but they are generally overlooked in the modeling work. Recently, Lajoinie et al. (2018) proposed a model for bubbles adherent to or in contact with a boundary, based on a segmentation of the bubble along the inclination angle and with the bubble dynamics solved on a per-segment basis. Despite its extreme simplicity, this model recovers the non-spherical shapes observed experimentally and quantifies the effect of the boundary on the bubble oscillation amplitude and subsequently on streaming (discussed later under "Modeling of higher-order effects"). Much work, however, remains to be done to achieve a rigorous description of non-spherical bubble oscillations near a boundary, especially in relation to resonance shift, and of non-linear emissions.

Bubbles in a tube
As blood-pool agents, contrast microbubbles circulate within the confines of the blood vessel network (i.e., within tubes of various sizes lined with compliant walls). The simplified 1-D case of a bubble of length l oscillating without surface tension or friction in the center of a rigid tube of length L yields (Sun et al. 2009): After linearization around lÀl 0 ð1 þ eÞ for small amplitude e ( 1, Owing to the incompressibility of the liquid, it is clear that significant pulsations are only possible if the tube has an open end. Oguz and Prosperetti (1998) proposed a more complete incompressible disk bubble model that includes the effect of surface tension and the position e of the bubble along the tube axis, with tube radius R tube (Fig. 8d). This model also provides a firstorder estimate for the thermal and viscous damping experienced by such a bubble.
In the same paper, Oguz and Prosperetti (1998) also proposed a model for small bubbles obtained by projecting the velocity potential in the liquid on to the Legendre polynomial base. According to this model, a bubble 10 times smaller than the tube is still strongly affected by the geometry, even when the tube in question is short. Going further, Leighton (2011) has shown, by changing the inertial term in the RP equation, that a long rigid tube 1000 times larger than the bubble can have a significant effect on its behavior.
Let us now consider a bubble trapped in a long tube, replicating a vessel, with a diameter of the same order as the bubble size. If we assume that the compressibility of the surrounding medium does not allow for large-amplitude oscillations (e.g., in water), significant bubble motion can only arise from large compliance of the vessel wall, which allows for bubble motion in three dimensions. The compliance of a blood vessel and its response to stress are complex to model in themselves (Misra and Singh 1983;Humphrey and Na 2002), in particular since blood vessels behave in a non-linear fashion, Qin and Ferrara (2006) have proposed a model that describes the axisymmetric bubble oscillations in a non-linear compliant vessel by adding a pressure jump across the vesselwall location. There, the properties of the blood vessel are fitted to stress measurements on frog mesenteric capillaries. An important conclusion of this model is that, on the one hand, the bubble behavior is modified by the presence of the wall, and on the other hand, the stress experienced by the vessel wall is far greater if a bubble is present within the system during insonation. This model was later amended by Qin and Ferrara (2007) and numerically solved to include the effect of soft tissue around the vessel, providing an extensive study of the impact of wall compliance and its geometric parameters. This line of modeling is of major importance for applications such as the controlled opening of the bloodÀbrain barrier. Owing to the complexity of such modeling, most investigations are performed numerically, but new biomedical applications of bubbles that require more precision and control would also benefit from simplified analytical models.

Bubble jetting
When a bubble undergoes collapse in an asymmetric environment, it forms a jet. The asymmetry can result from the presence of a nearby soft (or gaseous) interface (Krishnan et al. 2017), a rigid boundary, other bubbles in the vicinity (Han et al. 2015;Kwan et al. 2016) or even a non-linear acoustic wave (Ohl and Ikink 2003). Jet formation can have severe consequences; the high jet velocity is crucial in bubble-induced erosion and in bubble cleaning applications (Fernandez Rivas et al. 2012;Verhaagen and Fernandez Rivas 2016). Jetting can lead to the breakup of bubbles into numerous daughter bubbles, the destruction of the agent, or the deposition of the shell material onto the rigid boundary (Roovers et al. 2019). Jetting becomes more prominent as bubble size increases (owing to reduced stabilization by surface tension), as frequency decreases and as pressure increases (owing to increased inertia).
In terms of fluid dynamics, bubble jetting can be described using fairly strong assumptions: the flow is considered inviscid and irrotational, leading to a potential flow description. This approach was first presented by Plesset and Chapman (1971), who predicted jet velocities on the order of 100 m/s. The shapes that were computed based on this approach were later experimentally verified by Lauterborn and Bolle 1975 (Fig. 9a) with remarkable agreement. Going beyond this initial description, Blake et al. (1997) implemented the presence of the boundary using Green's function for the velocity potential to account for the source (the bubble) and its image. This Green's function enforces a no-flow condition across the boundary.
To better account for the large jet velocities involved and the presence of an acoustic wave, Wang and Blake (2010) extended these formulations to slightly compressible fluids, demonstrating a measurable effect of compressibility on bubble jetting and in particular on the orientation of the jet as compared with the propagation direction of the ultrasound wave. Next to bubble size, ultrasound frequency and pressure, jetting is also very sensitive to the initial distance between the bubble and the rigid boundary (Fig. 9bÀd) (Prentice et al. 2005;Wang et al. 2015). In particular, the problem was treated by Brujan (2017) in a weakly compressible theory, showing that jet velocity rapidly increases as the initial bubbleÀboundary distance increases, while the total amount of kinetic energy in the jet overall decreases. Biological tissues present a wide range of mechanical properties but are nearly all flexible in the sense that they can be deformed during bubble oscillations. This is true for cell membranes as well as for blood-vessel walls. It is therefore of prime importance to describe the effect of soft surfaces on jetting, which is expected to be a major cause of damage in vivo. It has been shown that the jet direction can be reversed when the substrate material becomes very soft, namely when it becomes easier to pull the elastic material than to overcome the water inertia (Gibson and Blake 1982;Brujan et al. 2001;Li et al. 2018). There have, however, been only a few studies on the transition regimes and the effect of excised biological tissues (Chen et al. 2011).

Modeling bubbles in tissue and soft matter
Given the complex nature of biological tissues and the interest in novel extravasating agents and bubble precursors, it is important to study the behavior of bubbles in viscoelastic materials. In particular, the bubble resonance frequency is known to change significantly with the elastic modulus of the enclosing medium. A number of mechanical models have been developed-initially based on KelvinÀVoigt models (Macosko 1994) and later on Maxwell models, both linear (Allen and Roy 2000a) and non-linear (Allen and Roy 2000b; Naude and Mendez 2008)-demonstrating the change in resonance frequency as a result of the elastic modulus of the bubble surroundings. De  investigated the dynamics of large bubbles (>50 mm) immersed in a yield-stress fluid using theory and numerical simulations. They quantified the viscoelastic and viscoplastic response of the material using the model developed by Saramito (2009). In another notable example, Yang and Church (2005) investigated the increased threshold for subharmonic emissions in these materials. Linear and non-linear viscoelastic models are more extensively reviewed by Dollet et al. (2019).

Streaming and shedding
Acoustic streaming is a phenomenon whereby an oscillating flow generates a non-linear, second-order steady flow (Fig. 10a, 10b). In particular, it occurs in the presence of walls, along which viscous effects are localized in a boundary layer of typical thickness with v the angular frequency. Rayleigh (1884) was the first to describe acoustic streaming theoretically. He showed that an oscillatory flow of velocity U(x)cos vt along a wall (x being a coordinate along the wall) causes a steady-state streaming flow of velocity u s such that u s ðxÞ ¼Àð3=2ÞUðxÞU 0 ðxÞ=v along the wall (e.g., Lighthill 1978); to be rigorous, this applies only if the streaming Reynolds number u s r L '/m L is much lower than 1 ( ' being the typical scale over which streaming takes place).
Acoustic streaming can be generated by an oscillating bubble close to a wall (Marmottant and Hilgenfeldt 2003) (Fig. 10c). It is thought to be a key transport mechanism for drug release from drug-loaded UCAs (Luan et al. 2014), and one of the main mechanisms of sonoporation through the induced shear stress on a cell layer. To theoretically assess the latter effect in the specific context of UCAs, Doinikov and Bouakaz (2010) computed the shear stress associated with streaming. Generally speaking, it is of order T xz % m L u s /d BL , z being the coordinate away from the wall. Hence, T xz % ffiffiffiffiffi r L p m L =vUðxÞU 0 ðxÞ. Doinikov and Bouakaz (2010) also computed the velocity U related to a microbubble located at a distance d from the wall and undergoing volumetric oscillations of radius R(t) forced by an acoustic field, using a modified RP equation with the model by de Jong et al. (1994) for the shell term, and obtained the corresponding shear stress. They showed that it displays a maximum at a finite radial distance along the wall, and estimated its potential implications for sonoporation. Chomas et al. (2001) presented one of the first systematic studies of UCA destruction under insonation. The most violent destruction mechanism is fragmentation into smaller bubbles at high amplitude, but at lower amplitude they evidenced a progressive deflation of the UCA, which they argue is owing to acoustically driven diffusion. Actually, when lipid monolayers get sufficiently compressed, they can shed some material in the surrounding liquid. This mechanism has been recognized as important in the context of the slow dissolution of lipid-coated bubbles (Borden and Longo 2002) (Fig. 1i) and has been proposed as key for the understanding of UCA deflation under prolonged insonation. Using fluorescence imaging, Borden et al. (2005) evidenced ultrasound-induced shedding and quantified the time evolution of bubble deflation toward an equilibrium (around 2 mm in radius) (Fig. 10f). They showed different mechanisms of shedding, depending on lipid composition, from submicron particle expulsion to budding, then to buildup and expulsion of larger lipid aggregates as the shell becomes more cohesive. Similar observations have been made by Kooiman et al. (2017) using ultrahigh-speed imaging. Kwan and Borden (2012) argue that the microscopic mechanisms initiating shedding are the local collapse of the lipid monolayer and the subsequent aggregation of new folds. Guidi et al. (2010) subjected UCAs to repeated low-amplitude ultrasonic pulses at a frequency of 2 MHz, and consistently observed deflation for bubbles with an initial radius below 3 mm toward a stable radius of 1.4 mm, but no deflation for larger bubbles.
The time evolution of the deflation process is qualitatively similar to the one reported by Borden et al. (2005), with the fastest deflation rate near resonance, where UCAs are more prone to dissolve or shed their (c) Simulation of the acoustic streaming generated by an oscillating bubble of radius a (the half-circle) in contact with a plane wall located at z = 0. Only half of the space is sketched, owing to axisymmetry about the x = 0 axis. The solid curves are the streamlines, and the direction of velocity is indicated with the arrows; in particular, the streaming flow goes away from the wall above the bubble (reproduced with permission from Marmottant et al. 2006). (d) Lipid shedding from a targeted bubble and transport by acoustic streaming. (Top) Schematics of the shedding process of the released material (red) with the calculated streamlines (solid black lines); (bottom) experimental recording of the release of the fluorescent material from a targeted oscillating microbubble driven at a pressure of 331 kPa, showing transport over a distance several times larger than the bubble size (reproduced with permission from Lajoinie et al. 2018). (E) Streak imaging of two ultrasound contrast agents of equilibrium radius 2.1 mm (top) and 2.0 mm (bottom), showing the displacement in the direction of propagation of the acoustic wave. The acoustic forcing represents 20 cycles at amplitude 380 kPa and frequency 2.25 MHz. These snapshots are overlaid with simulated curves of the time evolution of the bubble radius (reproduced with permission from Dayton et al. 2002). (F) Evolution of the radius of a lipid-coated microbubble subjected to repeated acoustic pulses (frequency 2.25 MHz, peak negative pressure 400 kPa, pulse repetition frequency 0.5 Hz; data kindly provided by Mark Borden).
coating owing to large volumetric oscillations. Viti et al. (2012) focused on megahertz bubble oscillations using ultrahigh-speed imaging while the bubbles underwent (slower) deflation. They showed that deflation occurs in discrete steps, and that the radius time oscillations usually evolve from being rather symmetric just after a deflation event to compression-only behavior after some time in a non-deflating state, thereby suggesting that the lipid shell reforms after each shedding event.
The idea that a sufficiently compressed lipid monolayer collapses and expels excess lipids is the basis of the surfactant-shedding model of O'Brien et al. (2011), based on phenomenological dependences of molecular diffusivity and surface viscosity on surface concentration of lipids. This model reproduced the deflation results of Borden et al. (2005). To clarify lipid shedding dynamics under insonation, Luan et al. (2014) used high-speed fluorescence imaging at 150,000 frames/s. They showed that shedding may occur already within a few ultrasound pulses, in contrast to the picture of a (relatively slow) build up of expelled material. Lajoinie et al. (2018) refined this setup by coupling it with ultrahigh-speed imaging. They showed that non-spherical oscillations of a lipid-coated bubble in contact with a wall lead to local oversaturation and shedding, and that the subsequent transport of the shed lipid is mediated by acoustic streaming (Fig. 10d).

Radiation forces
Any body of volume V in a space-dependent pressure field p experiences a force equal toÀV r p, provided that the size of the body is much smaller than the acoustic wavelength. If the body is a bubble in an acoustic field, this can lead to a strong time-averaged force, owing to the large gas compressibility. This force is termed acoustic radiation force or, in the bubble context, Bjerknes force (Bjerknes 1906;Crum 1975;Leighton 1994). Two radiation forces are usually distinguished: primary forces, when the bubble responds directly to the external applied pressure field (typically from a transducer), and secondary forces, when a given bubble reacts to the pressure field re-emitted by a second neighboring bubble or, in the case of a neighboring wall, by the image bubble. These forces are maximal at the resonance frequency of the bubble. For traveling waves, the primary radiation force always acts in the direction of ultrasound wave propagation. Dayton et al. (2002) characterized the influence of primary radiation forces on UCAs. They showed that bubbles with a radius of 1.5 mm can be displaced by as much as 5 mm by a 20-cycle pulse with a frequency of 2.25 MHz at an amplitude of 180 kPa (Fig. 10E). However, assessing radiation forces quantitatively turns out to be difficult because the force balance involves hydrodynamic forces, including added mass, drag and history force. Dayton et al. (1999), Zhao et al. (2004), Rychak et al. (2005), Lum et al. (2006) and Frinking et al. (2012) have all used primary radiation forces as a strategy to achieve targeted drug delivery. Segers and Versluis (2014) proposed a novel use of primary radiation forces: using the dependency of this force on bubble resonance, they were able to deflect and isolate resonant bubbles from the main flow of a polydisperse suspension of microbubbles and collect them using a dedicated microfluidic device. A simple force-balance model for the acoustic sorting showed very good agreement with the experiment. Thus, they were able to narrow down the size distribution and thereby to improve the sensitivity of these UCAs by more than 10 dB (Segers et al. 2018b).
Secondary radiation forces lead to attraction of bubbles of similar size and repulsion of bubbles of different size if the driving frequency is in between the two resonance frequencies (Leighton 1994). The effect of secondary radiation forces on UCAs has been studied by Garbin et al. (2011) and Kokhuis et al. (2013) for targeted microbubbles. They confirmed that UCAs of similar sizes attract. For bubbles isolated from a neighboring wall using optical tweezers, it has been shown that using the same kind of force balance as Dayton et al. (2002), one must include nontrivial added mass and history forces stemming from delayed viscous effects caused by the backand-forth motion of the bubble in the vorticity generated in its own wake (Garbin et al. 2009).

MODELS FOR NOVEL AGENTS
The creation of new classes of novel agents also triggers the need for novel models. First, bubbles can be loaded with functional nanoparticles as a method to enhance drug delivery. Such self-assembled systems include microbubbles that are conjugated with liposomes (Kheirolomoom et al. 2007;Lentacker et al. 2010;Luan et al. 2012;Peyman et al. 2012;Yan et al. 2013), echogenic liposomes (Paul et al. 2012;Raymond et al. 2015), colloids (Poulichet and Garbin 2015;Huerre et al. 2018) and polymeric nanoparticles (Borden and Rege 2012;Burke et al. 2012;Mørch et al. 2015;Z. Chen et al. 2019;Jamburidze et al. 2019) (Fig. 11a, 11b). Encapsulation of plasmonic, photoacoustic or magnetic nanoparticles allows for multimodality contrast imaging, with the aim of combining imaging techniques and thereby compensating for the inherent weaknesses of one modality by exploiting the strengths of another. Similarly, bubbles can be added to other constructs to add or enhance a specific functionality. For example, acoustic cluster therapy (Healey et al. 2016) utilizes an agent comprised of clusters of microdroplets conjugated with Sonazoid bubbles, where the microbubbles lower the activation threshold for the formation of gas bubbles, or emboli, from the microdroplets upon ultrasound exposure. This confines systemically administered chemotherapy drugs within the tumor to enable diffusion. To overcome the high washout rates of cardiac stem cells, Naaijkens et al. (2013) decorated adipose-derived stem cells with functionalized microbubbles (Fig. 11c), effectively creating echogenic complexes, dubbed StemBells, that are then susceptible to acoustic radiation force. The behavior of these complexes has been modeled by Kokhuis et al. (2017) based on a modified RP equation, showing both theoretically and experimentally how radiation forces can be used to push the StemBells out of the mainstream flow and toward the vessel wall. This approach leads to a significant increase in stem cell delivery (Kokhuis et al. 2015;Woudstra et al. 2016) to the damaged area of the heart. Finally, nano-sized agents, including nanobubbles, nanodroplets and gas vesicles, are a class of their own. Models describing the behavior of these novel contrast agents are highly needed to understand and further exploit their linear and non-linear contrast features. A more extensive discussion on these novel agents is provided in another review of this issue on future directions .

Modeling bubbles for multimodality imaging
Since microbubbles are unparalleled agents for ultrasound and display remarkable (non-linear) properties, there is great interest in driving them with other sources than ultrasound (Owen et al. 2018). In particular, researchers have combined microbubbles with plasmonic nanoparticles (Wilson et al. 2012), magnetic nanoparticles (Stride et al. 2009;Vlaskou et al. 2010) or an additional liquid or polymeric layer that contains a dye (Lajoinie et al. , 2017aVisscher et al. 2019) (Fig. 1e, 1f). These constructs aim at making the bubbles light-sensitive so that they can generate contrast in both ultrasound and photoacoustic imaging. Heat transfer in an acoustically driven bubble was investigated decades ago, in particular by Prosperetti (1977a), which led to the well-known thermal damping in the modified RP equation (discussed earlier in this article). The thermal exchanges in a bubble driven into oscillations using light are, however, different from the classical picture, and one can no longer assume polytropic behavior of the gas. Heat diffusion must be resolved in space and time, together with the momentum, to describe the bubble oscillations. Initial models have been proposed for bubbles coated with a dye-doped oil layer (Lajoinie et al. 2017b) that qualitatively explain the experimental observations. Further modeling is required to accurately describe and understand both linear and non-linear bubble oscillations for these systems. In particular, this modeling must include the partial pressures of the surrounding liquids for multilayered bubbles and a more detailed description of the optical and mechanical properties of the bubble shells once they are coated with plasmonic nanoparticles. Modeling bubbles formed from droplet vaporization A sufficiently intense burst of ultrasound can vaporize otherwise inert droplets, forming echogenic microbubbles, a process termed acoustic droplet vaporization (ADV) (Kripfgans et al. 2000;Rapoport et al. 2007;Williams et al. 2013;Sheeran et al. 2014) (Fig. 11d). The freshly formed bubbles are stabilized over timescales of seconds to minutes by the same coating material initially used to encapsulate their liquid precursors . There is great interest in the use of these agents for both diagnosis and therapy, for example, as boosters for high-intensity focused ultrasound therapy (Zhang and Porter 2010;Phillips et al. 2013), embolotherapy , dissolved oxygen scavengers (Radhakrishnan et al. 2016) and vehicles for drug delivery (Rapoport 2012).
Contrast microbubble size is a disadvantage for a number of applications because the bubbles are confined to the blood vasculature and cannot easily extravasate into the interstitium. Liquid perfluorocarbon nanodroplets, however, have been studied as novel extravascular contrast agents for ultrasound. With a typical size of a few hundred nanometers, these droplets have the ability to extravasate selectively into porous tumor tissue-presumably through the enhanced permeability and retention effect (Maeda et al. 2000)-where they can then be vaporized for imaging or drug delivery purposes.
The mechanisms underlying droplet nucleation have received considerable interest (Kripfgans et al. 2000;Sheeran et al. 2011;Shpak et al. 2014). This paragraph now focuses on bubble stabilization upon vaporization for the purpose of generating stable contrast on clinically relevant timescales (hundreds of milliseconds). Reznik et al. (2013) studied the pathways after initial vaporization, demonstrating the complex interaction of the vapor bubble with its surrounding, including coalescence, gas diffusion and rectified heat transfer (Shpak et al. 2013a;Doinikov et al. 2014). Shpak et al. (2013b) simulated the role of gas in ultrasonically driven vaporbubble growth and underlined the fundamental role of gas diffusion during acoustic forcing. They showed that an amount of air as low as 10% of the saturation concentration in the perfluorocarbon is enough to allow the bubble to survive the first collapse. Mercado-Shekhar et al.
(2019) exploited this rapid diffusion property for an ultrasound-mediated gas-scavenging technology for dissolved oxygen. There, the oxygen transfer is driven by the gas concentration gradients resulting from droplet vaporization. Outside the context of ADV, Sarkar et al. (2009), Kwan and Borden (2010) and Segers et al. (2016b) have also modeled the important effects of gas exchange for bubble stability, including the effects of lipid formulation (Sz ıjj art o et al. 2012;Garg et al. 2013;van Rooij et al. 2015).

Modeling nanobubbles and nanovesicles
In recent years, nanobubbles have attracted considerable attention. Their existence, however, seems to be paradoxical. A simple classical estimate based on the Laplace pressure inside such bubbles suggests that they should dissolve in microseconds. However, once stabilized by surface heterogeneities, pinning of the three-phase contact line can occur and these surface nanobubbles can be stable for days as the gas influx from the gas-oversaturated aqueous surrounding and the outflux owing to Laplace pressure balance (Lohse and Zhang 2015). One such system developed in the context of cavitation therapy is a novel ultrasound-responsive single-cavity polymeric nanoparticle, or nanocup (Fig. 11e), capable of trapping and stabilizing gas against dissolution in the bloodstream (Kwan et al. 2015). Ultrahigh-speed imaging of the nucleated bubble shows an initial growth phase followed by detachment from the nanocup and subsequent inertial collapse, confirming the underpinning hypothesis. These experiments were quantitatively predicted using a modified RP model (Kwan et al. 2016). Direct observation of bubble nucleation from nanopits (Borkent et al. 2009), which typically reside in mesoporous silica nanoparticle agents, is out of reach owing to the small spatial and temporal scales involved. Thus, for now we must rely on indirect measurements by scanning the vast parameter space of ambient pressure, surfactant concentration, wetting properties and gas saturation, going hand in hand with theory and numerical simulations, including finite-difference simulations and molecular-dynamics simulations.
Finally, naturally occuring nano-sized gas-filled vesicles can be bioengineered to form " a la carte" biomolecules to non-invasively image and control biological functions using ultrasound ). These gas vesicles (Fig. 11f) form a unique class of protein nanostructures that are detectable at subnanomolar concentrations and whose physical properties allow them to serve as highly sensitive imaging agents for ultrasound (Shapiro et al. 2014). Such gas vesicles also produce a non-linear response, enabling contrast imaging with these agents using amplitude modulation . While direct observation of their shell dynamics under ultrasound exposure is so far out of reach, a finite-element model can be used to predict shell buckling behavior for these vesicles, leading to non-linear ultrasound backscatter at both their fundamental and second harmonic frequencies of 10À20 MHz when exposed to ultrasound above a threshold pressure of 300 kPa.