Droplets, Bubbles and Ultrasound Interactions

Fig. 9.1

PFC-in-water, PFC-in-oil-in-water and oil-in-PFC-in-water emulsions under the microscope

PFC liquids are known for their use in medicine due to their biocompatibility and inertness (Biro et al. 1987). PFC nanodroplet emulsions can be utilized for selective extravasation in tumor regions (Long et al. 1978). Due to their biocompatibility and suggested ability to passively target regions of cancer growth, PFC droplets represent an attractive tool for cancer diagnosis. PFC droplets may also extravasate and be retained in the extravascular space due to the enhanced permeability and retention effect in tumors (Rapoport et al. 2007; Zhang and Porter 2010). Extravasated droplets may be acoustically converted into gas bubbles allowing for ultrasound tumor imaging. At the same time PFC droplets are rich in fluorine, which makes them potential candidates as a contrast agent for MRI imaging. The availability of both intravascular contrast agents (microbubbles), and tumor-specific extravascular contrast agents (nanodroplets), would significantly increase diagnostic and therapeutic capabilities. Moreover, the droplets may be used to deliver chemotherapeutic agents to tumor regions, and locally release them upon exposure to triggered ultrasound (Rapoport et al. 2009).

9.2 Nonlinear Propagation

The amplitude of the acoustic pressure that is required to nucleate droplets in ADV turns out to be very high (Kripfgans et al. 2000). To obtain a sufficiently high pressure, a focused ultrasound transducer is applied and the droplet is placed in the focal area of the emitted beam. Moreover, the frequency of the emitted ultrasound wave is several MHz. In a typical ADV experiment, the ultrasound wave travels a few centimeters (Kripfgans et al. 2000; Reznik et al. 2013; Shpak et al. 2013a, b; Giesecke and Hynynen 2003; Schad and Hynynen 2010; Williams et al. 2013) before impinging on the droplet. The high pressure, high frequency, applied focusing and long propagation distances are all factors that strengthen the nonlinear behavior of the ultrasound wave (Blackstock 1964; Bacon 1984). As a result, the wave that impinges on the droplet will be a highly deformed version of the one that is emitted by the transducer (Fig. 9.2). This has important consequences for the focusing inside the droplet, as will be demonstrated in Sect. 9.4.2.2.

Fig. 9.2

Schematics of nonlinear propagation of an ultrasound wave. T and f are the period of oscillation and frequency of an ultrasound wave, respectively. Two upper plots are the transmitted and propagated wave, and two lower plots are their frequency domains. –fft– stands for the Fast Fourier Transform

9.2.1 Basic Equations for the Nonlinear Ultrasound Beam

Similar to most cases involving nonlinear medical ultrasound, the description of the beam that hits the droplet can be based on the Westervelt equation (Westervelt 1963; Hamilton and Morfey 2008):

${\nabla}^2p-\frac{1}{c_0^2}\frac{\partial^2p}{\partial {t}^2}+\frac{d}{c_0^4}\frac{\partial^3p}{\partial {t}^3}=-\frac{b}{r_0{c}_0^4}\frac{\partial^2{p}^2}{\partial {t}^2}$

(9.1)

where ${\nabla}^2={\partial}^2/\partial {x}^2+{\partial}^2/\partial {y}^2+{\partial}^2/\partial {z}^2$ is the Laplace operator and $p=p\left(x,y,z,t\right)$ denotes the acoustic pressure. The medium in which the ultrasound wave propagates is characterized by the ambient speed of sound ${c}_0$ , the ambient density of mass ${\rho}_0$ , the diffusivity of sound $\delta$ and the coefficient of nonlinearity $\beta$ . Unfortunately, closed-form analytical solutions of this equation do not exist and its numerical solution generally requires considerable computational effort. However, in the present case of a narrow, focused beam and a homogeneous medium, some simplifying assumptions can be made. Firstly, it may be assumed that the predominant direction of propagation is along the transducer axis, which is taken in the z -direction. In this case, we can replace the ordinary time coordinate $$ t $$

by the retarded time coordinate $t=\left(t-{t}_0\right)-\left(z-{z}_0\right)/{c}_0$ , which keeps the same value when traveling along with the wave. Here, ${t}_0$ is the time at which the transducer emits the pressure wave, and ${z}_0$ is the axial position of the transducer. The equivalent of Eq. 9.1 in the co-moving time frame is:

