## Abstract

Ultrasound is extensively used in medical imaging, being safe and inexpensive and operating in real time. Its scope of applications has been widely broadened by the use of ultrasound contrast agents (UCAs) in the form of microscopic bubbles coated by a biocompatible shell. Their increased use has motivated a large amount of research to understand and characterize their physical properties as well as their interaction with the ultrasound field and their surrounding environment. Here we review the theoretical models that have been proposed to study and predict the behavior of UCAs. We begin with a brief introduction on the development of UCAs. We then present the basics of free-gas-bubble dynamics upon which UCA modeling is based. We review extensively the linear and non-linear models for shell elasticity and viscosity and present models for non-spherical and asymmetric bubble oscillations, especially in the presence of surrounding walls or tissue. Then, higher-order effects such as microstreaming, shedding and acoustic radiation forces are considered. We conclude this review with promising directions for the modeling and development of novel agents.

## Key Words

## Introduction

While in the 17

^{th}century Robert Hooke had already predicted that we would be able to image the human body with sound (), 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

; - Mulvagh S.L.
- Rakowski H.
- Vannan M.A.
- Abdelmoneim S.S.
- Becher H.
- Bierig S.M.
- Burns P.N.
- Castello R.
- Coon P.D.
- Hagen M.E.
- Jollis J.G.
- Kimball T.R.
- Kitzman D.W.
- Kronzon I.
- Labovitz A.J.
- Lang R.M.
- Mathew J.
- Moir W.S.
- Nagueh S.F.
- Pearlman A.S.
- Perez J.E.
- Porter T.R.
- Rosenbloom J.
- Strachan G.M.
- Thanigaraj S.
- Wei K.
- Woo A.
- Yu E.H.
- Zoghbi W.A.

American Society of Echocardiography consensus statement on the clinical applications of ultrasonic contrast agents in echocardiography.

*J Am Soc Echocardiogr.*2008; 21: 1179-1201

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 air-filled 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

; - Skyba D.M.
- Camarano G.
- Goodman N.C.
- Price R.J.
- Skalak T.C.
- Kaul S.

Hemodynamic characteristics, myocardial kinetics and microvascular rheology of FS-069, a second-generation echocardiographic contrast agent capable of producing myocardial opacification from a venous injection.

*J Am Coll Cardiol.*1996; 28: 1292-1300

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 (- Cohen J.L.
- Cheirif J.
- Segar D.S.
- Gillam L.D.
- Gottdiener J.S.
- Hausnerova E.
- Bruns D.E.

Improved left ventricular endocardial border delineation and opacification with OPTISON (FS069), a new echocardiographic contrast agent: Results of a phase III multicenter trial.

*J Am Coll Cardiol.*1998; 32: 746-752

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 μm, 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

).European Medicines Agency. SonoVue: EPAR—Product information. Available at: https://www.ema.europa.eu/en/documents/product-information/sonovue-epar-product-information_en.pdf. Accessed January 22, 2020.

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 bubbles 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 (- Pereno V.
- Aron M.
- Vince O.
- Mannaris C.
- Seth A.
- de Saint Victor M.
- Lajoinie G.
- Versluis M.
- Coussios C.
- Carugo D.
- Stride E.

Layered acoustofluidic resonators for the simultaneous optical and acoustic characterisation of cavitation dynamics, microstreaming, and biological effects.

*Biomicrofluidics.*2018; 12034109

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 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 lead to surprising contrast response characteristics, including thresholding behavior (

Emmer et al., 2007

), compression-only behavior (de Jong et al., 2007

; Sijl et al., 2011

) and subharmonics (Sijl et al., 2010

). 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 toThornycroft 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 19^{th}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 this was originally given by

Taylor and Davies, 1963

, although it is generally attributed to Poritsky, 1951

. The result was the Rayleigh-Plesset or Rayleigh–Plesset–Noltingk–Neppiras–Poritsky equation (Lauterborn, 1976

):${\rho}_{\mathrm{L}}\left(R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}\right)={P}_{\mathrm{B}}-{P}_{\infty}.$

