Effects of frequency on bubble-cloud behavior and ablation efficiency in intrinsic threshold histotripsy

A custom-designed, modular histotripsy transducer array with an aperture size of 120.5 mm, a geometric focus of 75 mm, and an f-number of 0.62 was used for all experiments in this study (figure 1). The transducer was comprised of three concentric rings of six, twelve, and fourteen ports, which could be populated by thirty-two interchangeable 500 kHz, 1 MHz, or 3 MHz elements.

Zoom In Zoom Out Reset image size

For all experiments, the transducer was positioned horizontally in a tank of degassed tap water with a dissolved oxygen concentration <28% of saturation. The individual elements were driven by a custom high-voltage amplifier designed to generate short therapy pulses of <2 cycles. The amplifier was populated with interchangeable amplifier boards to match the frequency of the elements being driven. A field-programmable gate array (FPGA) board (Altera DE0-Nano Terasic Technology, Hsinchu, Taiwan) was programmed with a custom microcontroller to synchronize the pulse arrivals by controlling the phase-delay and amplitude of each individual element. A computerized positioning system was used to guide the motion of agarose tissue phantoms or hydrophones used in the respective experiments. Microcontroller instructions and the positioner locations were serialized and communicated via a MATLAB (The MathWorks, Natick, MA, USA) script for each experiment (figure 1).

A fiber optic probe hydrophone (FOPH) measured pressure waveforms for 500 kHz, 1 MHz, and 3 MHz pulsing conditions (figure 2). The FOPH was cross-calibrated with a high-sensitivity reference rod hydrophone (HNR-0500, Onda Corporation, Sunnyvale, CA, USA) at nine points between 0 and 5 MPa peak-to-peak to ensure the FOPH accurately measured the pressures. The FOPH directly measured the waveforms in degassed water up to p− = 18.5 MPa. At higher pressures up to p− = 42 MPa, the focal pressures were estimated by summing the outputs of transducer sub-apertures of sixteen elements or eight elements to avoid cavitation damaging the hydrophone (Vlaisavljevich et al 2017). An oscilloscope (TBS2000 series, Tektronix, Beaverton, OR, USA) was used to process all waveforms at a sample rate of 500 MS s⁻¹. The waveforms were averaged over 512 pulses to minimize noise from the FOPH, and the final waveform was recorded to MATLAB. Finally, in addition to the focal waveforms, 1D beam profiles in the axial, transverse, and elevational directions were measured using the high-sensitivity rod hydrophone for all frequency conditions.

Zoom In Zoom Out Reset image size

The viscoelastic properties of soft tissue were replicated in this study by 1% (w/v) agarose tissue phantoms using methods similar to those in prior cavitation studies (Maxwell et al 2010b, Vlaisavljevich et al 2017, Edsall et al 2021). A prior study showed 1% w/v agarose to have a Young's modulus of 21.7 kPa (Vlaisavljevich et al 2015c), which is within the range of ∼6–25 kPa characteristic of relevant target tissues for histotripsy treatments (Yamada 1973, Duck 1990). The agarose gel tissue phantom was made by mixing 0.5% agarose powder (Type VII-A, Sigma Aldrich Corporation, St. Louis, MO, USA) with 99.5% deionized water at room temperature until evenly dispersed. The mixture was heated until boiling in a microwave and then stirred until the agarose powder was fully dissolved. The mixture was repeatedly heated to boiling and immediately stirred to induce flash boiling to release dissolved gas from the mixture. This process was repeated until 50% of the volume remained to obtain a degassed 1% (w/v) liquid agarose gel mixture. The mixture was placed in a vacuum desiccator and kept under a partial vacuum (∼33.62 kPa, absolute) for twenty-five minutes to remove remaining gas and minimize re-gassing as the agarose cooled. Any remaining gas was forced from solution by decreasing the pressure to 16.75 kPa for the final five minutes. Once the temperature of the agarose mixture dropped to 40 °C, 100 ml of the gel was slowly added by serological pipette down the wall of a rectangular silicone mold into a custom-designed poly-lactic-acid (PLA) phantom holder frame. The agarose-filled mold was stored in a refrigerator for one hour to solidify. Experiments were performed within two hours of each gel's creation to ensure a consistent agarose concentration and maintain the gel's degassed state during testing.

