The influence of medium elasticity on the prediction of histotripsy-induced bubble expansion and erythrocyte viability

The radial oscillations of cavitation nuclei in soft tissue were calculated by numerical integration of the Yang and Church model (2005). An adaptive fourth-order Runga-Kutta algorithm was implemented in MATLAB^® (The Mathworks, Natick, MA, USA), employing the function 'ode15s' to carry out the integration. The Yang/Church model has the form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(1-\frac{{\dot{R}}}{c} \right)R\ddot{R}+\left(1-\frac{{\dot{R}}}{c} \right)\frac{3}{2}{{\dot{R}}^{2}}=\left(1+\frac{{\dot{R}}}{c} \right)\frac{{{p}_{w}}\left(\dot{R},~R,~t \right)}{\rho }+\frac{R}{\rho c}\frac{\partial {{p}_{w}}}{\partial t}\nonumber \end{align} \tag{ 1 }$$

where R is the time dependent bubble radius, the diacritic dot denotes the temporal derivative, c is the sound speed of the viscoelastic medium (1540 m s⁻¹), and ρ is the medium density (1000 kg m⁻³). The pressure at the bubble wall, p_w, is defined in terms of the properties of the viscoelastic medium (Yang and Church 2005):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{w}}=\left({{P}_{0}}+\frac{2\sigma }{{{R}_{0}}} \right){{\left(\frac{{{R}_{0}}}{R} \right)}^{3\kappa }}-\frac{2\sigma }{R}-\frac{4\mu \dot{R}}{R}-\frac{4~G}{3}\left[ 1-{{\left(\frac{{{R}_{0}}}{R} \right)}^{3}} \right]-{{P}_{0}}-{{p}_{AC}}(t)\nonumber \end{align} \tag{ 2 }$$

where P₀ is the ambient pressure (0.1 MPa), R₀ is the initial radius of the air-filled cavitation nucleus, $\kappa$ is the ratio of specific heats (1.4), and p_AC is the time-dependent source term. The following values were used for the medium properties: surface tension, σ  =  0.056 N m⁻¹ (Holland and Apfel 1989, Church et al 2015), and dynamic viscosity, µ  =  0.005 kg (m · s)⁻¹ (Holland and Apfel 1989, Church et al 2015). Medium viscoelasticity was computed via the Kelvin–Voight model assuming an isotropic medium, which has been validated for soft tissues in the 3–13 MHz range (Madsen et al 1983). Bubble expansion during histotripsy is slow compared to the sound speed (Bader and Holland 2016), and the shear modulus G can be expressed in terms of the elastic modulus E for incompressible media as $G~\approx ~3E$. Elastic moduli between 1 kPa and 20 GPa were considered in this study, which spans the properties of soft and mineralized tissues (Duck 1990).

Shock-scattering histotripsy and microtripsy excitations were considered in this study. The maximum bubble diameter is independent of the nucleus diameter for cavitation nuclei larger than 5 nm (microtripsy) or 20 nm (shock-scattering histotripsy) (Bader and Holland 2016). Calculations in this study utilized these nuclei sizes for each respective form of histotripsy to provide an upper estimate to the maximum bubble size. The time-dependent source term in (2), p_AC, was approximated as a half cycle, peak negative only pulse for microtripsy excitations (Vlaisavljevich et al 2015b, Mancia et al 2017), with peak negative pressures of 25–50 MPa, and fundamental frequencies of 0.1–3 MHz. Shock-scattering histotripsy pulses were implemented with pressure waveforms measured experimentally of 1 MHz fundamental frequency source as described previously (Bader and Holland 2016).

The maximum size of a cavitation nucleus exposed to a single-cycle excitation was derived analytically via Holland and Apfel (1989):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{MAX}}~=~\left[ {{R}_{0}}+\sqrt{\frac{2{{P}_{0}}\xi }{9\rho }}\tau \right]{{\left[ \frac{\xi {{P}_{0}}}{3{{p}_{EFF}}}+1 \right]}^{1/3}}\nonumber \end{align} \tag{ 3 }$$

where $\xi$ and τ are defined in Holland and Apfel (1989), and P₀, R₀, and ρ are as defined in section 2.1. The effective pressure, p_EFF, acts as a 'Blake brake' to arrest bubble expansion (Bader and Holland 2016), and depends on the histotripsy insonation scheme. Microtripsy pulses do not have an appreciable compressional phase (Maxwell et al 2013, Lin et al 2014, Vlaisavljevich et al 2015a), and the effective pressure can be set to the ambient pressure, p_EFF  =  P₀ (Bader and Holland 2016). For shocked pulses, the bubble wall velocity is maximized at the shock front, and the effective pressure can be approximated as the magnitude of the negative pressure, p_EFF  =  |p_N|.

