The relationship between bubble concentration and the acoustic emission energy of separate frequency bands

In microbubble-mediated ultrasound therapies, cavitation—volumetric oscillations of a bubble—is responsible for producing both therapeutic and adverse biological effects. Cavitation imposes stress on the surrounding tissue, either by direct bubble-tissue contact (Church and Miller, 2016) or by secondary effects, such as microstreaming, microjetting, or shockwave emission (Chen et al., 2011; Song et al., 2016; Chowdhury et al., 2020). The concentration of bubbles experiencing cavitation determines how often these mechanical forces are applied. Studies using similar sonication protocols but with different bubble concentrations showed significantly different outcomes (Ward et al., 2000). Currently, knowledge of the concentration of acoustic cavitation events in the body is limited.

Acoustic-based methods can estimate the concentration of intravenously injected preformed microbubbles using acoustic pressures and frequencies relevant to imaging procedures (Marsh et al., 1998; Mari et al., 2007; Lampaskis and Averkiou, 2010; Sciallero et al., 2011; Leithem et al., 2012). However, an acoustic-based method that uses therapeutically relevant pressures and frequencies is not available. In such a regime, bubbles are expected to act very differently, such as fracturing, dissolution, and undergoing inertial cavitation (Chomas et al., 2001; Lindsey et al., 2015; Shi et al., 2000; Ilovitsh et al., 2018). Hence, the concentration of cavitation events during a treatment is often unknown.

In this letter, we describe how, under a therapeutically relevant ultrasound field, the concentration of microbubbles affects the spectra of the captured bubble emissions. This is based on the manner in which emissions from many different bubbles interfere at a sensor. For each bubble inside the focal region, sound will travel at the speed of sound (⁠$c$⁠) from the emitter (located at the location $x→t$⁠) to the bubble (at $x→b$⁠), and from the bubble to the sensor (at $x→s$⁠). The time required for sound to travel this distance (⁠$Tsig$⁠) can be calculated,

Since the travel time differs for each bubble, the emissions from different bubbles arrive at the sensor with time differences (Sujarittam and Choi, 2020). For each transmitter-bubble-sensor arrangement, the value of $Tsig$ is bounded by an upper and lower limit, defined by the locations of the transmitter, the sensor, and the region where the bubbles are being sonicated. Thus, the time differences between emissions from different bubbles are bounded by a maximum value, $ΔTsigmax$⁠. The time difference between a given pair of bubbles is random and lower than this maximum time difference,

This maximum time difference influences the signal the sensor captures. Microbubbles emit a wide range of harmonics, subharmonics, ultraharmonics, and broadband frequencies—each having different periods. Low frequencies whose periods are much longer than the maximum time difference will always interfere constructively at the sensor. In contrast, frequencies with shorter periods (i.e., higher frequencies) could interfere either constructively or destructively, depending on their periods compared to the time difference between any given pair of bubbles. Thus, as the number of bubbles undergoing cavitation increases, the amplitudes of the lower frequencies, which always interfere constructively, will increase at a faster rate than those of the higher frequencies. Thus, the ratio between low frequency energy and high frequency energy may correlate with the concentration of bubbles undergoing cavitation.

The purpose of this study was to evaluate whether the energy ratio of low to high frequency content of acoustic emissions from bubbles changed with bubble concentration. Since the maximum time difference (⁠$ΔTsigmax$⁠) can be adjusted by changing the angle of the sensor relative to the transmission pulse and microbubbles [Figs. 1(a) and 1(b)], the applicability of this method depends on the relative position of the sensor. Two experimental conditions were investigated. In the first condition [Fig. 1(a)], the sensor was placed at a position that minimised $ΔTsigmax$⁠. In this figure, the purple line represents the longest transmitter-bubble-sensor distance, while the green line shows the shortest counterpart. The minimum $ΔTsigmax$ is the difference between the lengths of the two lines divided by the speed of sound. In the second condition [Fig. 1(b)], the sensor was placed where a longer $ΔTsigmax$ would be observed. In this configuration, the longest transmitter-bubble-sensor distance is shown by the purple line, while the shortest counterpart is shown by the red line.

An experiment was setup to capture signals where the emissions from different bubbles would arrive at a sensor with small time differences and at another with larger time differences [Fig. 1(c)]. For transmission, a 0.5-MHz focused transducer (focal width FWHM: 4 mm, focal length FWHM: 37 mm; focal distance: 62.6 mm; model H107; Sonic Concepts, Bothell, Washington) was driven by a function generator (model 33500B; Keysight, Santa Rosa, CA), after amplification by a 50-dB RF amplifier (model 2100 L; Electronics & Innovation, Rochester, NY). For reception, two 1-mm needle hydrophones (Precision Acoustics, Dorset, UK) with similar sensitivities and frequency response curves were placed 40 mm from the focal centre of the transducer, one at a 90º angle from the axis of the transducer, and the other concentrically aligned with the transducer. Its signal was amplified by the manufacturer's preamplifier and was digitised simultaneously at a sampling frequency of 100 MS/s (GaGe model Octave Express, DynamicSignals, Lockport, IL). The experiment was conducted in a tank of degassed, deionised, and filtered water.