Optical imaging was used to characterize the bubble clouds generated by the 500 kHz, 1 MHz, and 3 MHz pulsing. A machine-vision camera (FLIR Blackfly S monochrome, BFS-U3-32S4M-C 3.2 MP, 118 FPS, Sony IMX252, Mono, FLIR Integrated Imaging Solutions, Richmond, BC, Canada) with a global shutter and a 100 mm F2.8 Macro lens (Tokina AT-X Pro, Kenko Tokina Co., LTD, Tokyo, Japan) captured images with a resolution of 3.25 μm per pixel and a depth of field of ∼0.7 mm. The samples were backlit using a custom high-speed LED strobe light. Both the camera and strobe light were individually controlled by triggers from the amplifier box with the camera shutter opening synchronized with the histotripsy pulse generation and the strobe acting as the shutter. Strobe duration was kept as short as possible (3 μs) to minimize the motion blur of the expanding bubbles while maintaining sufficient lighting. All exposures were centered to a delay at 7 μs after pulse convergence. This time point occurred prior to the time between 8 and 50 μs at which the bubbles have previously been reported to reach their maximum radius, R_max, at these tested frequencies (Vlaisavljevich et al 2015c, Mancia et al 2017) and was determined to be the optimum time delay to detect clearly visible bubbles while minimizing the overlap of adjacent bubbles within the cloud as became more prevalent with larger bubbles at later time points. The timing accuracy of the amplifier, camera trigger, and strobe trigger was accurate to the 10 ns clock speed of the FPGA board. A 0.5 Hz pulse repetition frequency (PRF) was used for bubble-cloud characterization experiments to minimize memory effect on subsequent pulses (Maxwell et al 2013). Since each image was taken very early in the life of the bubble cloud at a low PRF, it was possible to accurately identify individual bubbles within each image for most samples. The optical images for each sample were analyzed using a custom MATLAB script that identified bubble clouds for dimensional measurement by converting the captured grayscale image (figure 3(B)) into a negative binary image based on an intensity threshold determined by the difference from the background intensity as described in a previous study (figure 3(C)) (Maxwell et al 2013). For counting the number of bubbles in the cloud, individual bubbles were identified using this same thresholding technique with the additional use of distance transform watershed lines to identify curvature to separate overlapping bubbles. Only regions with >5 pixels were counted as bubbles and the output of each image was visually confirmed by at least two independent readers before analysis. Using these criteria, the minimum resolvable bubble diameter was ∼16 μm.

Zoom In Zoom Out Reset image size

Cavitation thresholds were obtained for each frequency condition by applying one-hundred histotripsy pulses to a single point in the tissue phantoms at p− levels ranging from 0 to 40 MPa in increments of 2 MPa. The possibility of an earlier pulse altering the probability of cavitation on the subsequent pulse was minimized by applying the pulses at a PRF of 0.5 Hz to allow sufficient time between pulses or the residual bubbles to dissolve before the subsequent pulse arrived (Maxwell et al 2013, Vlaisavljevich et al 2015a). All points were created at a uniform depth into the agarose phantom in the axial direction. A new point was used for every pressure separated by at least two times the predicted transverse/elevational dimensions of the cloud as determined by 1D beam profiles scaled to 40 MPa. The cavitation probability, p_cav, was determined for each pressure level as the fraction of the total one-hundred pulses for which cavitation was detected by high-speed optical imaging. The cavitation threshold, p_t , is defined as the p− corresponding to p_cav = 0.5.

The cavitation threshold for each frequency, p− _500kHz , p− _1MHz , p− _3MHz , and p− _Dual , was determined using a previously reported method (Maxwell et al 2013, Vlaisavljevich et al 2015b) whereby the probability of observing cavitation follow a sigmoid function given as

$$\begin{eqnarray}&&P\left({p}_{-}\right)=\displaystyle \frac{1}{2}\left[1+erf\left(\displaystyle \frac{{p}_{-}\,-\,{p}_{t}}{\sqrt{2{\sigma }^{2}}}\right)\right],\end{eqnarray} \tag{ 1 }$$

where erf is the error function, p_t is the negative pressure at which the probability p_cav = 0.5, and σ is a variable related to the width of the transition between p_cav = 0 and p_cav = 1, with ± σ giving the difference in pressure from about p_cav = 0.15 to p_cav = 0.85 for the fit (Maxwell et al 2013). MATLAB software was used for curve fitting all data sets.