Incorporation of elasticity into the analytic model requires two modifications: (1) a delay to bubble expansion (Apfel 1986, Holland and Apfel 1989), and (2) an additional overpressure as the bubble expands (Gent and Lindley 1959). Bubble expansion will lag behind the onset of tension due to conformational changes of the elastic medium (Vincent 2012). Following Apfel (1981b) and Holland and Apfel (1989), the delay to bubble expansion due to elasticity, τ_G, can be estimated by equating the elastic and driving terms in the Yang/Church model (fourth and sixth term in (2)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\tau }_{G~}}=~4G\frac{{{R}_{0}}}{3\xi {{P}_{0}}}\sqrt{\frac{9\rho }{2{{P}_{0}}\xi }}.\nonumber \end{align} \tag{ 4 }$$

Delays in bubble expansion due to elasticity can be integrated into the analytic model through modifying the term τ in (3) as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tau ={{\tau }_{HA}}-{{\tau }_{G~}}\nonumber \end{align} \tag{ 5 }$$

where τ_HA is as defined by equation (9) in Holland and Apfel (1989), which includes the effects of surface tension, viscosity, and inertia.

Elasticity also modifies the effective pressure (Gent and Lindley 1959). This is denoted in figure 1, where the pressure terms due to surface tension, viscosity, and elasticity are plotted below for a bubble exposed to the first cycle of a shocked histotripsy pulse (second, third, and fourth terms, respectively, in (2)). As noted in figure 1(B), the pressure contributions from surface tension and viscosity decrease rapidly during bubble expansion. The elastic pressure contribution asymptotes to a fixed value of 4G/3 for the finite-strain Kelvin–Voight elastic model (the limit of R₀/R $\ll$ 1 for the fourth term in (2))

Zoom In Zoom Out Reset image size

The additional overpressure modifies the parameters $~\xi$ and p_EFF in the analytic model. The former represents the time-averaged pressure as the bubble undergoes expansion. Substituting P₀  +  4G/3 for P₀ into equation (4) of Holland and Apfel (1989), $\xi$ can be estimated to second order as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \xi ~=~{{\xi }_{HA}}-~\frac{4G}{{{P}_{0}}}\nonumber \end{align} \tag{ 6 }$$

where ${{\xi }_{HA}}$ is as defined by equation (6) in Holland and Apfel (1989). Modification to the effective pressure term depends on the insonation type (Bader and Holland 2016). The elastic pressure term acts to increase the overall ambient pressure (figure 1(B)). For a microtripsy excitation, the additional overpressure acts to increase the effective pressure (p_EM):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{EM}}~=~{{P}_{0}}+\frac{4G}{3}.\nonumber \end{align} \tag{ 7 }$$

For a shock-scattering histotripsy excitation, the elastic pressure acts to decrease the overall wall velocity and the effective pressure takes the form (p_ES):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{p}_{ES}}~=~\left| {{p}_{N}} \right|~-~\frac{4G}{3}.\nonumber \end{align} \tag{ 8 }$$

The analytic model (3) predicts the bubble size after a single-cycle excitation. A shock-scattering histotripsy pulse is 3–20 cycles in duration (Maxwell et al 2012), during which time the bubble undergoes continuous expansion (Maxwell et al 2011b). The bubble expansion rate is largest during the first cycle of the histotripsy pulse. Beyond the first cycle, the expansion rate slows due to the shocked compressional pressure. The slowed bubble expansion rate after the first cycle of the histotripsy pulse $\left\langle {{{\dot{R}}}_{s}} \right\rangle$ can be approximated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\langle {{{\dot{R}}}_{S}} \right\rangle =~\sqrt{\frac{2\left(4G/3-{\langle{p}_{AC}\rangle} \right)}{3\rho }}.\nonumber \end{align} \tag{ 9 }$$

The value $\langle {{p}_{AC}}\rangle$ is the pressure of the histotripsy pulse, time-averaged from the second cycle to the end of the pulse. Equation (9) is the same form as equation (9) in Bader and Holland (2016), with the addition of an elastic term. The final extent of the bubble cloud, L_AX, can be computed utilizing (3) and (9) following Bader and Holland (2016).

The functional form of tissue elasticity is an active area of research (Schmitt et al 2011). For the massive bubble expansion characteristic of histotripsy-induced cavitation (R_MAX $\gg$ R₀), alternative forms of finite-strain elastic models can readily be incorporated into the analytic model using the same formalism described by (4)–(9). Table 1 denotes modifications to the effective pressure, $\xi$, and the time-averaged bubble wall velocity for Kelvin–Voight, linear, and neo-Hookean elastic models (Gent 2005b, Gaudron et al 2015).

Table 1. Modifications for the parameter $\xi$, the effective pressure for microtripsy (p_EM) and shock scattering histotripsy (p_ES), and bubble expansion rate ($\left\langle {\dot{R}} \right\rangle$) based on form of the elastic model.

