Numerical Simulations of the Nonlinear Interaction of a Bubble Cloud and a High Intensity Focused Ultrasound Field

: We studied the e ﬀ ects of a small bubble cloud located at the pre-focal area of a high-intensity focused ultrasound ﬁeld. Our objective is to show that bubbles can modify the bioe ﬀ ects of an ultrasound treatment in muscle tissue. We model a three-dimensional ultrasound ﬁeld in an idealized conﬁguration of real operating conditions. Simulations are performed using a combined method based on the Khokhlov-Zabolotskaya-Kuznetsov equation, describing the ultrasound propagation, and a Rayleigh-Plesset equation, modeling the bubble oscillations. The nonlinear interaction of the ultrasound ﬁeld and the bubble oscillations is considered. Results with and without bubbles for di ﬀ erent void fractions of the cloud and di ﬀ erent acoustic powers are compared. The cloud induces scattering, nonlinear distortion, and shielding of ultrasound, which increase the mechanical index in the pre-focal zone, shift the location, reduce the size, and modify the shape of the volume of tissue of high mechanical index values, and lower the pressure at the intended focus considerably. Although some hypothesis and parameters used in the models do not ﬁt the real HIFU situations, the simulation results suggest that the e ﬀ ects caused by a bubble cloud located in the pre-focal area should be considered and monitored to ensure the safety of high-intensity focused ultrasound treatments.


Introduction
Magnetic resonance imaging (MRI)-guided focused ultrasound (FUS) is a method for noninvasive surgery that provides temperature feedback to the operator for the control of the acoustic parameters of the treatment. This technique has successful and potentially therapeutic applications, such as the treatment of uterine fibroids [1,2], bone metastases [3], and prostate [4] and breast cancer [5]. Moreover, the therapy modality was used for the treatment of neurological conditions, such as chronic pain [6], and essential tremors [7,8]. During treatment, the induced thermal exposure is determined via thermal imaging using MRI [9] and is used to adjust the output power and sonication duration for the subsequent FUS sonications. Thus, the desired treatment can be achieved while reducing the need for prediction of exposure to ultrasound during the treatment planning phase.
Ultrasonic cavitation exists when a high amplitude ultrasound wave travels in a liquid or a tissue. As a result, bubbles can be generated, excited to oscillation, and even collapse at inertial cavitation [10]. This can give way to very high temperature elevations at microscopic scale [10]. The characterization of a cavitation field generated by a high amplitude ultrasound wave is complex and was the subject of many experimental and theoretical studies [10][11][12][13][14][15].
It was experimentally observed that cavitating bubbles can appear in a high-intensity focused ultrasound (HIFU) field used in therapy and it was shown that the spatial distribution of the bubbles