The predicted bubble-cloud dimensions generated by each frequency and pressure level were estimated by following a previously published method (Roberts et al 2014, Vlaisavljevich et al 2017). 1D beam profiles measured by the high-sensitivity rod hydrophone in the elevational, transverse, and axial directions of the transducer for all frequency conditions were normalized to the p− at the focus for each pressure level tested experimentally. The predicted bubble-cloud dimensions for each p− were then estimated by measuring the distance between the intersections of the experimentally measured cavitation thresholds and the normalized beam profile (figure 3(A)). These predicted bubble-cloud dimensions were compared to the cloud dimensions measured experimentally using optical imaging, as described in the next section (figure 3(D)).

The effect of frequency on the bubble-cloud size was observed by capturing bubble clouds for each frequency condition of 500 kHz, 1 MHz, and 3 MHz using optical imaging, as described above. For each frequency, 75 histotripsy pulses were applied to agarose phantoms for each frequency and pressure level combination. To minimize the effects of tissue phantom liquefaction on the results, 25 histotripsy pulses each were applied to three separate locations within the agarose tissue phantom, with the results pooled together to total 75 bubble-cloud images for each pressure level tested. The cloud length was measured as the distance from the closest edge of the nearest bubble to the transducer to the furthest edge of the furthest bubble from the transducer. Similarly, the cloud width was defined as the distance from the furthest edge of the bubble highest in the positive transverse direction to the furthest edge of the lowest bubble in the negative transverse direction (figure 3(D)).

A custom MATLAB cloud analysis program counted the number of bubbles in each binary image by identifying the centroids to determine the number of bubbles in each cloud. The area of the cloud was determined by an elliptical area approximation using the axial and transverse 1D beam profile measurements scaled to the respective pressure levels and the equation: ${\rm{A}}{\rm{r}}{\rm{e}}{\rm{a}}=\left(\tfrac{{\rm{A}}{\rm{x}}{\rm{i}}{\rm{a}}{\rm{l}}\,{\rm{D}}{\rm{i}}{\rm{a}}{\rm{m}}{\rm{e}}{\rm{t}}{\rm{e}}{\rm{r}}}{2}\right)\left(\tfrac{{\rm{T}}{\rm{r}}{\rm{a}}{\rm{n}}{\rm{s}}{\rm{v}}{\rm{e}}{\rm{r}}{\rm{s}}{\rm{e}}\,{\rm{D}}{\rm{i}}{\rm{a}}{\rm{m}}{\rm{e}}{\rm{t}}{\rm{e}}{\rm{r}}}{2}\right)\pi .$ The bubble density within the cloud was calculated by dividing the number of bubbles for each cloud by the corresponding predicted cross-sectional area of the cloud (figure 3(E)). Using the imaging techniques described in the cavitation detection using optical imaging section above, the maximum number of bubbles that could be theoretically counted at the minimum bubble diameter would be ∼3800 bubbles mm⁻² assuming the bubbles were aligned diameter-to-diameter in a uniform matrix across the focus. This resolution decreases as bubbles expand larger than the minimum resolvable bubble diameter. It should be noted that due to these limitations, in bubble cloud regions of high bubble overlap and coalescence, this method may not completely identify the entirety of every individual bubble but will still identify the centroid of the majority of the bubbles in these regions. The differences in mean cloud density across all pressure levels at each frequency and across all frequency conditions were compared in PRISM (GraphPad Software, San Francisco, CA, USA) using a two-way ANOVA analysis and Student's t-test comparisons. For all comparisons, p-values less than 0.05 were considered significant.

RBC phantoms comprised of three agarose layers, with the middle layer containing 5% (v/v) red blood cells were created for ablation experiments (Maxwell et al 2010b, Vlaisavljevich et al 2013a, Edsall et al 2021). A 1% agarose mixture was produced using the preparation described above but using degassed 0.9% saline in place of deionized water. Fresh porcine blood was obtained from subjects in an unrelated study and added to an anticoagulant solution of Citrate Phosphate Dextrose Anticoagulant (CPD, Sigma Aldrich Corporation, St. Louis, MO, USA) with a CPD-to-blood ratio of 1:9 ml. Whole blood was separated by centrifugation at 1250 RCF for 10 min. The plasma and white buffy coat were removed by transfer pipet, and the RBCs were saved for addition to the phantom. An initial layer of agarose was poured into the silicone mold at 45 °C to create the base of the RBC phantom. The mold was placed in a refrigerator at 4 °C to allow the base agarose layer to cool and solidify. The remaining agarose mixture was kept in the vacuum chamber until it had cooled to 38 °C. Once cooled, 9.5 ml of the agarose was combined with 0.5 ml of RBCs by gentle inversion. The liquid RBC-agarose solution was poured onto the set agarose layer and quickly tilted to coat the entire surface. The excess was then immediately discarded leaving behind a thin layer of the RBC-agarose solution. The phantom was replaced in the refrigerator for five minutes. Once the RBC-agarose layer was solidified, the remaining agarose without RBCs was poured on top of the first two layers to fill the silicone mold. The phantom was replaced in the refrigerator for one hour to fully solidify, producing a thin layer of RBCs, approximately a single cellular layer thick, suspended in the center of the agarose phantom.