As a target, a wall-less polyacrylamide gel channel was placed making a 45º angle from the axes of the transducer and of the sensor. The gel was casted around a 0.56-mm nylon wire, which was removed after the gel had set to form the channel. We used a cross-linker-monomer ratio that yielded the Young's modulus of approximately 8.73 ± 0.79 kPa. The gel composition and manufacturing protocol has been described by Tse and Engler (2010).

For microbubbles, dipalmitoylphosphatidylcholine (DPPC)-coated perfluorobutane microbubbles manufactured in-house were used. The manufacturing procedure was described elsewhere (Pouliopoulos et al., 2014). The bubble concentration and size distribution were measured using an optical microscope and an image processing software (Sennoga et al., 2010). In brief, A diluted solution of microbubble was pipetted onto a haemocytometer, covered with a clear slide, and placed under a microscope. A time of approximately 3 min was allowed for the bubbles to float into the focal plane, and an image was taken. A computer algorithm detected circular shapes of the microbubbles and recorded their diameters. The concentration of the undiluted microbubbles was (6.48 ± 2.46)×10⁹ MBs/ml. Less than 1% of the bubbles had diameters larger than 10 μm. The mean and standard deviation of the bubble diameters were 1.80 ± 1.96 μm (1.68 ± 1.40 μm if excluding those with diameters>10 μm). The bubble resuspension was diluted in deionised, filtered water, and drawn through the channel using a syringe pump (Harvard Apparatus, Holliston, MA) at a rate of 0.1 ml/min (approximate average velocity of 6.5 mm/s in the channel).

The transducer was excited using a 0.5-MHz sinusoidal pulse, creating an emitted waveform similar to a nine-cycle sinusoid. The peak-rarefactional pressures ranged from 0.46 to 0.92 MPa. A two-second pause was allowed between subsequent transmissions to let new bubbles replenish the channel. Three dilutions were used. The undiluted bubbles were diluted to factors of 1:10k, 1:100k, and 1:1000k microbubble-to-water; to nominal concentrations of (6.48 ± 2.46)×10⁵, (6.48 ± 2.46)×10⁴, and (6.48 ± 2.46)×10³ MBs/ml, respectively. Thirty-five sets of data were collected for each pressure and dilution. Control data were collected in the same manner but with bubble-free, pure water in the channel.

Background-subtracted signals were calculated and inspected as they provided insights into how emissions from different bubbles interfered. Signals captured from ultrasound-stimulated microbubbles were subtracted by an average of 10 control signals obtained under the same insonation pressures. This removed the background acoustic signal, leaving only the signal generated by the bubble emissions (Leighton et al., 2002).

We predicted that as the concentration of microbubble increased, the acoustic emissions would contain a higher proportion of low-frequency content relative to high-frequency content. We calculated the maximum time difference (⁠$ΔTsigmax$⁠) from the geometry of our experimental setup (Fig. 1). For the perpendicular sensor, the longest transmitter-bubble-sensor distance was from the bubble furthest from the transmitter; the shortest counterpart was from the bubble approximately at the proximal-centre of the channel. For the concentric sensors, the longest transmitter-bubble-sensor distance was from the bubble furthest from the transmitter; the shortest counterpart was from the bubble closest from the transmitter. The differences between the longest and shortest transmitter-bubble-sensor distances were 0.83 and 7.2 mm, for the two sensors, resulting in the $ΔTsigmax$ values of 0.55 and 4.8 μs, respectively (assuming c = 1498 m/s). This estimation assumed that the insonation pulse was approximately plane wave in the focus, that the channel was 0.56 mm in diameter, that the bubbles were distributed in the whole volume of the channel within the focused ultrasound beam, and that the receiver's position was a point source located at the center of the sensor surface. The size of the ultrasound beam for approximated by its −6 dB FWHM region. Thus, for the perpendicular sensor, as the bubble concentration increased, we expected that low frequency signals [lower than ((0.55 μs)⁻¹)/2 = 0.91 MHz] would increase at a faster rate than higher frequency signals. This trend should not be observed in the concentric sensor, where the low frequencies are below ((4.8 μs)⁻¹)/2 = 0.10 MHz. Inspecting the spectra of the perpendicular sensor's signals, it was found that the relative amplitudes of frequencies up to the third harmonic (1.5 MHz) increased relative to higher frequencies with the bubble concentration.