Materials and Methods
We consider the three-dimensional problem shown in Figure 1 written in cylindrical coordinates (z, r, θ), which are the axial coordinate, radial distance, and azimuth angle, respectively [32,41,44]. The HIFU field is excited by the focused transducer. Ultrasound travels in water and transmits into the muscle tissue. Cavitating bubbles are created and agglomerate in the pre-focal volume, forming a bubble cloud. During the propagation of ultrasound at high acoustic-pressure level through water and tissue, the wave is affected by nonlinear distortion and attenuation. When the sound wave interacts with the cloud located in the pre-focal zone, the bubbles exhibit nonlinear oscillations and the ultrasonic field undergoes changes due to bubble volume variations, since the bubbly medium is dispersive and nonlinear. The process can be described by coupling the pressure field and the oscillating bubble field. In addition, the bubble agglomerate provokes scattering and shielding of the acoustic field.
The ultrasonic field generated by the sonication at frequency f is defined by the acoustic pressure ( ) θ p z,r, ,t and the bubble oscillations field is defined by the bubble-volume variation ( ) θ v z,r, ,t , both being functions of time t . This three-dimensional problem is extremely complex to solve and would not be manageable for computations in practice. It is thus simplified by setting the variable θ as the angle of symmetry to reduce the mathematical problem to the ( ) z,r plane. During the propagation of ultrasound at high acoustic-pressure level through water and tissue, the wave is affected by nonlinear distortion and attenuation. When the sound wave interacts with the cloud located in the pre-focal zone, the bubbles exhibit nonlinear oscillations and the ultrasonic field undergoes changes due to bubble volume variations, since the bubbly medium is dispersive and nonlinear. The process can be described by coupling the pressure field and the oscillating bubble field. In addition, the bubble agglomerate provokes scattering and shielding of the acoustic field.
The ultrasonic field generated by the sonication at frequency f is defined by the acoustic pressure p(z, r, θ, t) and the bubble oscillations field is defined by the bubble-volume variation v(z, r, θ, t), both being functions of time t. This three-dimensional problem is extremely complex to solve and would not be manageable for computations in practice. It is thus simplified by setting the variable θ as the angle of symmetry to reduce the mathematical problem to the (z, r) plane. Acoustic pressure and bubble-volume variation are now functions of three variables: p(z, r, t) and v(z, r, t).
The procedure developed for solving the problem is structured by splitting the simulation model into two coupled parts ( Figure 1a): The acoustical problem is simplified by using the KZK approximation to compute the sound propagation between the transducer and a zone close to the bubble cloud (denoted by Model 1), whereas the bubble-sound interaction is simplified by using the RP approximation in the vicinity of the focus (denoted by Model 2). Model 1 simulates the nonlinear propagation of the focused ultrasound in the water and muscle tissue from the source z = 0 at time t = 0 up to the interface between both models at the distance z s at time t s [32], whereupon the simulated pressure field is used as an input for Model 2, which simulates the nonlinear interaction of the ultrasound wave and bubble oscillation field in the volume of soft tissue surrounding the bubble cloud [41].
Model 1 computes the propagation of nonlinear ultrasound in the domain 0 ≤ z ≤ z s , 0 ≤ r ≤ R from t = 0 to t = t s by simulating the propagation with the KZK equation: where ρ t,w and c t,w are the equilibrium density and the small-amplitude sound speed of the homogeneous medium, respectively, t = t − z/c t,w is the retarded time, δ t,w is the acoustic diffusivity related to the absorption, β t,w is the coefficient of nonlinearity of the medium, and ∇ 2 ⊥ denotes the Laplacian operator in the perpendicular plane to the propagation (index t refers to muscle tissue and index w refers to water). In both media, the attenuation follows the power-law α t,w (f) = α 0t,w f α 0yt,w , where α 0t,w and α 0yt,w are specific parameters of the propagation medium. Details of the implementation of the model can be found in the previous work [32]. For parameter values used in Model 1, see Figure 1a and Section 3.
Model 2 assumes that the spherical gas bubbles are contained in a small volume of muscle next to the focus of the HIFU field. The void fraction in that volume is V f . The time-dependent acoustic-pressure output g s (r, t) = p(z s , r, t) from Model 1 serves as input function for Model 2, which simulates the domain z s ≤ z ≤ Z, 0 ≤ r ≤ R from t = t s to the last instant of the study t = T. The following coupled system of two differential equations, composed of the RP equation (written in variable v) and the wave equation [41][42][43], is solved after adding the adequate initial and boundary conditions (system at rest at time t s , imposition of the incident ultrasound field from Model 1 at z = z s , axial symmetry, and non-reflecting boundary conditions at z = Z and r = R) [41,45]: In these Equations, ω 0 is the resonance frequency of the bubbles (f 0 ), δ = 4µ k /ω 0 R 2 0 is the viscous damping coefficient of the bubbly medium (µ k is the kinematic viscosity of the tissue), N(z, r) is the bubble density in the liquid, η = 4πR 0 /ρ t , the nonlinear coefficient associated to the adiabatic gas law is a = (γ + 1)ω 2 0 /2v 0 , the nonlinear coefficient associated to the dynamic response of bubbles is b = 1/6v 0 , R 0 and v 0 are the initial bubble radius and volume, respectively, and γ is the specific heat ratio of the gas. The difference between the current total volume and v 0 is v. Details of the implementation of the model can be found in the previous works [41,46]. For parameter values used in Model 2, see Figure 1a and Section 3.
Model 1 takes the effects of diffraction, attenuation, and nonlinear propagation of ultrasound into account. The parabolic approximation for diffraction used to derive the model limits its validity to waves that mainly propagate in one preferred direction and sound beams with a moderate angular spectrum [32,47,48]. It is thus adapted to the focused propagation of ultrasound in the homogeneous media considered here (water and tissue without bubbles). Model 2 takes into account the nonlinear interaction of both the acoustic field and the bubble oscillations by considering p and v as unknown variables of the differential problem [41]. This matches the requirements to simulate the coupling of both fields in the pre-focal zone with bubbles. The model assumes restrictions that limit its validity to monodisperse spherical bubbles in populations with low void fraction and, since the nonlinear parameter of the bubbly medium is several orders of magnitude higher than that of the tissue [42,43], it does not take into account the nonlinear propagation effects in the medium where no bubbles are present [10,42,43,48,49]. We thus set z s close to the bubble cloud.
The coupling of both numerical models at z s requires the synchronization of the signal at the discretization points in the radial direction. This process does not add any numerical imprecisions. This task is made easier in Section 3 since the point number per wavelength in the radial direction is the same for both numerical models.