${\nabla}^2\overline{p}-\frac{2}{c_0}\frac{\partial^2\overline{p}}{\partial z\partial \tau }+\frac{d}{c_0^3}\frac{\partial^3\overline{p}}{\partial {\tau}^3}=-\frac{b}{r_0{c}_0^4}\frac{\partial^2{\overline{p}}^2}{\partial {\tau}^2}$

(9.2)

with $\overline{p}=\overline{p}\left(x,y,z,t\right)$ denoting the acoustic pressure in the co-moving time frame. Secondly, it may be assumed that in the retarded time frame the axial derivative ${\partial}^2\overline{p}/\partial {z}^2$ is much smaller than the lateral derivatives ${\partial}^2\overline{p}/\partial {x}^2$ and ${\partial}^2\overline{p}/\partial {y}^2$ . This motivates the use of the parabolic approximation ${\nabla}^2\overline{p}\approx {\nabla}_{\perp}^2\overline{p}$ , where ${\nabla}_{\perp}^2={\partial}^2/\partial {x}^2+{\partial}^2/\partial {y}^2$ is the Laplace operator in the lateral plane. This approximation is valid for waves propagating under at most 20° of the transducer axis (Lee and Pierce 1995). Applying the parabolic approximation to Eq. 9.2 and rearranging terms results in the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation (Zabolotskaya and Khokhlov 1969; Kuznetsov 1971):

$\frac{\partial^2\overline{p}}{\partial z\partial \tau }=\frac{c_0}{2}{\nabla}_{\perp}^2\overline{p}+\frac{\delta }{2{c}_0^3}\frac{\partial^3\overline{p}}{\partial {\tau}^3}+\frac{b}{2{r}_0{c}_0^3}\frac{\partial^2{\overline{p}}^2}{\partial {\tau}^2}$

(9.3)

Dedicated coordinate transformations may be applied to improve the numerical solution in the far field (Hamilton et al. 1985; Hart and Hamilton 1988) or to adapt to specific forms of focused beams (Kamakura et al. 2000), but these will not be discussed here.

9.2.2 Numerical Solution for the Nonlinear Ultrasound Beam

We will follow a well-known numerical solution strategy (Lee and Hamilton 1995; Cleveland et al. 1996) that is based on the time integrated version of Eq. 9.3:

$\frac{\partial \overline{p}}{\partial z}=\frac{c_0}{2}\underset{-\infty }{\overset{t}{{\displaystyle \int }}}{\nabla}_{\perp}^2\overline{p}\left({\tau}^{\prime}\right)d{\tau}^{\prime }+\frac{\delta }{2{c}_0^3}\frac{\partial^2\overline{p}}{\partial {\tau}^2}+\frac{b\overline{p}}{r_0{c}_0^3}\frac{\partial \overline{p}}{\partial \tau }$

(9.4)

The first term at the right-hand side of this equation accounts for the diffraction of the beam, the second term for its attenuation, and the third term for its nonlinear distortion. Further, the solution strategy is based on the split-step approach. This means that the field $\overline{p}$ is stepped forward over a succession of parallel planes with mutual distance $\varDelta \kern0.1em z$ , where the field $\overline{p}\left(x,y,{z}_0,\tau \right)$ in the transducer plane acts as the starting plane. The stepsize $\varDelta z$ is taken sufficiently small, allowing that each of the above phenomena may be accounted for in separate sub steps (Varslot and Taraldsen 2005). Therefore, the total step $z\to z+\varDelta z$ involves the numerical solution of the separate equations:

$\frac{\partial \overline{p}}{\partial z}=\frac{c_0}{2}\underset{-\infty }{\overset{t}{{\displaystyle \int }}}{\nabla}_{\perp}^2\overline{p}\left(\tau^{\prime}\right)d\tau^{\prime }$

(9.5)

$\frac{\partial \overline{p}}{\partial z}=\frac{\delta }{2{c}_0^3}\frac{\partial^2\overline{p}}{\partial {\tau}^2}$

(9.6)

$\frac{\partial \overline{p}}{\partial z}=\frac{b\overline{p}}{\rho_0{c}_0^3}\frac{\partial \overline{p}}{\partial \tau }$