We applied Tukey windows (0.125 MHz stop-bandwidth; 0.75 rectangular window width) to the raw bubble signals (before background subtraction) in the frequency domain to isolate their harmonics. We categorized low-frequency signals as the 2nd and 3rd harmonics (1.0 and 1.5 MHz centre frequencies) and high-frequency signals as the 7th to 12th harmonics (3.5, 4.0,…, 6.0 MHz). The energy contained in each harmonic was quantified by integrating the windowed spectra in the frequency domain. The harmonic energy ratio (HER) was then calculated as the ratio of low-frequency harmonic energy to high-frequency harmonic energy. Despite the HER-concentration exhibiting either quadratic or linear relationships (see Sec. 3), we applied linear regressions to all data for the sake of simplicity and uniformity, rather than selectively applying different regression types to different data sets. R² values were calculated to assess the goodness of fit of the linear regression for each sonication pressure. We neglected the effects of directivity of the hydrophone because they were expected to be minimal, as the hydrophone was facing the bubble population, had a small diameter, and was placed several centimeters away from the bubbles. We calculated the directivity of the sensor by modelling the needle hydrophone as a circular piston in a rigid planar baffle (Morse and Ingard, 1986); it was found that the sensitivity to signals from bubbles within the focal region would not drop below 93% for all relevant frequencies (up to 10 MHz).

The background-subtracted time-domain waveform captured with the perpendicular sensor qualitatively changed with the bubble concentration [Figs. 2(a)–2(c)]. At the lowest bubble concentration tested, both low-frequency oscillations and sharp compressional impulses were observed [Fig. 2(a)]. The low-frequency oscillation cycle had a period equal to the excitation frequency. Each cycle contained multiple sharp compressional peaks but only one rarefactional peak. In the frequency domain, the second and third harmonics were prominent. Higher order harmonics and broadband content could also be observed at low magnitudes, especially with acoustic pressures of 0.55 MPa and below. At higher pressures, more broadband components were observed. As the bubble concentration increased, the low-frequency oscillation became more prominent [Fig. 2(c)]. The compressional peaks were less apparent, appearing smoother. At the highest concentration, a prominent compressional peak could be observed in each cycle. Most cycles contained only one secondary compressional peak with a very small amplitude. Overall, this resulted in a periodic waveform where two compressional peaks and one rarefactional peak appear on every acoustic period. In the frequency domain, the second and third harmonics were significantly higher than other frequencies. Higher order harmonics and broadband components were relatively small.

For the perpendicular sensor, where $ΔTsigmax$ was low, HER increased with the bubble concentration [Figs. 3(a)–3(c)]. Applying lower excitation pressures resulted in a ratio that changed approximately quadratically with the logarithm of bubble concentration. As the pressures increased, the relationship became approximately linear, shown by the increased R² values with the linear-fitted curve at pressures ≥0.69 MPa [Fig. 4(b)]. The slope of the linear-fitted HER-concentration curve also remained constant at pressures above this level [Fig. 4(a)]. On the other hand, for the concentric sensor, HER did not vary with bubble concentration.

Inspecting each frequency range, we note the different increase in harmonic energy with bubble concentration. For the perpendicular sensor, as the dilution ratio decreased from 1:1000k to 1:10k bubble-to-water, the 2nd harmonic energy increased by 6.1 ± 1.6 times, the 3rd harmonic energy increased by 4.3 ± 1.0 times, while the 7th–12th harmonic energy increased by only 1.6 ± 0.3 times (mean ± standard deviation from all sonication pressures). For the concentric sensor, with the same increase in bubble concentration, the 2nd harmonic energy increased by 1.1 ± 0.1 times, the 3rd harmonic energy increased by 1.4 ± 0.2 times, while the 7th–12th harmonic energy increased by 1.7 ± 0.4 times.

The relationship between HER and bubble concentration was demonstrated at two sensor locations. For the perpendicular sensor, HER correlated with bubble concentration. For this sensor, the observed time-domain signals were a superposition of many individual microbubble emissions [Figs. 2(a)–2(c)]. In each oscillation cycle, each bubble emitted a low-frequency emission associated with its relatively slow expansion, and a broadband impulse associated with its relatively rapid collapse-and–rebound action (Song et al., 2016). The emissions of different bubbles travelled to the sensor with time differences. The maximum time difference was expected to be 0.55 μs—longer than the duration of the collapse impulses but smaller than the wavelength of the expansion signal. Thus, as the bubble concentration increased, the collapse impulses did not add up coherently, and the narrow impulses observed at lower concentrations appeared as broader peaks with lower relative amplitudes at higher concentrations (leading compressional peak in each acoustic cycle). In contrast, the expansion signal added up coherently, resulting in low-frequency oscillations—the trailing compressional peaks and the rarefactional peaks—whose relative amplitudes increased at a faster rate with the bubble concentration. The higher rate of increase in lower frequencies with concentration was also shown in the frequency domain [Figs. 2(d)–2(f)].