Tissue ablation for histotripsy bubble clouds was characterized using RBC phantoms. Ablation of the RBCs is visualized through successive pulses lysing the cells turning a region of the RBC layer from opaque to translucent (Maxwell et al 2010b, Wang et al 2012, Roberts et al 2014, Vlaisavljevich et al 2017). The RBC phantom was oriented at the focus of the transducer such that sound propagated parallel to the RBC layer, and 2000 histotripsy pulses were applied to the RBC layer at a PRF of 1 Hz and p− = 38 MPa for all frequency conditions. The layer was aligned with the center of the cloud using 0.1 mm increments on the positioning system to ensure the center of the cloud aligned with the RBC layer. The bubble cloud from each pulse and the resulting ablation after each pulse were recorded by high-speed optical imaging. A custom MATLAB script was used to separate the cloud images and the ablation images, and measure the ablated area for each pulse using image thresholding. Finally, ablation efficiency was determined as the number of pulses needed to completely ablate the central cross-sectional plane of the focus, which was determined by the mean area of the bubble cloud at the applied p−. The ablated area of each image was normalized by dividing each frame's measured ablated area by the average area of the bubble cloud to allow for a comparison of the relative ablative efficiency of histotripsy across all frequency conditions.

Cavitation threshold experiments showed cavitation bubbles were observed for all frequencies when the pressure was increased above ∼24–26 MPa, with the number of bubbles and the area covered by the cavitation clouds increasing as p− was increased (figure 4). For all conditions, well-defined and sharply demarcated bubble clouds were observed, with larger bubble clouds containing more bubbles seen at higher pressure levels, confirming the results suggested by previous intrinsic threshold histotripsy studies (Maxwell et al 2013, Roberts et al 2014, Vlaisavljevich et al 2015c, 2017). The intrinsic thresholds were measured to be p− _{500kHz =} 26.4 MPa, p− _{1MHz =} 28.7 MPa, and p− _{3MHz =} 24.8 MPa for 500 kHz, 1 MHz, and 3 MHz, respectively. A distinct threshold behavior was observed as expected for all experimental conditions with σ _mean values of σ _500kHz = 0.71, σ _1MHz  = 1.02, and σ _3MHz = 0.65_. Results show all frequencies produced histotripsy cavitation with distinct thresholds at similar pressures_.

Zoom In Zoom Out Reset image size

The mean individual bubble size for each frequency condition was also quantified at pressures directly above the respective cavitation thresholds with values of R_500kHz = 146.0 ± 16.7 μm, R_1MHz = 70.2 ± 7.7 μm, and R_3MHz = 29.9 ± 4.5 μm (n = 20). These measurements were obtained from images taken at the same time point (7 μs after the initial convergence of the separate pulses at the focus). Therefore, these measurements are predictably smaller than the values previously seen in modeling and in experimental studies where R_max was the primary variable (Vlaisavljevich et al 2015c, Bader and Holland 2016). However, although these measurements are not the maximum bubble radius the measured bubble sizes at this time point during the growth phase of bubble expansion are consistent with those reported in prior experimental and theoretical studies of bubble expansion in intrinsic threshold histotripsy (Vlaisavljevich et al 2015c, Mancia et al 2017). Supplementary figure 1 (available online at stacks.iop.org/PMB/66/225009/mmedia) shows magnified images of the individual bubbles measured for each frequency condition.