(9.7)

over the same interval, where the result of solving one equation is used as the input for solving the next one. A numerical implementation of the above process is used to step the acoustic pressure from the transducer to the focus of the beam, i.e. the location of the droplet. For convenience, it is now assumed that the droplet is located at the origin of the coordinate system and that the source emits the pressure wave at ${t}_0={z}_0/{c}_0$ . This makes $\tau =t$ at the position of the droplet. For ease of notation, the bar and the coordinates of the droplet will be suppressed, and the pressure at the droplet position, as obtained from the numerical solution of the KZK-equation, will simply be indicated by ${p}_{KZK}(t)$ .

9.2.3 Nonlinear Pressure Field at the Focus of the Beam

The nonlinear pressure field at the focus of the beam can be expanded in a Fourier series:

$\begin{array}{c}{p}_{KZK}(t)={\displaystyle \sum}_{n=0}^{\infty }{a}_n \cos \left(nwt+{f}_n\right)\\ {}=Re\left[{\displaystyle \sum}_{n=0}^{\infty }{a}_n{e}^{i\left(nwt+{f}_n\right)}\right]\end{array}$

(9.8)

where ${a}_n$ and ${f}_n$ are the amplitudes and the phases of the $$ n-th $$

harmonic component of the ultrasound wave. For convenience, all the subsequent derivations will be given in the complex representation, so we will omit taking the real part and simply write:

${p}_{KZK}(t)={\displaystyle \sum}_{n=0}^{\infty }{a}_n{e}^{i\left(n\omega t+{\phi}_n\right)}$

(9.9)

Given that nonlinear deformation of the waveform builds up over distance and the droplet is four orders of magnitude smaller in size than the distance to the transducer, the additional nonlinear distortion inside the droplet is neglected. This implies that wave propagation inside the droplet is considered linear, so the superposition theorem holds, and the focusing of each harmonic component in the droplet may be analyzed on an individual basis, as will be done in Sect. 9.4.2.2.

9.3 Bubble Dynamics

9.3.1 Dynamics of a Gas Bubble

Bubble radial oscillations are governed by the Rayleigh-Plesset equation:

$\ddot{R}R+\frac{3}{2}{\dot{R}}^2=\frac{\varDelta P}{\rho }$

(9.10)

where

, $\dot{R}$ , and $\ddot{R}$ are the radius, the velocity and the acceleration of the bubble wall, respectively, and $\rho$ is the density of the liquid. $\varDelta P={P}_L(R)-{P}_{\infty }$ is the pressure difference between the liquid at the bubble wall ${P}_L(R)$ and the external pressure infinitely far from the bubble ${p}_{\infty }$ . Equation 9.10 was first described by Lord Rayleigh (1917) for the case $\varDelta P=0$ and was later refined (Plesset 1949; Noltingk and Neppiras 1950; Neppiras and Noltingk 1951; Poritsky 1952). It is derived for a spherically symmetric bubble, and follows from the Bernoulli’s equation and the continuity equation (Leighton 1994). Equation 9.10 assumes spherical symmetry of the bubble, and the motion of the liquid around the bubbles is considered to be spherically symmetric. The liquid is incompressible.

The bubble is assumed to be much smaller than the acoustic wavelength, such that acoustic pressure is considered to be uniform. Thus, the pressure at infinity is the sum of the acoustic forcing $$ P(t) $$

and the ambient pressure ${P}_0$ :

${p}_{\infty }=P(t)+{P}_0$

(9.11)

The interfacial pressure acting on the liquid at the bubble wall consists of the Laplace pressure $2\sigma /{R}_0$ , viscous pressure $4\mu \dot{R}/R$ and the gas pressure ${P}_g$ . Neglecting the vapor pressure of the liquid, the gas pressure inside the bubble as a function of the bubble radius $$ R $$

can be described by the ideal gas relation ${P}_g{V}^{\gamma }= const$ , where $\gamma$ is the polytropic constant and $V\propto {R}^3$ is the bubble volume. For this derivation we first neglect the gas diffusion. Thus, the total number of gas molecules inside the bubble is constant. In equilibrium, the pressure inside the bubble ${P}_{eq}$ is equal to the sum of the ambient pressure ${P}_0$ and the Laplace pressure:

${P}_{eq}={P}_0+\frac{2\sigma }{R_0}$

(9.12)

where $\sigma$ and ${R}_0$ are the surface tension and the equilibrium radius, respectively. In combination with the ideal gas law, the dependence of the gas pressure as a function of the bubble radius can be written as ${P}_g=\left({P}_0+\frac{2\sigma }{R_0}\right){\left(\frac{R_0}{R}\right)}^{3\gamma }$ . The right-hand side of Eq. 9.10 can then be written as:

$\begin{array}{c}DP=\left({P}_0+\frac{2s}{R_0}\right){\left(\frac{R_0}{R}\right)}^{3\gamma }-{P}_0-\frac{2s}{R}\\ {}-4\mu \frac{\dot{R}}{R}-P(t)\end{array}$

(9.13)

which gives the final form of the bubble dynamic equation:

$\begin{array}{c}r\left(\ddot{R}R+\frac{3}{2}{\dot{R}}^2\right)=\left({P}_0+\frac{2s}{R_0}\right){\left(\frac{R_0}{R}\right)}^{3g}\\ {}-{P}_0-\frac{2s}{R}-4m\frac{\dot{R}}{R}-P(t)\end{array}$

(9.14)

The microbubbles in the ultrasound contrast agents can be encapsulated with a phospholipid, protein or polymer coating, thus preventing bubbles from dissolution. For more details please see (Marmottant et al. 2005; de Jong et al. 2007; Church 1995). The viscoelastic coating also contributes to an increased stiffness and to additional viscous damping (Overvelde et al. 2010).

9.3.2 Linearization

The acoustic pressure typically has the form of a sinusoidal oscillation $P(t)={P}_A \sin \left(\omega t\right)$ , with ${P}_A$ being the driving pressure amplitude, and $\omega$ the driving pressure angular frequency. With relatively small oscillation amplitudes Eq. 9.10 can be linearized. To rewrite Eq. 9.10 in linear terms we express the bubble radius $$ R $$

as:

$R={R}_0\left(1+x\right)$

(9.15)

with ${R}_0$ the equilibrium radius, as before, and $x\ll 1$ a small dimensionless perturbation to the radius. Substituting Eq. 9.15 into Eq. 9.10 and retaining only first-order terms, $$ x $$

, $\dot{x}$ and $\ddot{x}$ gives:

$\ddot{x}+2b\dot{x}+{w}_0^2x=\frac{P_A}{r{R}_0^2} \sin \left(w\kern0.1em t\right)$

(9.16)

where

${w}_0=\sqrt{\frac{1}{r{R}_0^2}\left[3g\left({P}_0+\frac{2s}{R_0}\right)-\frac{2s}{R_0}\right]}$

(9.17)

the eigenfrequency of bubble oscillations, and:

$b=\frac{2m}{r{R}_0^2}$

(9.18)

the damping due to viscosity. The damping has the dimensions of the reversed time $\left[{s}^{-1}\right]$ and represents how fast the amplitude of oscillations is decaying in time due to the energy loss.

The solution to the equation Eq. 9.16 is:

$x(t)={X}_t{e}^{-bt} \cos \left({w}_1t\right)+{X}_s \cos \left( wt+{f}_1\right)$

(9.19)

with the ${\phi}_1$ being the phase shift between the two terms:

${f}_1= \arctan \left(\frac{w_0^2-{w}^2}{2bw}\right)$

(9.20)

The first term of Eq. 9.19 is the transient solution. Its amplitude dampens out in time as ${X}_t{e}^{-bt}$ , where ${X}_t$ is the amplitude of transient oscillations at time ${t}_0=0$ . Not only the viscosity of water can contribute to the damping, but also the acoustic reradiation and the viscosity of the coating shell and thermal damping. For more details please see (Overvelde et al. 2010). The frequency of the transient solution is equal to ${w}_1=\sqrt{w_0^2-{b}^2}$ . The amplitude of the transient solution ${X}_t$ depends strongly on the initial conditions.

The second term of Eq. 9.19 is the steady-state solution. The amplitude of the steady-state response depends on the driving frequency as:

${X}_s=\frac{P_A}{r\kern0.1em {R}_0^2}\frac{1}{\sqrt{{\left({w}_0^2-{w}^2\right)}^2+4{w}^2{b}^2}}$