Results
The simulations were performed using a focused transducer with F-number of 2, diameter D of 5 cm, and radius of curvature of 10 cm, operating at f = 1 MHz that was positioned to deliver a focus at a depth of 3.6 cm in the muscle tissue (c t = 1609 m s −1 , ρ t = 1050 kg m −3 , α 0t = 1.0048 × 10 −6 Np Hz −1.1 , α 0yt = 1.1, and β t = 4 kg 3 m −4 s −2 ) through a water path (c w = 1522 m s −1 , ρ w = 1000 kg m −3 , α 0w = 2.88 × 10 −16 Np Hz −2 , α 0yw = 2, and β w = 3.5 kg 3 m −4 s −2 ) [32]. The acoustic power was varied as P = 25-300 W. Following Section 1, an idealized air-bubble cloud of radius 0.15 cm was located at 0.91 cm in front of the focus, in which bubbles such that R 0 = 35 µm were considered with a void fraction in the range V f = 0-0.45% (N = 0-2.5 × 10 10 m −3 ) [18,21]. The other parameter values of the cloud were: Sound speed in the air c g = 340 ms −1 , air density ρ g = 1.29 kgm −3 , γ = 1.4, δ = 0.0056, η = 5.98 × 10 −7 m 4 kg −1 , a = 3.91 × 10 23 Hz 2 m −3 , and b = 3.18 × 10 11 m −3 [41]. The geometry of the simulation is shown in Figure 1a. Note that this configuration (geometry, ultrasound, cloud, and media) is set to define an idealized model simulation able to show the effects a bubbly cloud present in the pre-focal region induces on a HIFU field. However, other parameter values could be employed. For the simulations we used 64 points per wavelength radially for Model 1 and Model 2, 2 and 64 points per wavelength axially for Model 1 and Model 2, respectively, 100 and 99 points per cycle for Model 1 and Model 2, respectively, and 19 periods were computed with Model 2. The effect of nonlinear propagation in water and tissue on pressure is seen in Figure 1b by displaying the pressure waveform situated at z = 6.48 cm on the symmetry axis in front of the cloud at low and high powers (P = 25 W and P = 300 W) before transmitting into the bubble cloud. The nonlinear distortion undergone by the signal due to the propagation through the water and muscle tissue at high power is clearly observed by comparing with the sine signal obtained at low power. The linear wave obtained at low power and the nonlinear wave obtained at high power will affect the bubble oscillations in the cloud in different ways. This cloud is seen in Figure 1c, which shows the void fraction distribution in the focal area of the axial (z, r) plane. The MI defined as in which f is given in MHz and p n , the minimum amplitude of the negative pressure, is given in MPa [50], was studied under different values of P and V f . The impact of the bubble cloud with high void fraction on the pressure field at high power is presented in Figure 2b. The same case without cloud is shown in Figure 2a. We observe that harmonics are focused in different ways without bubbles, when the nonlinearity is due to the water and tissue only (Figure 2a), and with bubbles, when scattering due to the bubble cloud exists and the nonlinearity is due to the water, tissue, and bubbles (Figure 2b). The volumes of highest energy at each frequency strongly differ from one case to another. With bubbles, the distribution of energy mainly concentrates in front of the cloud and inside the cloud, whereas the harmonics produce an elongated volume of high energy around the axis without bubbles. Note that the depth of the pressure peak of each harmonic depends on all the parameters of the problem (ultrasonic field, cloud, nonlinear interaction of oscillating bubbles and ultrasound, scattering, and dispersion).
fraction distribution in the focal area of the axial ( ) z,r plane. The MI defined as in which f is given in MHz and n p , the minimum amplitude of the negative pressure, is given in MPa [50], was studied under different values of P and f V .
The impact of the bubble cloud with high void fraction on the pressure field at high power is presented in Figure 2b. The same case without cloud is shown in Figure 2a. We observe that harmonics are focused in different ways without bubbles, when the nonlinearity is due to the water and tissue only (Figure 2a), and with bubbles, when scattering due to the bubble cloud exists and the nonlinearity is due to the water, tissue, and bubbles (Figure 2b). The volumes of highest energy at each frequency strongly differ from one case to another. With bubbles, the distribution of energy mainly concentrates in front of the cloud and inside the cloud, whereas the harmonics produce an elongated volume of high energy around the axis without bubbles. Note that the depth of the pressure peak of each harmonic depends on all the parameters of the problem (ultrasonic field, cloud, nonlinear interaction of oscillating bubbles and ultrasound, scattering, and dispersion). The effects produced at high power by the bubble cloud with high void fraction on the MI are evidenced in Figure 3b by comparing to the same case without bubbles shown in Figure 3a. On one hand, we observe that the volume of tissue in which MI values are the highest is much smaller and The effects produced at high power by the bubble cloud with high void fraction on the MI are evidenced in Figure 3b by comparing to the same case without bubbles shown in Figure 3a. On one hand, we observe that the volume of tissue in which MI values are the highest is much smaller and is slightly shifted toward the front of the cloud. This change is provoked by the scattering due to the bubble cloud, as shown by the fringes seen in Figure 3b, and by the creation of highest pressure harmonics in this zone due to the nonlinearity of the oscillating bubbles ( Figure 2). MI reaches higher values in the muscle with the cloud than without bubbles. On the other hand, a large volume of low MI values is observed beyond the cloud. The ultrasonic signal is drastically lowered there because of the shielding effect of the cloud on the incident focused ultrasound. The peak-to-peak pressure value at z = 7.84 cm exhibits a loss of 69%. The effects produced at high power by the bubble cloud with high void fraction on the MI are evidenced in Figure 3b by comparing to the same case without bubbles shown in Figure 3a. On one hand, we observe that the volume of tissue in which MI values are the highest is much smaller and is slightly shifted toward the front of the cloud. This change is provoked by the scattering due to the bubble cloud, as shown by the fringes seen in Figure 3b, and by the creation of highest pressure harmonics in this zone due to the nonlinearity of the oscillating bubbles ( Figure 2). MI reaches higher values in the muscle with the cloud than without bubbles. On the other hand, a large volume of low MI values is observed beyond the cloud. The ultrasonic signal is drastically lowered there because of the shielding effect of the cloud on the incident focused ultrasound. The peak-to-peak pressure value at = z 7.84 cm exhibits a loss of 69%.
(a) (b) The influence of acoustic power on MI values without bubbles and when the bubble cloud is set at high void fraction is displayed in Figure 4a. The difference of maximum MI values obtained with and without the bubble cloud increases as the sonication power is raised. MI reaches 6.04 at 300 W with bubbles, 20% above the value obtained without cloud at the same power. The position at which MI reaches its highest value at 300 W with bubbles is located at = n z 7.57 cm. It is shifted 1.82 mm toward the front of the cloud from its position obtained without bubbles. The effect of void fraction on the MI at high acoustic power is represented in Figure 4b by modifying the bubble density N in The influence of acoustic power on MI values without bubbles and when the bubble cloud is set at high void fraction is displayed in Figure 4a. The difference of maximum MI values obtained with and without the bubble cloud increases as the sonication power is raised. MI reaches 6.04 at 300 W with bubbles, 20% above the value obtained without cloud at the same power. The position at which MI reaches its highest value at 300 W with bubbles is located at z n = 7.57 cm. It is shifted 1.82 mm toward the front of the cloud from its position obtained without bubbles. The effect of void fraction on the MI at high acoustic power is represented in Figure 4b by modifying the bubble density N in the cloud. The MI values remain almost unperturbed at small V f < 0.24% and they increase rapidly as N rises above this threshold. the cloud. The MI values remain almost unperturbed at small < f V 0.24 % and they increase rapidly as N rises above this threshold.