Images of histotripsy bubble clouds generated at p− ranging from 26 to 40 MPa (figure 4) showed larger bubble clouds were generated at higher p−. Close agreement between the measured and predicted cloud dimensions was observed for all frequency and pressure conditions (figure 5). The measured axial cloud dimensions shown in figure 5 ranged from 0.13 ± 0.32 mm to 5.78 ± 0.27 mm for 500 kHz, 0.23 ± 0.05 mm to 3.34 ± 0.22 mm for 1 MHz, and 0.14 ± 0.22 mm to 1.48 ± 0.06 mm for 3 MHz. The measured transverse cloud dimensions ranged from 0.10 ± 0.21 mm to 2.05 ± 0.10 mm for 500 kHz, 0.20 ± 0.33 mm to 1.50 ± 0.25 mm for 1 MHz, and 0.11 ± 0.16 mm to 0.53 ± 0.02 mm for 3 MHz (figure 5).

Zoom In Zoom Out Reset image size

Optical images (figures 4 and 6) of bubble clouds formed at all frequency and pressure conditions were further analyzed to quantify the number of bubbles in the cloud and to determine the density of the bubbles within the bubble cloud. Results show the number of bubbles significantly increased with increasing p− for all frequencies (figure 7). More specifically, the mean number of bubbles measured within the bubble clouds for between p− of 30–40 MPa ranged from 21.7 ± 5.8 to 191.8 ± 34.5 bubbles at 500 kHz, 4.4 ± 2.1 to 131.8 ± 23.2 bubbles at 1 MHz, and 37.3 ± 5.7 to 135.1 ± 12.7 bubbles at 3 MHz (figure 7(A)). It is worth noting that the thresholding, bubble curvature, and centroid analysis may not account for all of the bubbles in regions with significant overlap of several bubbles. This was especially limiting at higher p− when the density at the center of the bubble cloud became increasingly dense with expanding bubbles, likely resulting in an underestimate of bubbles (and density) at higher p− values. An example of this limitation is provided in supplementary figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bubble-cloud density at each of the frequency conditions was evaluated for p− levels ranging from 30 to 40 MPa by dividing the number of bubbles in each cloud by the area of the predicted focal region at that frequency and pressure condition (figure 7(B)). This p− range was chosen to compare the effects of frequency on bubble-cloud density for conditions at which fully-formed bubble clouds were generated for each frequency condition. Only pressures for which complete clouds were formed across all three frequency conditions were considered. At lower p− values closer to the intrinsic threshold (24–28 MPa), the cavitation clouds contained only a few bubbles inside smaller focal regions, resulting in less consistent cloud density measurements. Comparing the results between frequency groups showed increasing cloud densities with increasing frequency (figure 7(B)). Results showed significant differences in the cloud density between all frequency conditions at each examined pressure level (p < .05), with 3 MHz generating the highest densities and 500 kHz the lowest densities at each examined p−. 1 MHz was intermediate to 3 MHz and 500 kHz conditions (figure 7(B)).

Histotripsy ablation experiments showed that the extent of the damaged area in the RBC phantoms increased with increasing pulse number across all conditions (figure 8). Complete ablation of the entire focal region was seen within 2000 pulses for all frequency conditions (figure 8). For all frequencies, ablation created well-defined lesions with the ablated area of the RBC layer localized to the region of the focus in which cavitation was observed. Consistent with prior studies (Maxwell et al 2010b, Lin et al 2014b, Vlaisavljevich et al 2017), lesions formed with sharply delineated boundaries between the treated and untreated regions of the RBC layer. The total ablated area increased with each consecutive pulse for all conditions, with the total ablated area decreasing with increasing frequency (figure 9(A)). For all tested frequencies, the final ablated area of the RBC layer closely matched the bubble-cloud area for the respective frequency (figures 8 and 9(A)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ablation efficiencies for each condition were quantitatively compared by plotting the ablation area normalized to the mean measured histotripsy bubble cloud formed in the RBC phantom as a function of pulse number (figure 9(B)). The normalized lesion areas for all frequency conditions increased with pulse number, but showed a more rapid increase at the lower frequencies, corresponding to an increased ablation efficiency (figure 9(B) and table 1). More specifically, the 500 kHz treatments had the most rapid ablation of all conditions, achieving >75% ablation after a mean of 74.7 ± 61 pulses (table 1). The 1 MHz was the next most efficient, removing >75% of the focus in 200.7 ± 115.9 pulses (table 1). Finally, the 3 MHz treatments proved to be the least efficient at ablating the RBC phantoms, requiring 325 ± 72.7 pulses to remove >75% of the focal volume (table 1). Figure 9(C) shows a magnified view of the RBC ablation efficiency for the first 100 pulses of the respective treatments to better visualize the differences in ablation efficiency during the initial pulses.