(9.21)

The resonance frequency ${\omega}_{\mathrm{res}}$ of the system, by definition, corresponds to the maximal amplitude of the steady-state solution. ${X}_s$ is at maximum, when the denominator in the Eq. 9.21 is at minimum. Thus, the resonance frequency relates to the eigenfrequency ${\omega}_0$ as:

${w}_{\mathrm{res}}=\sqrt{w_0^2-2{b}^2}$

(9.22)

The smaller the damping $$ b, $$

the closer the resonance frequency to the eigenfrequency of the bubble oscillations. Additionally, for large bubbles, when the Laplace pressure is small compared to the ambient pressure Eq. 9.17 simplifies to:

${w}_M=2p{f}_M=\sqrt{3g\frac{P_0}{r{R}_0^2}}$

(9.23)

with ${f}_M$ the Minnaert eigenfrequency, resonance frequency of the bubble (Minnaert 1933). Relation Eq. 9.23 tells us that the resonance frequency can be estimated directly from the bubble radius ${R}_0$ . For a bubble in water at standard pressure ( ${P}_0=100kPa,\;r=1000 kg/{m}^3$ ), the equation becomes ${f}_M{R}_0\approx 3.26\;m\kern0.1em m. MHz$ . The smaller the bubble radius, the higher the resonance frequency becomes.

It is insightful to make the analogy to the classical mass-spring system. The dynamics of the classical mass spring system is governed by the equation:

$\ddot{x}+2\frac{b^{\prime }}{m}\dot{x}+{w}_0^{\prime 2}x=\frac{F_0}{m} \sin (wt)$

(9.24)

where $w{\prime}_0=\sqrt{k/m}$ is the eigenfrequency, ${\beta}^{\prime }$ is the damping constant, $$ k $$

is the spring stiffness, ${F}_0$ is the driving force and $$ m $$

is the mass. Equation 9.24 has the same form as Eq. 9.16. Thus, a gas inside the bubble, represented by the polytropic constant $$ g $$

, acts as the restoring force, the liquid around the bubble acts as a mass $\left(4p{R}_0^3r\right)$ , and ultrasound is acting as a driving force $\left(12 pg\kern0.1em {R}_0{P}_0\right)$ .

9.3.3 Pressure Emitted by the Bubble

Far from the bubble wall, at a distance $$ r $$

, the velocity of the liquid ${v}_r$ can be calculated from the continuity equation (Prosperetti 2011):

$4p{r}^2{v}_r=4p{R}^2\dot{R}$

(9.25)

${v}_r=\frac{R^2}{r^2}\dot{R}$

(9.26)

The liquid is incompressible and the bubble wall and the liquid motion around the bubble are spherically symmetric. $\dot{R}$ is the velocity of the bubble wall in the radial direction, as before.

The pressure field, generated by the radial bubble wall oscillations, can be calculated from the Euler equation (Prosperetti 2011):

$r\frac{\partial v}{\partial t}+\frac{\partial p}{\partial r}=0$

(9.27)

where

is the pressure emitted by the bubble. In Eq. 9.27 we omit the nonlinear convective term. Substituting the expression of the velocity field (Eq. 9.26) into Eq. 9.27 gives the pressure gradient:

$\frac{\partial p}{\partial r}=-\frac{r}{r^2}\frac{d}{dt}\left({R}^2\dot{R}\right)$

(9.28)

and the pressure emitted by the bubble:

$p=\frac{r}{r}\frac{d}{dt}\left({R}^2\dot{R}\right)=r\left(\frac{R^2\ddot{R}+2R{\ddot{R}}^2}{r}\right)$

(9.29)

9.3.4 Secondary Bjerknes Force

Let us now consider two interacting gas bubbles, separated by a distance $$ l $$

. The distance between the bubbles 1 and 2 is much larger that their radii ${R}_1(t)$ and ${R}_2(t)$ . Thus, we can consider the motion of the liquid around the bubbles to be spherically symmetric. Bubble 2 with volume ${V}_2=\frac{4}{3}p{R}_2^3$ experiences a force ${F}_{12}$ as a result of the pressure emitted by bubble 1, ${p}_1$ (Leighton 1994):

${F}_{12}=-{V}_2\nabla {p}_1$