Discussion
Although some hypothesis and parameters used in the simulations do not fit exactly the ones in real HIFU situations, the results presented in this paper suggest that a bubble cloud located in the pre-focal area of a HIFU field can have a considerable effect on the ultrasound field. The interface between the soft tissue and the bubble cloud can generate scattering of the ultrasonic field, which can also be affected by nonlinear distortion. These phenomena can increase the MI values, modify the shape, change the location, and reduce the size of the volume of tissue of high MI values. Ultrasonic pressure is lowered beyond the bubble cloud. These conclusions are in concordance with experimental and numerical results available in the literature [15,18,21]. The increase in the MI value by 20% due to the bubble cloud observed here could have strong bioeffects, since the risk of damage in cell tissue from ultrasound (higher thermal and mechanical effects) would rise. Despite the restrictions of the theoretical models and their differences with experimental studies, the main results obtained here concerning the phenomenon of shielding and shifting of the intense volume of treatment agree with experimental data. Bubbles at a pre-focal position close to the theoretical focus induce intense rarefaction pressure in front of the bubbles and act as a barrier for ultrasound propagation. These results are important for HIFU treatments and suggest that bubble generation by cavitation should be monitored since the modification of the location, size, and shape of the highest MI volume observed in simulations suggest possible issues for safety purpose due to unintended damage of surrounding pre-focal tissue.
It must be noted that we chose to illustrate the use of the model with muscle tissue [32] in this idealized configuration, but the model can be applied to other tissues (liver tissue, other soft tissues, etc.) in which malignancies are treated by HIFU [50 (Section 7, Chapter VIII, and citations therein)]. Also, other configurations different from the idealized case presented here need to be studied. It must be noted that the complexity of the problem, due in particular to the dispersion and nonlinearity of bubbly media, makes it possible to consider an infinity of different situations by using other bubble sizes, densities, geometries, source frequency to bubble resonance ratios, etc. Here we wanted to show the effects of a bubbly cloud on a HIFU field and to stress the necessity of taking these phenomena into account, but many parameters can be modified and many other varieties of situations can be studied. This specific model simulation is only the first step to understand the effect of the presence of a bubble cloud on the HIFU field. Other transducers with lower and higher Fnumbers and other clouds with different geometries and bubble size distributions should be analyzed. The current simulation model includes restrictions (such as diffraction approximation of the KZK equation, monodisperse bubble requirement of the RP equation, no shear waves in tissues, etc.) that limit the validity of the simulations compared with real HIFU situations. It could be