Interestingly, for the perpendicular sensor, we observed that the frequencies up to 1.5 MHz increased at a fast rate with concentration despite $ΔTsigmax$ being at 0.55 μs, which implies that only frequencies below 0.9 MHz would be guaranteed to interfere in a constructive manner. This is likely because at frequencies close to or slightly higher than the low frequency range estimate calculated from $ΔTsigmax$⁠, emissions from a large proportion of the bubbles around the centre of the bubble cloud could still interfere constructively. In this light the low frequency bound calculated from $ΔTsigmax$ should be viewed as an estimate rather than a hard limit.

Since the increase in HER was caused by constructive interference of low frequencies, as opposed to random interferences of high frequencies, it was primarily driven by the numerator. Quantitatively, the rates of increase in the energy at different frequency bands depend on multiple factors, however, and to date we are not aware of a theoretical framework that could account for all relevant contributions, such as bubble instability or coalescence, which could predict the rate of these increases precisely. We also note that here we calculated the HER using the ratio of harmonics, but conceptually, ratio of harmonic-to-broadband signals could be used as well, provided that appropriate frequency ranges were used in the numerator and the denominator.

The HER-concentration relationship was pressure-dependent. At low pressures, the relationship was pressure-dependent and increased approximately quadratically with the log of bubble concentration. At the pressure of 0.83 MPa and above, the relationship was pressure-independent, exhibiting a linear relationship with the same slope regardless of the applied pressure (Figs. 3 and 4). It was unclear why the HER-bubble concentration relationship exhibited these trends. More investigations are needed to understand how bubble interactions may affect the emissions they produce, and how different frequencies emitted by different bubbles quantitatively add up. Nonetheless, given that the HER-concentration saturated at the sonication pressures of 0.83 MPa and above in our results, we speculate that this method has less variation for applications which use moderate to high acoustic pressures, such as HIFU, where a large portion of the focal volume experiences sufficiently high pressures that would make the HER-concentration relationship the same for the entire focal volume.

On the other hand, for the concentric sensor, the time differences of emissions from the bubbles were expected to be up to 4.8 μs. This corresponded to the frequency of 0.10 MHz. Thus, we expected the constructively interfering frequencies (≪ 0.10 MHz) to be outside the relevant frequencies emitted by the microbubbles. Indeed, HER calculated using the same frequency range as the perpendicular sensor showed no correlation with bubble concentration. Since no signal within the constructively interfering frequencies were available, the HER method could not be used to estimate the bubble concentration when the sensor was placed at this position. Thus, the definition of low and high frequencies used to calculate the HER and the utility of this method for estimating the concentration of active bubbles is dependent on the insonation centre frequency and $ΔTsigmax$ calculated for each setup.

With the appropriate experimental setup and conditions, the method above could have several uses. The HER-concentration relationship could be calibrated for each setup and used in subsequent experiments, or the HER values from different experiments performed under the same conditions could be compared to account for the potential differences in the concentration of cavitation events being produced. In animal applications or in vitro, there should be opportunities to use the method when microbubbles are confined in small regions, such as a vessel model, and when there is flexibility on where to place the acoustic sensor. Our method may also be useful for estimating the number of cavitation events produced by solid or liquid nucleation methods, such as nanodroplets (Kripfgans et al., 2000) and nanocups (Kwan et al., 2015), allowing a better comparison between these nucleation methods and the use of pre-formed microbubbles.

Cavitation is responsible for the many therapeutic and harmful bioeffects produced in ultrasound therapies. However, the concentration of cavitation bubbles produced during sonication has not yet been extractable from the acoustic emissions captured during the treatment. This information is important, because knowing whether an increase in acoustic emission energy is due to an increased number of cavitation bubbles or an increased energy of the cavitation bubbles would facilitate treatment monitoring. Using a theoretical framework of how emissions from different bubbles accumulate at a sensor, we described how a sensor placed at certain locations could capture signals where the magnitude ratio between the low and high frequency components correlates with the concentration of cavitation events. The ratio can be used for estimating the concentration of bubbles undergoing cavitation within a therapeutic ultrasound beam. However, it is important to note that the HER-concentration relationship is setup-dependent, as other factors, such as frequency-dependent attenuation or the shape and size of the bubble cloud, may affect the relative amplitudes of different frequencies as well.

This work was supported by Alzheimer's Research UK (ARUK-IRG2017A-7) and the Royal Thai Government Scholarship Program.