Elasticity model	$\xi$	pEM	pES	$\left\langle {\dot{R}} \right\rangle$
Linear	ξHA  −  6G/P0	${{P}_{0}}\left(1+\frac{2G}{{{P}_{0}}} \right)$	$\left| {{p}_{N}} \right|-2G$	$\sqrt{\frac{2\left(2G-{\langle{p}_{AC}\rangle} \right)}{3\rho }}$
Kelvin-Voight	ξHA  −  4G/P0	${{P}_{0}}\left(1+\frac{4G}{3{{P}_{0}}} \right)$	$\left| {{p}_{N}} \right|-\frac{4G}{3}$	$\sqrt{\frac{2\left(4G/3-{\langle{p}_{AC}\rangle} \right)}{3\rho }}$
Neo-Hookean	ξHA  −  15G/2P0	${{P}_{0}}\left(1+\frac{5G}{2~{{P}_{0}}} \right)$	$\left| {{p}_{N}} \right|-\frac{5G}{2}$	$\sqrt{\frac{2\left(5G/2-{\langle{p}_{AC}\rangle} \right)}{3\rho }}$
Funga	ξHA  −  58.2G/P0	${{P}_{0}}\left(1+\frac{58.2G}{{{P}_{0}}} \right)$	$\left| {{p}_{N}} \right|-58.15G$	$\sqrt{\frac{2\left(58.2G-{\langle{p}_{AC}\rangle} \right)}{3\rho }}$
Gentb	ξHA  −  17.5G/P0	${{P}_{0}}\left(1+\frac{17.5G}{{{P}_{0}}} \right)$	$\left| {{p}_{N}} \right|-17.5G$	$\sqrt{\frac{2\left(17.5G-{\langle{p}_{AC}\rangle} \right)}{3\rho }}$

^aThe solution for the Fung model was computing assuming fracture of an agar medium beyond an expansion ratio of 2.2 (i.e. solution of (10) for λ  =  2.2) (Movahed et al 2016). ^bThe solution for the Gent model was computed assuming fracture of an agar medium beyond an expansion ratio of 2.93 (i.e. solution of (11) for λ  =  2.93).

In addition to the consideration of bubble behavior in soft tissues, in vitro studies describing histotripsy-induced bubble behavior are frequently conducted in agar-based phantoms (Maxwell et al 2010, Vlaisavljevich et al 2013b, Bader et al 2018). Agar behaves as a strain-hardening material (Barrangou et al 2006), and the elastic behavior can be modeled via Fung (Fung et al 1979, Movahed et al 2016) or Gent (Gent 1996). In the Fung model, the elastic pressure term (fourth term in (2)) can be computed as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{G}}~=~2G\int_{1}^{\lambda }{\left[ {{\beta }^{-5}}+~{{\beta }^{-2}} \right]{\rm exp}\left[ \alpha \left(2{{\beta }^{2}}+{{\beta }^{-4}}-3 \right) \right]{\rm d}\beta }\nonumber \end{align} \tag{ 10 }$$

where λ  =  R/R₀, and α ~ 1 (Movahed et al 2016). Gent elasticity can be expressed as (Gent 2005a):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{G}}~=~2G\int_{1}^{\lambda }{\left[ {{\beta }^{-5}}+~{{\beta }^{-2}} \right]1/\left[ 1-\left(2{{\beta }^{2}}+{{\beta }^{-4}}-3 \right)/{{J}_{M}} \right]{\rm d}\beta }.\nonumber \end{align} \tag{ 11 }$$

The term J_M in (11) represents the effect of finite chain extensibility (Roland 2011). For agar phantoms, the Gent model accurately predicts pure strain measurements of agar phantoms (1%–2.5% w/v) when J_M ~ 17.5 (appendix) (Barrangou et al 2006).

Unlike finite-strain elastic models, the elastic pressure contribution for a strain-hardening elastic model does not asymptote to a fixed value for massive bubble expansion (figure 2). Agar is a brittle material, and will fracture under sufficient strain (Madsen et al 2003, Wilks et al 2015). The elastic energy released by bubble expansion will then equate to the fracture energy required for the increased bubble size (Gent 1999). Thus, P_G will not increase beyond a critical fracture expansion λ_F for strain-hardening elastic models (Griffith 1921, Hutchens et al 2016). The onset of medium fracture depends on the form of the elastic model, and can be estimated using Griffith's criterion (Gent and Wang 1991). For the Fung elastic model, Griffith's criterion is satisfied when λ_F ~ 2.2 (Movahed et al 2016), or P_G ~ 58.2G as calculated via (10). Griffith's criterion is satisfied when λ_F ~ 2.93, or P_G ~ 17.5G for the Gent elastic model (Horgan and Saccomandi 1999).

Zoom In Zoom Out Reset image size

The viability of red blood cells due to transient strain from laser-induced cavitation can be expressed as Li et al (2013):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm Viability }~ =\frac{100\%}{1+{\rm exp}[-(S-{{S}_{0}})/{{S}_{1}}]}.\nonumber \end{align} \tag{ 12 }$$

where S is the percent strain, S₀  =  105.8%, and S₁  =  0.18% (Li et al 2013). The normal Hencky (or true) strain in an isotropic medium surrounding an expanded or collapsed bubble was computed in these studies as Mancia et al (2017):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S~=\frac{2}{3}{\rm log}\left(1-\frac{R_{i}^{3}-R_{0}^{3}}{{{r}^{3}}} \right)\nonumber \end{align} \tag{ 13 }$$