Discussion
Although some hypothesis and parameters used in the simulations do not fit exactly the ones in real HIFU situations, the results presented in this paper suggest that a bubble cloud located in the pre-focal area of a HIFU field can have a considerable effect on the ultrasound field. The interface between the soft tissue and the bubble cloud can generate scattering of the ultrasonic field, which can also be affected by nonlinear distortion. These phenomena can increase the MI values, modify the shape, change the location, and reduce the size of the volume of tissue of high MI values. Ultrasonic pressure is lowered beyond the bubble cloud. These conclusions are in concordance with experimental and numerical results available in the literature [15,18,21]. The increase in the MI value by 20% due to the bubble cloud observed here could have strong bioeffects, since the risk of damage in cell tissue from ultrasound (higher thermal and mechanical effects) would rise. Despite the restrictions of the theoretical models and their differences with experimental studies, the main results obtained here concerning the phenomenon of shielding and shifting of the intense volume of treatment agree with experimental data. Bubbles at a pre-focal position close to the theoretical focus induce intense rarefaction pressure in front of the bubbles and act as a barrier for ultrasound propagation. These results are important for HIFU treatments and suggest that bubble generation by cavitation should be monitored since the modification of the location, size, and shape of the highest MI volume observed in simulations suggest possible issues for safety purpose due to unintended damage of surrounding pre-focal tissue.
It must be noted that we chose to illustrate the use of the model with muscle tissue [32] in this idealized configuration, but the model can be applied to other tissues (liver tissue, other soft tissues, etc.) in which malignancies are treated by HIFU [50]. Also, other configurations different from the idealized case presented here need to be studied. It must be noted that the complexity of the problem, due in particular to the dispersion and nonlinearity of bubbly media, makes it possible to consider an infinity of different situations by using other bubble sizes, densities, geometries, source frequency to bubble resonance ratios, etc. Here we wanted to show the effects of a bubbly cloud on a HIFU field and to stress the necessity of taking these phenomena into account, but many parameters can be modified and many other varieties of situations can be studied. This specific model simulation is only the first step to understand the effect of the presence of a bubble cloud on the HIFU field. Other transducers with lower and higher F-numbers and other clouds with different geometries and bubble size distributions should be analyzed. The current simulation model includes restrictions (such as diffraction approximation of the KZK equation, monodisperse bubble requirement of the RP equation, no shear waves in tissues, etc.) that limit the validity of the simulations compared with real HIFU situations. It could be improved by considering, e.g., other nonlinear models for the propagation of ultrasound and bubble oscillations that would allow us to simulate the high-amplitude fields in more realistic bubbly media, with, in particular, lower bubble sizes. Further experimental research is required in real treatment cases. The effects caused by bubbles located in the water between the transducer and the skin should also be studied.

Conclusions
Numerical simulations have been performed to provide information on pre-focal effects due to a bubble cloud on a HIFU field. The numerical tool used here combines two existing numerical models and considers a nonlinear ultrasonic signal interacting with a population of nonlinearly oscillating bubbles confined in a cloud. The acoustic power of the focused transducer and the void fraction of the bubble cloud have been increased in the analysis. The numerical results suggest that a bubble cloud located in the pre-focal area can induce scattering and nonlinear distortion of the ultrasonic field that have important effects: The MI values are increased; the shape and location of the volume of tissue of high MI values are modified and its size is reduced; the pressure at the intended focus is lowered considerably. Although some hypothesis and parameters of the model and simulations do not fit exactly the real HIFU situations, this information is potentially of interest for HIFU treatments since the presence of a bubble cloud may have implications for the tissue in the pre-focal zone. Funding: This research was funded by the National Agency for Research of Spain and the European Regional Development Fund, grant number DPI2017-84758-P.