(1)

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}; $\dot{R}$ and $\ddot{R}$ are the velocity and acceleration, respectively, of the bubble wall;

*ρ*

_{L}is the density of the surrounding fluid;

*P*

_{B}is the pressure in the liquid at the bubble wall; and

*P*

_{∞}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.

Table 1List of symbols

Symbol | Meaning |
---|---|

A | Bubble surface area |

A_{0} | Bubble surface area at rest |

an | Amplitude of surface mode number n |

α | Total attenuation |

c_{L} | Speed of sound in the liquid |

d | Acoustic path length |

dA | Bubble surface-area dilatation |

δ_{BL} | Boundary-layer thickness |

δ_{tot} | Total damping |

δ_{rad} | Radiative damping |

δ_{vis} | Liquid viscous damping |

δ_{shell} | Shell viscous damping |

δ_{th} | Thermal damping |

e | Bubble location in the tube |

ϵ | Small oscillation amplitude (dimensionless) |

ϕ | Azimuthal angle of the spherical harmonics |

η | Backscatter coefficient |

f | Ultrasound frequency |

f_{T} | Transmit frequency |

f_{res} | Resonance frequency |

f_{0} | Bubble eigenfrequency |

κ | Polytropic exponent |

κ_{s} | Shell viscosity |

l | One-dimensional bubble length |

l_{0} | One-dimensional bubble length at rest |

L | One-dimensional tube length |

ℓ | Typical streaming length scale |

μ_{L} | Liquid viscosity |

n | Spherical harmonic mode number |

n(R) | Bubble concentration of bubbles with radius R |

P_{0} | Hydrostatic pressure |

P_{∞} | Pressure in the liquid far from the bubble |

P_{A} | Driving pressure |

P_{B} | Pressure in the liquid at the bubble wall |

p_{v} | Vapor pressure |

P_{s} | Scattered pressure |

P_{T} | Transmit pressure |

P_{vis} | Viscous pressure term |

P_{elas} | Elastic pressure term |

R | Bubble radius |

$\dot{R}$ | Bubble wall velocity |

$\ddot{R}$ | Bubble wall acceleration |

R_{0} | Bubble radius at rest |

R_{b} | Bubble buckling radius |

R_{r} | Bubble rupture radius |

R_{tube} | Tube radius |

ρ_{L} | Liquid density |

σ | Surface tension |

σ(R) | Size-dependent surface tension |

σ(R_{0}) | Surface tension at rest |

σ_{s} | Scattering cross section |

t | Time |

θ | Polar angle of the spherical harmonics |

τxz | Shear stress along the wall |

U | Oscillatory flow velocity |

u_{s} | Streaming velocity |

V | Bubble volume |

ω | Angular ultrasound frequency |

ω_{0} | Angular bubble eigenfrequency |

ω_{0}n | Angular bubble eigenfrequency of surface mode n |

Ω_{s} | Linearized scattering cross section |

x | Coordinate along the wall |

χ | Shell elasticity |

Yn | Spherical harmonic of mode number n |

z | Coordinate away from the wall |

When the gas in the bubble is an ideal gas, from the pressure balance at the bubble interface

following the ideal gas law, with

*P*_{B}can be expressed as${P}_{\mathrm{B}}=\left({P}_{0}+\frac{2\sigma}{{R}_{0}}-{p}_{\mathrm{v}}\left({R}_{0}\right)\right){\left(\frac{{R}_{0}}{R}\right)}^{3\kappa}-\frac{2\sigma}{R}-\frac{4{\mu}_{\mathrm{L}}\dot{R}}{R}+{p}_{\mathrm{v}},$

(2)

following the ideal gas law, with