where r is the radial distance from the center of the bubble and R_i corresponds to the maximum or minimum bubble radius. The maximum bubble radius was computed analytically via (3) for a Kelvin–Voight elastic medium (table 1). The minimum bubble size was computed assuming the bubble behaved as a Rayleigh cavity (Rayleigh 1917, Leighton 1994) with an additional contribution due to Kelvin–Voight elasticity:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{R}_{MIN}}~=~{{\left[ \frac{{{R}_{MAX}}}{\kappa -1}\frac{{{P}_{0}}+\frac{2\sigma }{{{R}_{0}}}}{{{P}_{0}}+\frac{4G}{3}}{{\left(\frac{{{R}_{0}}}{{{R}_{MAX}}} \right)}^{3}} \right]}^{1/3\kappa -3}}.\nonumber \end{align} \tag{ 14 }$$

Histotripsy-induced cavitation expands to sizes much larger than that of the original nucleus (Maxwell et al 2013, Bader and Holland 2016), meaning the bubble will behave as a Rayleigh cavity until the final stages of bubble collapse (Lastman and Wentzell 1981). Dampening will soften the bubble collapse, and the predictions of minimum bubble radius by (14) will be an underestimate compared to the true size (Prosperetti and Lezzi 1986). The large bubble expansion induced by histotripsy pulses necessitates that the bubble size at the completion of the collapse will be much smaller than the size of the initial nucleus (Gaudron et al 2015). Errors in the strain surrounding the bubble at its minimum radius as calculated via (13) will be less than 0.1%.

Delays to bubble expansion due to surface tension, viscosity, elasticity, and medium inertia are shown in figure 3. Surface tension dominated the onset of bubble expansion for nuclei smaller than 5 nm (figure 3(A)) and for peak negative pressures less than 8 MPa (figure 3(B)). For the operational range of histotripsy (peak negative pressures greater than 15 MPa) in soft tissues (elastic modulus less than 1 MPa (Duck 1990)), viscous and inertial effects dominated bubble expansion (figure 3(B)). As the elastic modulus increased beyond 105 MPa, medium elasticity was the largest contributor to delaying bubble expansion (figure 3(C)).

Zoom In Zoom Out Reset image size

Maximum bubble sizes based on numerical computations (Yang/Church model) and analytic calculations are shown in figure 4 as a function of the medium elasticity. Numerical computations were conducted with microtripsy or single-cycle shock-scattering histotripsy excitations, and a Kelvin–Voight elastic model. Analytic calculations were conducted without the contribution of elasticity (i.e. G  =  0), and with Kelvin–Voight elasticity incorporated into the analytic model (table 1). For elastic moduli less than 1 MPa, the numeric and analytic calculations agreed within 1.5% for shock-scattering histotripsy and microtripsy excitations. The analytic model underestimated the maximum bubble size relative to the numerical computation over the elastic modulus range of 1.8 and 20 MPa for microtripsy excitations.

Zoom In Zoom Out Reset image size

The parameters $\xi$ and p_ES become negative when the elastic moduli exceeds approximately 20 MPa for both shock-scattering histotripsy and microtripsy. Mathematically, this means that the predictions of the analytic model were complex. Physically, this means the tension of the histotripsy pulse is insufficient to cause expansion of the cavitation nucleus.

The analytic prediction for the bubble expansion rate over the duration of a shock-scattering histotripsy pulse is shown in comparison the numerical computation in figure 5(A). The analytic model agreed within 2.5% of the numerical computations for tissue with an elastic modulus less than 875 kPa and greater than 2.5 MPa. For tissue with an elastic modulus between 875 kPa and 2.5 MPa, the analytic model overestimated the bubble expansion rate compared to estimates based on numerical computations. The bubble expansion rate was only weakly dependent on medium elasticity when the elastic modulus was less than 115 kPa for both numeric and analytic calculations. For stiffer media, the bubble growth rate decreased rapidly to 0 µm/cycle.

Zoom In Zoom Out Reset image size

The analytic computation of the bubble cloud axial length is shown in figure 5(B) for a 5 µs duration shock-scattering histotripsy pulse (1 MHz fundamental frequency). For elastic moduli less than 1.5–2.5 MPa, the axial extent of the bubble cloud was found to depend weakly on the medium elasticity. In stiffer media, the incidental nucleus did not grow sufficiently large over the duration of the histotripsy pulse to initiate shock scattering. Consequently, bubble cloud formation was suppressed and cavitation activity was limited to the expansion of the initial nucleus.

The dependence of the maximum bubble size on the elastic model is shown in figure 6. For shock-scattering histotripsy excitations, predictions of bubble expansion based on finite-strain elastic models (classical, Kelvin–Voight, and neo-Hookean) agreed within 11.9% when the medium elastic modulus was less than 1 MPa. Predictions of the maximum bubble size for strain-hardening elasticity (modified Fung and Gent) and Kelvin–Voight elasticity were within 10% when the elastic modulus was less than 100 kPa. Beyond 100 kPa, bubble expansion was suppressed in strain-hardening media compared to finite-strain elasticity.

Zoom In Zoom Out Reset image size

For microtripsy excitations, the predicted maximum bubble size was within 18.8% for finite-strain elastic modulus when the elastic modulus was less than 1 MPa. The predicted maximum bubble size was on average 62% smaller for the strain-hardening models compared to the Kelvin–Voight model when the elastic modulus was less than 100 kPa. For all the elastic models and insonation schemes (frequency and peak negative pressure), bubble expansion was suppressed in a medium with an elastic modulus greater than 20 MPa.

Previous measurements of the maximum size of microtripsy induced cavitation as a function of frequency and elastic modulus dependence in agar phantoms (Vlaisavljevich et al 2015b) and predictions of the analytic model are shown in figure 7. Kelvin–Voight, modified Fung, modified Gent, or no elastic model (i.e. G  =  0) were considered in the analytic calculations (table 1). A comparison of the measured bubble size and predictions of the analytic model with Kelvin–Voight elasticity agreed within 9.4  ±  6.3 µm (3.8%  ±  1.4%) when the elastic modulus was 1.13 kPa. With increasing elastic modulus, the predictions of the model with Kelvin–Voight elasticity and the measurements diverged, with a 123.9  ±  95.7 µm (236.6%  ±  66.5%) difference at 570 kPa. Over the same range of elasticities, the analytic model prediction agreed with the measured bubble diameter within 38.8  ±  52.4 µm (14.3%  ±  13.0%) for modified Fung elasticity, and within 27.1  ±  36.6 µm (19.4%  ±  23.6%) for modified Gent elasticity. The modified Gent model appeared to be accurate within the standard deviation of the experimental measurements for frequencies 345 kHz and 1.5 MHz. The modified Fung model under predicted the maximum bubble size for frequencies 345 kHz and 1.5 MHz, but was within experimental error for 500 kHz and 3 MHz.

Zoom In Zoom Out Reset image size

In addition to the elastic models in table 1, an elastic model optimized to the experimental data was also utilized in the analytic calculation. The optimized elastic model denoted that the elastic pressure contribution approached 30G for large bubble expansion, or ξ  =  ξ_HA  −  30G/P₀ and p_EM  =  ${{P}_{0}}\left(1+\frac{30G}{{{P}_{0}}} \right)$. Predictions of the maximum bubble size with the optimized elastic model were within experimental error for all conditions considered, except for 3 MHz insonation of agar phantoms with an elastic modulus of 21.7 kPa (figure 8). For all conditions, the predicted maximum bubble diameter based on the best fit elastic model was within 9.5  ±  11.0 µm (10.2%  ±  10.0%) of the measured maximum bubble diameter.

Zoom In Zoom Out Reset image size

The analytic prediction of the maximum bubble size appeared to be a well-behaved function of the medium elasticity and the histotripsy insonation parameters (figure 8). For a given set of insonation conditions (microtripsy versus shock-scattering histotripsy, fundamental frequency, and peak negative pressure), a critical elastic modulus was observed. For tissues stiffer than the critical elastic modulus, the maximum bubble diameter decreased exponentially with increasing medium elasticity. To reflect this critical elastic behavior, the maximum bubble size predicted by the analytic model was fit in the least squares sense to a function of the form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{D}_{MAX}}~=~{{a}_{1}}{{\left| {{p}_{N}} \right|}^{{{a}_{2}}}}f{{~}^{{{a}_{3}}}}{\rm exp}\left[ -(E-{{a}_{4}})/{{a}_{5}} \right]\nonumber \end{align} \tag{ 15 }$$

using the 'nlinfit' function in MATLAB, where the maximum bubble diameter in micrometers is denoted by D_MAX, |p_N| is the magnitude of the peak negative pressure in MPa, f is the fundamental frequency of the histotripsy pulse in MHz, E is the elastic modulus in kPa, and a₁, a₂, a₃, a₄, and a₅ are fitting parameters. The critical elastic modulus is modeled in (15) via the parameter a₄, and the exponential decay via parameter a₅. The fitting parameters are listed in table 2 in terms of the type of histotripsy excitation, along with mean square error of the fit.

The maximum bubble diameter was inversely proportional to frequency (a₃  =  −1 for both shock scattering and microtripsy) and directly proportional to the peak negative pressure (a₂  =  0.56 for shock scattering histotripsy and 0.93 for microtripsy). The parameters a₄ and a₅ indicate tissue elasticity has a stronger influence on microtripsy excitation of cavitation nuclei compared to shock-scattering histotripsy when the elastic modulus is less than 2.95 MPa. Beyond 2.95 MPa, shock scattering histotripsy-induced cavitation will be more strongly influenced by the tissue elasticity compared to microtripsy.

The viability of red blood cells is shown as a function of distance from the bubble wall for microtripsy excitations in figure 9, and for a single-cycle shock-scattering histotripsy excitation in figure 10. Strains were computed analytically at the maximum bubble diameter for Kelvin–Voight elasticity. The red blood cell viability increased rapidly from 0% at the bubble wall to 100% within 50 µm for all insonation conditions. The red blood cell viability increased sigmoidally as a function of distance away from the bubble wall, and was fit to a function of the form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle V\left(d \right)~=~\frac{100\%}{1+{\rm exp}\left[ -\frac{\left(d-{{d}_{0}} \right)}{{{d}_{s}}} \right]}\nonumber \end{align} \tag{ 16 }$$