*σ*the surface tension at the gas–liquid interface,*κ*the polytropic constant and*μ*_{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 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 $\left(|\dot{R}|/{c}_{\mathrm{L}}\sim 1\right)$. 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 liquid and the local speed of sound at the bubble wall.In 1980, building upon the earlier models,

with

For an exhaustive review of the derivation of the different RP-type equations, see

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:$\left(1-\frac{\dot{R}}{{c}_{\mathrm{L}}}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{\dot{R}}{3{c}_{\mathrm{L}}}\right){\dot{R}}^{2}=\frac{1}{{\rho}_{\mathrm{L}}}\left(1+\frac{\dot{R}}{{c}_{\mathrm{L}}}+\frac{R}{{c}_{\mathrm{L}}}\frac{d}{dt}\right)\left({P}_{\mathrm{B}}-{P}_{\infty}\right),$

(3)

with

*P*_{B}and*P*_{∞}as defined before. This expression was subsequently refined byProsperetti 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 $|\dot{R}|/{c}_{\mathrm{L}}\sim 1$. Therefore, a popular form in the context of sonoluminescence and UCA modeling is given by (Löfstedt et al., 1995

; Barber et al., 1997

; Brenner et al., 2002

; Marmottant et al., 2005

; Paul et al., 2010

):${\rho}_{\mathrm{L}}\left(R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}\right)=\left({P}_{0}+\frac{2\sigma}{{R}_{0}}\right){\left(\frac{{R}_{0}}{R}\right)}^{3\kappa}\left(1-\frac{3\kappa \dot{R}}{{c}_{\mathrm{L}}}\right)-\frac{2\sigma}{R}-\frac{4{\mu}_{\mathrm{L}}\dot{R}}{R}-{P}_{0}-{P}_{\mathrm{A}}\left(t\right).$

(4)

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 *hard-shelled 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 semi-empirical 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 Doinikov and Bouakaz, 2011

.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 lipid-coated 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 (de Jong et al., 2007

; 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., 1992

, de Jong 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 Doinikov and Bouakaz, 2011

. For a more extensive review, see Faez et al., 2013

.The viscous pressure term

with

*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

):${P}_{\text{vis}}=4{\kappa}_{\mathrm{s}}\frac{\dot{R}}{{R}^{2}},$

(5)

with