where d is the distance from the center of the bubble, and d₀ and d_s are fitting parameters. The parameter d₀ was used to express the distance from the bubble wall over which the predicted red blood cell viability decreased by 50%:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{r}_{50}}~=~{{d}_{0}}-{{R}_{MAX}}.\nonumber \end{align} \tag{ 17 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The frequency, pressure, and elasticity dependence of r₅₀ had a similar form to that for the maximum bubble size, as indicated in figure 11. For a medium with elastic modulus greater than 20 MPa, r₅₀ approached 0 µm as bubble expansion was completely suppressed. For fixed medium elasticity and insonation conditions (peak negative pressure and fundamental frequency), r₅₀ was greater for microtripsy pulses than shock scattering histotripsy pulses (table 3). Thus, the distance over which bubble expansion effects red blood cell viability is greater for microtripsy compared to shock scattering histotripsy. For fundamental frequencies greater than 1 MHz, changes in red blood cell viability did not extend beyond 5 µm, roughly the size of a red blood cell (Lewis 1996).

Zoom In Zoom Out Reset image size

Changes in cell viability induced by strain from the bubble collapse are shown in figure 12. While strain during the expansion pushes red blood cells away from the bubble, strain during collapse pulls the red blood cells toward the bubble. For conditions of massive bubble expansion (i.e. medium elasticity less than 20 MPa), the cell viability at a fixed distance from the bubble varied by less than 0.1% over the range of media elasticity or insonation conditions considered in this study. Bubble collapse-induced changes in cell viability appear to be restricted within 2 nm for microtripsy excitations, and 10 nm for shock scattering histotripsy.

Zoom In Zoom Out Reset image size

An analytic model to predict histotripsy-induced bubble expansion (Bader and Holland 2016) was extended to include the effects of medium elasticity. Good agreement was observed between the predictions of the analytic model and numerical calculations based on integration of the Yang/Church model (Yang and Church 2005) for medium elasticities that span soft tissue (elastic modulus less than 1 MPa (Duck 1990)), figure 4. For a medium with an elastic modulus between 1.8 and 20 MPa, the analytic model underestimates the maximum bubble size compared to numerical computations for microtripsy excitations. For a medium with an elastic modulus less than 1.8 MPa, the analytic model and numerical computations agree within 1.5% for both forms of histotripsy considered in this study.

The maximum size of a nucleus in a viscoelastic medium exposed to a single-cycle excitation can be computed analytically via (3) in combination with table 1. Microtripsy insonations rely on a single tensile pulse, and bubble expansion is sufficiently described by the analytic model (3). For shock scattering histotripsy, bubble clouds are generated over the course of several cycles (Maxwell et al 2012). The onset of bubble cloud formation occurs once the initial bubble grows sufficiently large to scatter the incident shock wave geometrically (Maxwell et al 2011b, Bader and Holland 2016). Analytic predictions of bubble expansion over the duration of the shock scattering histotripsy pulse are in good agreement with numerical computations for mediums with an elastic modulus less than 875 kPa, as shown in figure 5. The analytic model overestimates bubble expansion in a medium with elastic modulus between 875 kPa and 2.5 MPa. For tissues over this range of stiffness, the axial extent of the bubble cloud, and therefore the extent of tissue liquefaction, will be overestimated for the analytic formulation.

Analytic predictions of maximum bubble size were also compared to experimental measurements of microtripsy-induced cavitation (figure 7). Analytic predictions based on strain-hardening elastic models (Fung and Gent) were the closest to measured bubble sizes in agar phantoms, consistent with the observation that agar is a strain-hardening medium (Hall et al 1997, Erkamp et al 2000). While there was some discrepancy between the measured values and the strain-hardening models, incorporation of a 'best-fit' elastic model improved the accuracy of the analytic predictions 8.7% and 3.3% compared to the use of the Gent and Fung models, respectively. The agreement between the analytic prediction and measured bubble sizes indicates the analytic model captures the behavior of histotripsy-induced cavitation, but requires the proper medium-dependent elastic model for accurate prediction of the bubble size (Schmitt et al 2011). Alternatively, there may be a strain rate dependence (and therefore frequency dependence) contribution to the bubble dynamics not properly modeled here.

Medium elasticity was found to influence the cavitation dynamics through two mechanisms: a retardation to the onset of bubble expansion, and a modification to the effective pressure that halts bubble expansion. The longest delays in bubble expansion are dictated by surface tension compare to other mechanisms for peak negative pressures less than 5 MPa (figure 3(B)). Histotripsy operates at peak negative pressures greater than 15 MPa (Maxwell et al 2012), where inertial and viscous effects are the major retarding effects to bubble expansion for soft tissues and muscle (figure 3(B)). As the elastic modulus exceeds 105 MPa, bubble expansion will be dictated by the medium elasticity (figure 3(C)). Thus, histotripsy-induced bubble dynamics in tissues such as hair, nail, collagen, and mineralized tissue (e.g. bone, tooth, kidney stones, etc) will be dictated by the medium elasticity (Duck 1990). This also has interesting implications for histotripsy studies with fixed tissues. Fixation of tissues in formalin has been shown to increase the elastic modulus by two orders of magnitude compared to fresh tissues (Braun et al 2017). Such stiffening effects of the tissue will cause additional delay to the onset of bubble expansion compared to fresh tissue, potentially altering the dominant mechanisms of bubble dynamics (Hoffmeister et al 1996).

The medium elasticity also influences the effective pressure, or 'Blake brake' (equations (7) and (8)). The form of the effective pressure depends on the insonation type (Bader and Holland 2016), and therefore the means by which elasticity effects the maximum bubble diameter. A critical elastic modulus was identified for both forms of histotripsy, a₄ in table 3. For a given set of insonation conditions, the maximum bubble size was found to decrease exponentially for an elastic modulus greater than the critical modulus. For shock-scattering histotripsy, which has a large compressional component in the excitation waveform, the best-fit for a₄ was found to be 2.95 MPa. The critical elastic modulus for microtripsy, which has no compressional component, was found to be reduced by a factor of seven compared to shock-scattering histotripsy (425 kPa). The best-fit values for a₄ indicate bubble expansion is more resistant to tissue elasticity for shock-scattering histotripsy compared to microtripsy over the elastic modulus range 425 kPa to 2.95 MPa.

The axial extent of the bubble cloud induced by shock-scattering histotripsy also has a critical elastic modulus, decreasing rapidly for a medium stiffness greater than 1.5–2.5 MPa (figure 5(B)). Beyond the critical elastic modulus, the initial cavitation nucleus does not expand sufficiently to scatter the incident shock wave geometrically. Consequently, the efficacy of shock-scattering histotripsy will be significantly reduced for a medium with elasticity greater than 2.5 MPa. The erosion of tissue via shock-scattering histotripsy has been demonstrated to be significantly reduced for tissues with elastic moduli ranging from several MPa to GPa (Duck 1990), including vascular adventitia, skin, cartilage, tendon, bone, and tooth (Vlaisavljevich et al 2013b).

The critical elastic moduli denoted in these studies suggest the insonation scheme should be chosen to induce a specific bubble dynamics (shock-induced bubble cloud or microtripsy bubble expansion) based on a priori knowledge of the tissue elasticity. Microtripsy can nucleate cavitation in a medium with elasticity up to 425 kPa, beyond which the efficacy of the therapy will be reduced. The formation of shock-scattering bubble clouds is impervious to the medium elasticity up to 2.95 MPa, beyond which the cloud expansion is suppressed and the efficacy of shock-scattering histotripsy is mitigated. The explosive bubble expansion necessary for microtripsy tissue liquefaction can still be nucleated in a medium with elasticity greater than 2.95 MPa (figure 4), but at a reduced efficacy.

For both forms of histotripsy, bubble expansion was completely suppressed for a medium with elastic modulus greater than 20 MPa (figures 9 and 10). The elastic modulus for arterial tissue can reach 20 MPa under systolic pressure (Nichol et al 2011), indicating histotripsy liquefaction would be suppressed. Histotripsy has been noted to act as a vessel sparing therapy (Vlaisavljevich et al 2013a, 2013b, 2017, Zhang et al 2017), consistent with this conclusion. Bubble expansion was also suppressed for the largest elastic moduli considered in this study, which are representative of mineralized tissues such as cortical bone and kidney stones (5–80 GPa), enamel (20–80 GPa), and dentine (11–21 GPa) (Duck 1990, Heimbach et al 2000). The largest tension considered in this study was based on the current state of histotripsy technology (Sukovich et al 2016). In order to generate cavitation in mineralized tissue (elastic modulus ~ 10 GPa), a peak negative pressure in excess of 1 TPa would be required based on analytic model predictions. For perspective, the side lobes of a highly focused source (f-number of 0.5, 6.5 cm focal length and 10 cm diameter, 1 MHz fundamental frequency) would be in excess of 2 GPa to generate 1 TPa in the focal zone (Mast et al 2005). Surface erosion thus appears the most feasible means for histotripsy treatment of mineralized tissue (Xu et al 2005, Duryea et al 2015), restricting the targets to superficial disease.

These studies focused primarily on the use of the Kelvin–Voight elastic model because of its applicability for soft tissues (Frizzell et al 1976, Madsen et al 1983). The predicted maximum bubble size was similar for all the finite-strain models (Kelvin–Voight, linear, and Neo-Hookean in figure 6), indicating that the large bubble expansion for soft tissues is only mildly dependent on the form of the elastic model. The best agreement between measured bubble sizes in an agar phantom and the analytic model were for strain-hardening elastic models (Fung and Gent models in figure 7). Agar displays behaviors similar to a strain hardening medium, with a stronger increase in elasticity with strains than ex vivo soft tissue (Erkamp et al 2000, Hall et al 1997) and fracturing at moderate strain (Madsen et al 2003). Strain-hardening models were developed to capture the response of tissue with high content of collagen fibers, such as skin and arterial walls (Tong and Fung 1976, Deng et al 1994, Horgan and Saccomandi 2003). Interestingly, this suggests that bubble behavior observed in agar phantoms may not be totally reflective of that occurring in soft tissues, but in tissue that derives its structure from extracellular components (i.e. collagen).

Red blood cell viability surrounding histotripsy cavitation activity was modeled based on measurements of impulsive stretching of erythrocytes (Li et al 2013). Strain fields surrounding the bubble were rapidly attenuated within 50 µm from the bubble wall, consistent with previous calculations of histotripsy-induced cavitation (Mancia et al 2017). Only red blood cells within a few micrometers from the bubble wall experienced lethal strain (table 3), and the erythrocyte viability increased sigmoidally with distance away from the bubble. On a macroscopic level, the limited distance over which red blood cell viability is reduced observed in this study is consistent with the sharp boundaries observed in histotripsy lesions (Parsons et al 2006, Kieran et al 2007, Vlaisavljevich et al 2013a, Darnell et al 2015). In vitro studies have demonstrated histotripsy-induced breast cancer cell bisection occur only within a few micrometers to a nucleated bubble (Vlaisavljevich et al 2016). Under the similar insonation conditions, the calculations in this study indicate changes in red blood cell viability would be restricted within 1 µm of the bubble (figures 11(C) and (E)).

Tissue elasticity was observed to influence the degree of bubble expansion, and the subsequent strain surrounding the bubble. Erythrocyte viability for a given histotripsy exposure will be restricted to the bubble wall in highly elastic tissue, reducing the volume of red blood cell lysis. This has implications for the age of histotripsy-treated thrombus. Retraction of the thrombus over time changes its elasticity (Xie et al 2005) and composition (Sutton et al 2013). Zhang et al (2016) noted a change in in elastic modulus from 2.31 kPa for unretracted clots to 11.12 kPa for retracted clots, with a significant reduction in the efficacy of 1 MHz microtripsy pulses. Interestingly, the calculations here suggest a 1.1% decrease in bubble expansion and 3.5% reduction in red blood cell lysis over the same range of elastic moduli. The change in thrombus structure may alter the acousto-mechanical properties, increasing the attenuation of the histotripsy field (Nahirnyak et al 2006). Alternatively, anisotropies in fibrin-rich clot may not be properly modeled by the normal strains computed in this study (Weiss et al 2013). Up to 43% of thrombi are multifaceted, with both acute and chronic portions (Mewissen et al 1999, Gagne et al 2015). To ensure sufficient mechanical action throughout the thrombus, histotripsy image guidance should thus focus on not only qualifying the presence of cavitation, but also spatially quantifying the degree of cavitation activity (Macoskey et al 2017, Bader et al 2018).

In media stiffer than 20 MPa, calculations indicate the strain on erythrocytes is insufficient to induce lysis for the histotripsy insonation parameters considered (figure 11). In softer media, the extent of damage beyond the bubble volume was found to depend on the insonation type. Microtripsy-nucleated bubbles induce larger strains at greater distances from the bubble compared to shock scattering cavitation for a given set of insonation parameters. Depending on the necessity for precision relative to the observed spatial location of the bubble cloud, the choice of insonation scheme should be considered.

There are several aspects of this study that limit the generalizability of these findings. The amplitude of the measured shock scattering histotripsy pulses used as the source term in (1) may be underestimated due to spatial averaging of the fiber width and the bandwidth of the system (Canney et al 2008). The size of the bubble is assumed to be much smaller than the wavelength of the excitation pressure for the Yang/Church model, which is not the case for shockwave excitation. The analytic model does not consider the effects of bubble-bubble interactions, which are likely for a highly dense bubble cloud. Tissue damage was modeled here due purely to strain in the medium surrounding the bubble (Li et al 2013), but the type of cell death can be dictated by the strain rate (Bar-Kochba et al 2016). The strain-dependent cell viability was based on impulsive measurements of erythrocyte stretching (Li et al 2013), and may not be applicable across all cell types (Barbee 2005, Ofek et al 2009). The use of (17) to predict the minimum bubble radius is based on the collapse of a Rayleigh cavity in an elastic medium, and will underestimate the actual minimum bubble size due to damping (Gaudron et al 2015) and gas diffusion (Church 1989). The goal of this study was to estimate the strain during collapse, which depends on the difference in the cube of the minimum bubble size and the size of the initial bubble nucleus via (16). For large bubble expansion, such as histotripsy-induced cavitation (Bader and Holland 2016, Mancia et al 2017), the minimum bubble size will be much smaller than the size of the nucleus, (13) will be largely independent of the minimum bubble size. Normal strains were computed in this study, which would be indicative of isotropic tissues. The strain will have an angular dependence in anisotropic tissues such as tendon and myocardium (Hoffmeister et al 1996). The analytic model does not include the contribution of gas diffusion in the prediction of the maximum bubble size. Bubble expansion is controlled by the inertia of the medium, with minimal contributions from diffusion (Church 1989).