*κ*_{s}the surface dilatational viscosity of the lipid monolayer shell. Typically, it is assumed that*κ*_{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

Assuming a constant and surface-area–independent shell elasticity, eqn (6) can be integrated to find the effective interfacial tension: $\sigma \left(R\right)=\chi \left({R}^{2}/{R}_{0}^{2}-1\right)$, with

Equation (7) is the most common form of the elastic pressure term (

where

*P*_{elas}accounts for the pressure increase inside the bubble owing to the stiffness of the shell. Shell elasticity*χ*can be described to result from interfacial tension gradients*dσ*induced by (acoustically-driven) bubble surface-area variations*dA*(Glazman, 1983

; ; Marmottant et al., 2005

; Sarkar et al., 2005

):$\chi =A\frac{d\sigma}{dA}.$

(6)

Assuming a constant and surface-area–independent shell elasticity, eqn (6) can be integrated to find the effective interfacial tension: $\sigma \left(R\right)=\chi \left({R}^{2}/{R}_{0}^{2}-1\right)$, with

*R*_{0}the equilibrium bubble radius. With the effective interfacial tension, the elastic pressure term reads${P}_{\text{elas}}=\frac{2\sigma \left(R\right)}{R}=\frac{2\chi}{R}\left(\frac{{R}^{2}}{{R}_{0}^{2}}-1\right).$

(7)

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:${\rho}_{\mathrm{L}}\left(R\ddot{R}+\frac{3}{2}{\dot{R}}^{2}\right)=\left({P}_{0}+\frac{2\sigma \left({R}_{0}\right)}{{R}_{0}}\right){\left(\frac{{R}_{0}}{R}\right)}^{3\kappa}\left(1-\frac{3\kappa \dot{R}}{{c}_{\mathrm{L}}}\right)-{P}_{0}-{P}_{\mathrm{A}}\left(t\right)-\frac{4{\mu}_{\mathrm{L}}\dot{R}}{R}-\frac{2\chi}{R}\left(\frac{{R}^{2}}{{R}_{0}^{2}}-1\right)-\frac{4{\kappa}_{\mathrm{s}}\dot{R}}{{R}^{2}},$

(8)

where

*σ*(*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

Thus, the scattered pressure

*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

):${P}_{\mathrm{s}}=\frac{{\rho}_{\mathrm{L}}}{4\pi r}\ddot{V}=\frac{{\rho}_{\mathrm{L}}}{r}\left({R}^{2}\ddot{R}+2R{\dot{R}}^{2}\right).$

(9)

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 μm. Three driving pressures were used here—5, 50 and 500 kPa—corresponding 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 $\dot{R}$ and $\ddot{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 μm 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 (

with $\dot{\epsilon}$ and $\ddot{\epsilon}$ the first- and second-order derivatives of ε with respect to time,

In these equations,

*i.e.*, $R={R}_{0}\left(1+\epsilon \left(t\right)\right)$, with ε (*t*) ≪ 1), and by neglecting higher-order terms to yield the classic equation of a driven harmonic oscillator:$\ddot{\epsilon}+{\delta}_{\text{tot}}{\omega}_{0}\dot{\epsilon}+{\omega}_{0}^{2}x=\frac{{P}_{\mathrm{A}}}{\rho {R}_{0}}sin\omega t,$

(10)

with $\dot{\epsilon}$ and $\ddot{\epsilon}$ the first- and second-order derivatives of ε with respect to time,

*δ*_{tot}the total damping of the system,*ω*the angular ultrasound frequency and*ω*_{0}the angular eigenfrequency of the oscillator:${\omega}_{0}=2\pi {f}_{0}=\frac{1}{{R}_{0}}\sqrt{\frac{1}{{\rho}_{\mathrm{L}}}\left(3\kappa {P}_{0}+\left(3\kappa -1\right)\frac{2\sigma \left({R}_{0}\right)}{{R}_{0}}+\frac{4\chi}{{R}_{0}}\right).}$

(11)

In these equations,

*σ*(*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*χ*of the coating increases the eigenfrequency of the coated bubble compared with the free gas bubble.The total damping

This expression for radiation damping is in approximate agreement with the exact form ${\delta}_{\text{rad}}=\omega {R}_{0}/{c}_{\mathrm{L}}$ (). A microbubble also experiences thermal damping

*δ*_{tot}is a summation of the damping contributions that arise from reradiation of sound, viscous dissipation in the surrounding liquid and shell viscous dissipation: ${\delta}_{\text{tot}}={\delta}_{\text{rad}}+{\delta}_{\text{vis}}+{\delta}_{\text{shell}}$. The linearization of eqn (8) results in the following expressions for the damping terms at*ω*=*ω*_{0}:${\delta}_{rad}=\frac{3\kappa {P}_{0}}{{\rho}_{L}{\omega}_{0}{c}_{L}{R}_{0}},{\delta}_{vis}=\frac{4{\mu}_{L}}{{\rho}_{L}{\omega}_{0}{R}_{0}^{2}},{\delta}_{shell}=\frac{4{\kappa}_{s}}{{\rho}_{L}{\omega}_{0}{R}_{0}^{3}}.$

(12)

This expression for radiation damping is in approximate agreement with the exact form ${\delta}_{\text{rad}}=\omega {R}_{0}/{c}_{\mathrm{L}}$ (). A microbubble also experiences thermal damping

*δ*_{th}owing to heat diffusion (Prosperetti, 1976

; ). 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 (van der Meer et al., 2007

; 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 (*χ*= 0.5 N/m). Shell damping was calculated using a fixed shell viscosity*κ*_{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 μm. The total damping*δ*_{tot}of the system results in a slightly reduced microbubble resonance frequency*f*_{res}as compared with its eigenfrequency*f*_{0}: ${f}_{\text{res}}={f}_{0}\sqrt{1-{\delta}_{\text{tot}}^{2}/2}$. 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 μm) 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*κ*_{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).#### 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*χ*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 (de Jong et al., 2007

) and the resonance-frequency shift with increasing acoustic pressure amplitudes for lipid-coated bubbles (Overvelde et al., 2010

).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}:$\sigma \left(R\right)=\{\begin{array}{cc}0\hfill & \text{if}\phantom{\rule{1em}{0ex}}R\le {R}_{\mathrm{b}}\\ \chi \left(\frac{{R}^{2}}{{R}_{\mathrm{b}}^{2}}-1\right)\hfill & \text{if}\phantom{\rule{1em}{0ex}}{R}_{\mathrm{b}}\le R\le {R}_{\mathrm{r}}\\ {\sigma}_{\text{water}}\hfill & \text{if}\phantom{\rule{1em}{0ex}}R\ge {R}_{\mathrm{r}}.\end{array}$

(13)

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 $\left(\chi =Ad\sigma /dA=0\right)$; 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 [Borden et al., 2006

]) 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 *σ*_{water}of the surrounding medium is reached, again leading to a tensionless state $\left(\chi =Ad\sigma /dA=0\right)$. 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*χ*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

*σ*(*R*_{0}) of the bubble. At*σ*(*R*_{0}) <*σ*_{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*σ*(*R*_{0}) >*σ*_{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 $\sigma \left({R}_{0}\right)={\sigma}_{\text{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 (Emmer et al., 2007

). A detailed explanation of thresholding behavior is given by Overvelde et al., 2010

.Next to compression-only and driving-pressure–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., 2010

, 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}*σ*/*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

*σ*(*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 $\chi \left(A/{A}_{0}\right)$ curve, Figure 6f, has been integrated to find the $\sigma \left(A/{A}_{0}\right)$ 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 $\chi \left(A/{A}_{0}\right)$ 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

*κ*_{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 *κ*_{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 (van der Meer et al., 2007

), which indicates shear-thinning behavior. Based on those experimental data, Doinikov et al., 2009a

included shear-thinning behavior in the RP equation in an empirical way. Furthermore, Segers et al., 2018a

have measured a dependence of *κ*_{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 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., 2016b

, Segers et al., 2017

, Segers 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 blood-pressure 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

with

*σ*_{s}:${\sigma}_{\mathrm{s}}=4\pi {r}^{2}\frac{{\left|{P}_{\mathrm{s}}\left(r\right)\right|}^{2}}{{\left|{P}_{\mathrm{T}}\right|}^{2}},$

(14)

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 Ω

Equation (15) shows that for a bubble driven at resonance $\left(\omega ={\omega}_{0}\right)$ 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.

_{s}can be expressed as follows (de Jong et al., 1992

; ):${\Omega}_{\mathrm{s}}=\frac{4\pi {R}_{0}^{2}}{{\left({\left({\omega}_{0}/\omega \right)}^{2}-1\right)}^{2}+{\delta}_{\text{tot}}^{2}}.$

(15)

Equation (15) shows that for a bubble driven at resonance $\left(\omega ={\omega}_{0}\right)$ 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.

In addition to scattering, energy of the insonifying ultrasound wave is dissipated through 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*σ*_{e}. It is found by multiplying the scattering cross section*σ*_{s}and the ratio of the total damping*δ*_{tot}to the radiation damping*δ*_{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

with

*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}Ω_{s}≪ 1 (Commander and Prosperetti, 1989

; Goertz et al., 2007

). The total attenuation *α*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$\alpha =-\frac{10}{ln\left(10\right)}\sum _{R}\frac{{\delta}_{\text{tot}}}{{\delta}_{\text{rad}}}n\left(R\right){\Omega}_{\mathrm{s}}d,$

(16)

with

*d*the acoustic path length over which the sound wave interacts with the bubbles. Similarly, again for linear bubble oscillations, the backscatter coefficient*η*can be calculated as follows (Marsh et al., 1999

; Gorce et al., 2000

):$\eta =-\frac{1}{4\pi}\sum _{R}n\left(R\right){\Omega}_{\mathrm{s}}.$

(17)

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}. 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 interaction 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

). 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

with

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:

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:$R\left(t,\theta ,\varphi \right)=R\left(t\right)+{a}_{n}\left(t\right){Y}_{n}\left(\theta ,\varphi \right),$

(18)

with

*Yn*the spherical harmonic mode of order*n*and*an*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${\ddot{a}}_{n}+\frac{3\dot{R}}{R}{\dot{a}}_{n}-\left(\left(n-1\right)\frac{\ddot{R}}{R}-\frac{\left(n-1\right)\left(n+1\right)\left(n+2\right)\sigma}{{\rho}_{\mathrm{L}}{R}^{3}}\right){a}_{n}=0.$

(19)

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:

${\omega}_{0n}=\frac{\left(n-1\right)\left(n+1\right)\left(n+2\right)\sigma}{{\rho}_{\mathrm{L}}{R}_{0}^{3}}.$

(20)

In this potential flow approach, the velocity field is described by the gradient of the velocity potential

*∇ϕ*, 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 $\dot{R}=0$. This model was later picked up byCeschia 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 *δ*_{BL}depends on the bubble size and the mode number:${\delta}_{\text{BL}}=min\left(\sqrt{\frac{{\mu}_{\mathrm{L}}}{{\rho}_{\mathrm{L}}\omega}},\frac{{R}_{0}}{2n}\right).$

(21)

Including this short-range vorticity modifies the second term of eqn (19), where $3\dot{R}/R$ must be replaced by the damping term

*Bn*(*t*) given by eqn (13) ofHilgenfeldt 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 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 $\sqrt{2/3}$. Much later,

Doinikov et al., 2011

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 non-spherical 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

):${\rho}_{\mathrm{L}}\ddot{l}\left(L-l\right)={P}_{0}{\left(\frac{{l}_{0}}{l}\right)}^{\kappa}-P\left(L\right).$

(22)

After linearization around $l={l}_{0}\left(1+\epsilon \right)$ for small amplitude ε ≪ 1,

with eigenfrequency

$\left[{\rho}_{\mathrm{L}}{l}_{0}\left(L-{l}_{0}\right)\right]\ddot{\epsilon}+\kappa {P}_{0}\epsilon =-{P}_{\mathrm{A}},$

(23)

with eigenfrequency

${\omega}_{0}=\sqrt{\frac{\kappa {P}_{0}}{{\rho}_{\mathrm{L}}{l}_{0}\left(L-{l}_{0}\right)}}.$

(24)

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 effects of surface tension, the position *e*of the bubble along the tube axis, and tube radius*R*_{tube}(Fig. 8d). This model also provides a first-order 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 vessel-wall 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 a 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

using 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 on its effect in 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 Corato et al., 2019

investigated the dynamics of large bubbles (>50 μm) 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

.## Modeling higher-order effects

### 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 ${\delta}_{\text{BL}}=\sqrt{{\mu}_{\mathrm{L}}/{\rho}_{\mathrm{L}}\omega}$, with ω the angular frequency.

Rayleigh, 1884

was the first to describe acoustic streaming theoretically. He showed that an oscillatory flow of velocity *U*(*x*)cos*ωt*along a wall (*x*being a coordinate along the wall) causes a steady-state streaming flow of velocity*u*_{s}such that ${u}_{\mathrm{s}}\left(x\right)=-\left(3/4\right)U\left(x\right){U}^{\prime}\left(x\right)/\omega $ along the wall (*e.g.*, ); to be rigorous, this applies only if the streaming Reynolds number*u*_{s}*ρ*_{L}*ℓ*/*μ*_{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 *τxz*≈*μ*_{L}*u*_{s}/*δ*_{BL},*z*being the coordinate away from the wall. Hence, ${\tau}_{xz}\approx \sqrt{{\rho}_{\mathrm{L}}{\mu}_{\mathrm{L}}/\omega}\phantom{\rule{0.25em}{0ex}}U\left(x\right){U}^{\prime}\left(x\right)$.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 byde 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 μm in radius) (Fig. 10f). They showed different mechanisms of shedding, depending on lipid composition, from submicrometer 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 μm toward a stable radius of 1.4 μm, 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 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