 Research
 Open Access
 Published:
Simulation of highintensity focused ultrasound lesions in presence of boiling
Journal of Therapeutic Ultrasoundvolume 4, Article number: 11 (2016)
Abstract
Background
The lesions induced by highintensity focused ultrasound (HIFU) thermal ablations are particularly difficult to simulate due to the complexity of the involved phenomena. In particular, boiling has a strong influence on the lesion shape. Thus, it must be accounted for if it happens during the pulses to be modeled. However, no acoustic model enables the simulation of the resulting wave scattering. Therefore, we propose an equivalent model for the heat deposition pattern in the presence of boiling.
Methods
Firstly, the acoustic field is simulated with kWave and the heat source term is calculated. Then, a thermal model is designed, including the equivalent model for boiling. It is rigorously calibrated and validated through the use of diversified ex vivo and in vivo data, including usually unexploited data types related to the bubble clouds.
Results
The proposed model enabled to efficiently simulate unitary pulses properties, including the sizes of the lesions, their morphology, the boiling onset time, and the influence of the boiling onset time on the lesions sizes.
Conclusions
In this article, the whole procedure of model design, calibration, and validation is discussed. In addition to depicting the creative use of data, our modeling approach focuses on the understanding of the mechanisms influencing the shape of the lesion. Further work is required to study the influence of the remaining bubble clouds in the context of pulse groups.
Background
In many industrial domains, numerical simulation has been adopted as a primary tool for optimizing processes and devices. However, the simulation of highintensity focused ultrasound (HIFU) treatments presents specific challenges. Thus, modeling the lesion creation process remains complex in many clinically relevant situations, where cavitation or boiling occur [1]. In these cases, where the influence of the bubble cloud cannot be directly simulated, equivalent models have to be developed [1, 2]. Beyond their predictive capability, which is by nature restricted, they aim at producing a “useful understanding” [1] by identifying the main physical phenomena shaping the lesion. In this context, this paper describes a rigorous approach for designing, calibrating, and validating an equivalent model to simulate the in situ temperature field.
Acoustic simulation of HIFU beams has been greatly studied during the past decades, and several wave equations have been derived, such as the widely used KhokhlovZabolotskayaKuznetsov (KZK) [3] equation and the Westervelt equation [4]. The choice of the acoustical model is not discussed in details in this paper, but the interested reader can for example refer to the book of Hamilton and Blackstock [5]. The key point is that the use of nonlinear models is recommended, as the presence of shock waves strongly enhances the heat deposition near the focal spot [6].
Depending on the acoustic power and the frequency, cavitation may also occur and strongly influence the lesion shape. It can mechanically damage the tissue, and it locally increases the heating [7, 8]. Cavitation mainly occurs near the focus, where the peak rarefactional pressure is bigger, or at the interfaces between tissues. An elegant work was reported by Chavrier et al. [2] for modeling this effect. The model is incorporated into the acoustic model as a local modification of the attenuation coefficient based on the density of microbubbles in the tissue.
It has also been reported in [9] that the nonuniform vibration of the transducer, due to surface waves, can increase the prefocal peak. If prefocal damages are under consideration, it may therefore be worth considering using hydrophone measurements and holographic reconstruction methods to properly predict the acoustic field.
The validation of acoustic simulations in clinically relevant situations is complex as hydrophones are delicate and not designed to work in biological tissues. Hence, validation is often performed in water, using needle hydrophones at low power or fiber optic hydrophones at higher power levels [10–12]. Two studies have also reported measurements respectively in oil [13] and in a gel phantom [9]. These measurement methods result in reliable validation of the numerical schemes. Nevertheless, validating the simulations in such wellcontrolled media circumvents the modeling of local heterogeneities and the uncertainties on tissue properties.
Many different thermal models have also been proposed during the last decades (see the review of Bhowmik et al. [14]) but the most widely used is Pennes bioheat equation [15]. In this model, perfusion is considered homogeneous in the tissue, and it is therefore illsuited to the thermal simulations near large blood vessels, especially if the temperature increase is moderate. Detailed discussion about the underlying hypotheses can be found in [16, 17], but, despite strong simplifications, Pennes’ model is generally considered to be accurate enough for a wide range of practical situations.
Thermal damages are usually computed using the thermal dose defined in [18] which is expressed in seconds and represents the equivalent time which would produce the same biological effects at a temperature of 43 °C. Its use for HIFU treatments raises major questions [19], notably its validity above 47 °C, but it is widely used and commonly accepted. A threshold of 14.4×10^{3} s is usually considered for thermal destruction but tissuedependent thresholds can be found in the review of Yarmolenko et al. [20].
When the temperature in the tissue is high enough, boiling may occur, resulting in the growth of a bubble cloud into the target. It was observed that the presence of boiling bubbles into the tissue induces a high impedance contrast which reflects most of the beam [21]. Therefore the heat deposition pattern is dramatically modified as the energy does not reach the postfocal area anymore. In particular, it is commonly reported that this causes a preferential growth of the lesion towards the transducer. This results in tadpoleshaped lesions for transducers with an fnumber close to unity of higher [8, 21–23]. Consequently, if boiling occurs during a treatment to be simulated, it is fundamental to consider it in the model.
To date, to the best of our knowledge, no acoustic model enables for simulating the resulting wave scattering. However, equivalent models have been proposed to model the changes in the heat deposition pattern [1, 21]. Wojcik et al. proposed to artificially modify the acoustical impedance of the tissue in the bubbly regions during the acoustical simulations [21]. By contrast, Meaney et al. accounted for the effect of boiling directly in the thermal model by considering that the bubble cloud shields the postfocal area [1]. All the energy which would have been deposited after the boiling zone was artificially deposited uniformly over a 0.5cmdiameter spherical volume centered on the most proximal part of the calculated boiling volume. The resulting lesions were found to be in good agreement with the experimental data.
Based on all these considerations, we designed a complete model to simulate the lesions resulting from HIFU pulses. The acoustical simulations are based on the kWave toolbox, and we use a “layerbylayer” approach to overcome the memory limitations associated with nonlinear 3D simulations. The heat source term calculation is also discussed in terms of energy conservation. A new thermal model, adapted from [1], is then proposed to account for the effect of boiling on heat deposition. Finally, the available experimental data are presented, which include original measurements related to the bubble clouds. Model calibration and validation are presented in the “Results and discussion” section, and the validity domain is subsequently discussed.
Methods
Materials
This modeling approach is applied to the particular example of thermal ablations using the Echopulse (Theraclion, Malakoff, France), which is used for the treatment of breast fibroadenomas [24] and thyroid nodules [25]. Benefitting from the important amount and diversity of data collected during preclinical ex vivo and in vivo experiments, the current paper focuses on the modeling of pulses on a rabbit liver.
The acoustic field is generated by a treatment head comprising a singleelement therapy transducer and a rectangular confocal imaging probe which enables the monitoring of the treatment using Bmode imaging. The therapy transducer is a spherical cap with a curvature radius of 38 mm and a diameter of 56 mm. The rectangular hole for the imaging probe however breaks the axial symmetry of the generated beam: consequently, threedimensional simulations are necessary.
Hereinafter, z denotes the main propagation axis, pointing from the transducer to the focus, x is parallel to the imaging probe, and y is orthogonal to it. The frequency used for the treatments is 3 MHz, which corresponds to a wavelength of about 0.5 mm.
At this frequency and with the acoustical powers used in the clinical configuration, which are below 150 W, no cavitation has been detected using a passive cavitation detector. Therefore, the effect of cavitation on the lesion creation is not considered in this work.
Acoustical simulations
At 3 MHz, considering a 60×60×50 mm domain and the simulation of the nonlinear propagation up to the tenth harmonic with the minimum of two points per wavelength, the pressure field at one given time step in single precision requires more than 40 GB of RAM, making the computations practically intractable on classical hardware.
Hence, we used a layerbylayer approach, taking advantage of the beam convergence. Actually, harmonics generation is cumulative, which implies that harmonics are of increasing importance as the beam width decreases. Therefore, we divided the computational grid into a few millimeters thick layers of decreasing width and increasing spatial and temporal resolution (see Fig. 1). The finer grid enabled for the propagation of 10 harmonics which guaranteed a reasonable accuracy according to previous fiberoptic hydrophone measurements in water. The layerbylayer approach introduces a quasioneway approximation which was considered acceptable as we did not want to simulate the imaging system.
All the acoustic simulations were performed using the compiled C version of kWave 1.1 [26]. Instead of considering a wave equation, it directly solves the coupled set of equations including mass conservation, momentum conservation and medium constitutive relation using a pseudospectral kspace firstorder method. In particular, it computes the velocity field, which is convenient for heat source computations. The properties of the different media are listed in Table 1. An example of acoustic field in water is reported in Fig. 2.
As discussed in the “Background” section, the validation is important but complex. In order to confront the kWave simulations to experimental data with the real tissues, we excised tissue layers comprising the skin, subcutaneous fat, and muscle from rabbits. They were then maintained on a rigid metallic frame and put into a degassed water tank. The acoustic field through these tissue layers was finally measured at low power using a needle hydrophone and compared to numerical simulations. We found a very good agreement between measurements and simulations over seven biological samples and 55 measurements for different tissue geometry although the biological tissue was modeled as a single homogeneous medium. An article reporting the details of the validation process is currently under review.
Heat source term
Heat source term is generally computed using the planewave approximation. Under the latter hypothesis, the temporalaveraged acoustic intensity can be defined as a scalar quantity:
with p _{A} the pressure amplitude, ρ _{0} the equilibrium density of the medium, and c _{s} the speed of sound. Then, in the linear case, the local expression for the power lost by the wave is
with α _{att} as the attenuation coefficient of the medium at the considered frequency.
In the general case, the temporalaveraged intensity is defined as
with p as the pressure field, v as the (vectorial) velocity field, and <.> denoting the temporal average over one period. An energy balance over an elementary volume combined with the divergence theorem leads to
In practical situations, however, when applied to kWave simulation results, the numerical computation of the divergence in the last formula using classical finite difference schemes leads to unphysical results, with a locally negative heat source term. Several numerical schemes have been tested, including highorder centered finite difference schemes and derivation in the spectral domain, without solving this issue. However, the cartesian grid is illsuited to the geometry of our beam, as the angular aperture is greater than $\frac {\pi }{2}$ , what might explain these results.
Therefore, as a tradeoff, we computed the acoustic power lost by the wave over each xy plane as $\frac {dP(z)}{dz} dz$ , with:
where the two last terms compensate for the power losses through the sides of the computational domain delimited by the planes x=x _{min}, x=x _{max}, y=y _{min}, y=y _{max}, z=z _{min}, and z=z _{max}. The power lost within each voxel layer was then spread according to the repartition of the heat source computed under the linear planewave approximation:
A comparison with the linear planewave approximation is provided in the “Results and discussion” section.
In reality, the energy lost by the wave is partly converted into heat (“absorbed”) and partly scattered by the tissue heterogeneities. Let α _{abs} and α _{scat}, respectively, denote the absorption and scattering coefficients:
Let A denote the ratio between absorption and attenuation:
The previously described expressions for the power lost by the wave must therefore be multiplied by A in order to compute the heat source term. The value of A is of great importance as it is a scaling factor on the heat source. Unfortunately, the values reported in the literature are very variable. For example, Pohlhammer et al. report that A is between 0.45 and 0.77 in the liver [27] while Goss et al. report a value of 0.31 [28]. As an uncertainty of factor 2 on the heat source term is not acceptable, A must be properly calibrated and not set according to the literature.
Thermal simulations
In this study, we chose to use the Pennes’ bioheat equation [15] due to its simplicity and the fact that perfusion values are available for most tissues of interest. Moreover, all in vivo pulses were delivered far from the main arteries. The equation was generalized to account for the energy consumed by protein denaturation and boiling based on published data. A simple dependance of the perfusion term with the thermal dose was also implemented. Finally, in the presence of boiling, an equivalent model was used to modify the heat source term: several parameters are introduced which aims at modeling the scattering of the ultrasound beam by the bubble cloud. As they are not directly related to tabulated physical values, they were calibrated. The details of this generalization are described below.
The proposed equation is reported below:
where T designates the tissue temperature, T _{b} the blood temperature, and ω _{b} the perfusion rate. ρ and C, respectively, refer to density and specific heat while the subscripts _{t} and _{b}, respectively, refer to tissue and blood. k _{t} is the tissue heat conductivity, D is the thermal dose in tissue, and Q is the heat source term. For the sake of concision, we note x=(x,y,z) as the dependence in space.
The equation was solved with Matlab using finite differences with a firstorder Euler explicit scheme in time and a centered scheme in space.
In order not to overestimate the tissue temperature, we modeled the energy consumed by the protein denaturation and water vaporization using an additional temperaturedependent specific heat.
with C _{t} as the temperaturedependent specific heat, C _{0}=3700 J.kg^{−1}.K^{−1}, C _{denat} as the profile reported in [29], C _{boil} as the profile adapted from [30].
The profiles used for C _{denat} and C _{boil} are based on differential scanning calorimetry (DSC) data. DSC consists in heating a sample of interest and a reference material following a predefined temperature profile. The differences in the power to supply to both samples reveal the reactions occurring in the sample of interest.
To the best of our knowledge, protein denaturation had not been considered in previously reported thermal simulations. However, according to the data reported by Lepock et al. [29] about a rat liver homogenate, it represents 22 kJ.kg^{−1}, which corresponds to about 10 % of the energy which is necessary to raise the temperature of the liver from 37 to 85 °C considering a constant specific heat.
For boiling, we considered that the tissue was made up of 75 % of water and that the vaporization was not performed at a constant temperature but progressively, following the profile reported in [30]. As the unit was not specified in the article, we normalized the data as follows:
with H _{vap}=2260 kJ.kg^{−1} as the enthalpy of the vaporization of water, C _{Rama} as the profile reported by Ramachandran et al. [30], and Θ as the integration variable corresponding to the temperature.
The temperature above which boiling induces additional heating is noted as T _{boil} and was set to 85 ° C to be consistant with the specific heat variations used.
When no boiling occurs, we note the heat source term Q _{0}. It is computed as previously described and convoluted with a 2D Gaussian kernel of standard deviation σ _{defoc}>0. This convolution aims at modeling the defocusing due to the medium heterogeneity and irregularities in the shape of the interfaces:
where ⊗ denotes the spatial convolution. If needed, the movement of the transducer is modeled by a simple translation of the heat source term, which is responsible for the dependence in time of Q and Q _{cons}.
During the simulations, when the temperature reached T _{boil}, the heat deposition profile Q was modified at each time step.
Let $\mathcal {B}(t)$ denote the set of grid points where T>T _{boil} at a given time t:

The heating is null within the bubble cloud itself. This aims at modeling the energy consumed by the additional mechanical and thermal tissue damages induced by the bubble cloud.
$$\forall \mathbf{x} \in \mathcal{B}(t), Q(\mathbf{x},t)=0 $$ 
A shielding coefficient r _{shield} is computed as the maximum proportion of the deposited power which is intercepted by the beam within all the xy planes:
$$r_{\text{shield}} = \max\limits_{\zeta}^{} \frac{\sum\limits_{(x,y,\zeta)\in\mathcal{B}}^{}Q(x,y,\zeta)}{\sum\limits_{x=x_{\text{min}}}^{x_{\text{max}}} \sum\limits_{y=y_{\text{min}}}^{y_{\text{max}}} Q(x,y,\zeta)} $$Let z _{shield} denote the z coordinate where the maximum is reached.

The zone $\mathcal {H}(t)$ where enhanced heating takes place is computed as a morphological dilatation of $\mathcal {B}(t)$ using Matlab imdilate() function with a spherical structuring element $\mathcal {S}$ of radius R _{ SE }:
$$\mathcal{H}(t) = \mathcal{B}(t) \oplus \mathcal{S}(R_{SE}) $$R _{ SE } is a parameter of the model to be calibrated.

The total power $P_{\mathcal {H}}$ to be distributed within $\mathcal {H}(t)$ is set as a fraction of the acoustical power reaching the plane z _{shield}:
$$P_{\mathcal{H}} = \eta_{\text{intercept}} \, r_{\text{shield}} \, P(z_{\text{shield}}) $$with η _{intercept} as the coefficient characterizing the efficiency of the power absorption, which has to be calibrated.

This power is then distributed according to the product of three weighting factors W _{1}, W _{2}, and W _{3}:
$${{}{\begin{aligned} Q_{\text{boil}}(\mathbf{x},t) = P_{\mathcal{H}}(t) \frac{W_{1}(\mathbf{x},t)\ W_{2}(\mathbf{x},t)\ W_{3}(\mathbf{x},t)}{\sum\limits_{x=x_{\text{min}}}^{x_{\text{max}}} \sum\limits_{y=y_{\text{min}}}^{y_{\text{max}}} \sum\limits_{z=z_{\text{min}}}^{z_{\text{max}}} W_{1}(x,y,z,t)\ W_{2}(x,y,z,t)\ W_{3}(x,y,z,t)} \end{aligned}}} $$
W _{1} is such that the deposited power decreases as the inverse of the distance to $\mathcal {B}$ .
$${{} {\begin{aligned} \forall h = (x,y,z) \in \mathcal{H}(t), W_{1}(x,y,z,t) = \frac{1}{min\left\{d(h,b) \colon b \in \mathcal{B}(t) \right\}} \end{aligned}}} $$ 
W _{2} is such that, beyond z _{shield}, the deposited power is multiplied by r _{shield}, which is smaller than 1 and models the shielding due to the bubble cloud.
$$W_{2}(x,y,z) = \left\{ \begin{array}{ll} 1 & \text{if} \ z\leq z_{\text{shield}} \\ r_{\text{shield}} & \text{if} \ z>z_{\text{shield}} \end{array} \right. $$ 
W _{3} is such that, if the top of the bubble cloud $\mathcal {B}$ is locally convex, heating is proportional to the source term in the absence of boiling. If the bubble cloud is locally concave, W _{3} is uniform but increased by a factor W _{+} within a cone having a θ _{cone} angular aperture and the apex at geometrical focus. This intends to model the preferential growing of the bubble cloud towards the transducer. Figure 3 illustrates the two types of deposition patterns, and the implementation details are presented below.
This part of the model is probably closely related to our transducer geometry and should not be considered for a generalization. We decided to consider these two different deposition patterns because each of them was not able to separately account for the observed lesion morphology. They were both qualitatively confirmed by kWave simulations when changing the properties of the medium around the focal spot to model the presence of different bubble cloud shapes. However, as kWave is designed for weakly heterogeneous media, we do not present these simulations. They were only a source of inspiration and do not fully justify our model, which must be considered as an equivalent model.
The top of the bubble cloud $\mathcal {B}(t)$ was considered locally convex at a given time t if :
$${{\kern17.1pt} {\begin{aligned} \min\left(z\mid(x,y,z)\in\mathcal{B}(t)\right) &= \min\left(z\mid(x,y,z)\in\mathcal{B}(t) \ \text{and}\right.\\ &\qquad\qquad\left.(x,y) = (x_{\text{foc}}(t),y_{\text{foc}}(t))\right) \end{aligned}}} $$where the subscript._{foc} refers to the coordinates of the actual focus. Then,
with $ \mathcal {C}(t) = {\fontsize {9.3pt}{6pt}{} {\begin {aligned}\left \{ (x,y,z) \mid z<z_{\text {foc}} \ \text {and}\ \text {atan}\left (\frac {\sqrt {(x^{2}+y^{2})}}{z}\right)\leq \theta _{\text {cone}}\right \} \end {aligned}}}$

Even if this equivalent model takes complex phenomena into consideration, it only introduces four additional parameters of simple signification:

R _{ SE }, the radius of the zone benefiting from additional heating.

η _{intercept}, the proportion of the intercepted power which contributes to the additional heating (the rest is supposed to be simply scattered).

W _{+}, the factor characterizing the additional heating due to the refocusing upstream to the bubble cloud when it tends to become concave.

θ _{cone}, the angular aperture of the cone where refocusing induces additional heating when the bubble cloud tends to become concave.
In the case where the HIFU transducer is a spherical cap without imaging probe, W _{+} and θ _{cone} are not necessary, which reduces the number of additional parameters to 2. The calibration and the relative importance of these parameters is assessed in the “Results and discussion” section.
Finally, thermal damages were assessed using the usual thermal dose formula [18]:
where τ is the time integration variable. We considered a threshold of 14.4×10^{3} s for tissue destruction. When the thermal dose reached this threshold, we considered that the blood vessels responsible for the perfusion were also coagulated. We modeled this dependence by a simple linear relation:
Experimental data
During preclinical studies, HIFU pulses of 4 s were performed ex vivo on a bovine liver. After degassing, the samples were immersed in degassed water at 37 °C. One hundred twenty individual pulses were delivered with a median acoustic power of 43.3 W at 3 MHz at a median depth of 16 mm.
The samples were then placed in a cassette which was vacuum sealed in a plastic bag in order for the liver to stick to the cassette walls and immediately frozen at −20 °C overnight. The samples were then sliced orthogonally to z, and blockface pictures were taken. Each lesion was then segmented by one operator in order to assess the lesion sizes. Due to the vacuum sealing, the lesion sizes along z were overestimated: therefore, we only relied on the lesions dimensions along x and y. The median boiling onset time was 3 s, and we considered the median lesion sizes over the 40 pulses having a boiling onset time between 2.5 and 3.5 s.
Thirty treatments of 19 pulses were also delivered in vivo on a rabbit liver in the same treatment conditions. The cutting protocol was modified, and no vacuum sealing was performed. Consequently the median size of the global lesion along the three dimensions was considered.
Treatments made up of several 8 and 12s pulses were also performed in vivo on the rabbit liver with the Echopulse. The pulses were delivered at a median depth of 14 mm with a median acoustic power of 49 W. During these pulses, the focus was describing a circular trajectory in the xy plane with diameters of, respectively, 1.1 mm for the 8s pulses and 1.3 mm for the 12s pulses.
In order to have threedimensional representations of the lesions morphology, the liver samples were frozen and sliced orthogonally to z. Each slice was subsequently photographed with a qualitatively constant lighting. On the photographs, the unitary pulses were distinguishable and the hypothesis was made that, due to the long cooling time between the pulses (approximately 40 s), the pulses were independent. Therefore, 10 unitary lesions were manually segmented, 5 of which resulting from 8s pulses and 5 from 12s pulses. The segmented contours were then stacked to obtain threedimensional representations of the unitary lesions.
Finally, an archetypal lesion was computed for each pulse type based on the individual 3D reconstructions: the five voxelized lesions of each type were then aligned based on their center of mass and each voxel was retained in the archetypal lesion if it belonged to at least four of the five individual lesions. The resulting 3D volume dimensions were found to be within 200 μm of the median lesion sizes computed over respectively 69 and 72 in vivo pulses, thus validating the representativeness of the segmented lesions.
As boiling occurred during the majority of the pulses, the boiling onset times were estimated based on videos of the Bmode imaging system. The hyperechoic marks (HEM) were also manually segmented right after each 12s pulse based on the Bmode images. This enabled us to evaluate the growth rate of the bubble cloud itself and to compare it to the simulations.
Calibration and validation
Firstly, the parameters for which typical values can be found in the literature were set to these values (see Table 1). The grid size was 71×71×111 with a spatial step of 200 μm. The time step was set to 0.05 s. The acoustic power was set to the average power used in the different experiments.
Then, we considered the 4s pulses and 12 s and tried to defocus the beam in order to obtain a boiling onset time of 3 s during the short pulses and 8 s during the long pulses by modifying the values of A (the ratio between absorption and attenuation) and σ _{defoc} (characterizing the defocusing).
Finally, the values of R _{ SE }, η _{intercept}, W _{+}, and θ _{cone} were set so that the simulated lesion matched the 12s lesion.
Results and discussion
Heat source term calculation
Figure 4 shows the evolution with z of the total deposited power per meter based on the acoustic simulation used for the 4s pulses. It was calculated as $\int _{x_{\text {min}}}^{x_{\text {max}}} \int _{y_{\text {min}}}^{y_{\text {max}}} Q(x,y,z) dx dy$ , with Q being computed as described in the “Methods” section for the linear planewave approximation and the proposed approach.
It was found that the differences between the two formulas were negligible. The main difference in the profile appears from 4 mm upstream to the focus, where the harmonics become significant and increase the absorbed power with the nonlinear formula. Besides, the curve obtained with the linear approximation is much smoother and less sensitive to the artefactual oscillations arising at the domain boundaries with the layerbylayer method. Finally, as the difference never exceeds a few percents, the choice of the formula should not have any influence on the simulated lesion. The method we describe could however be more accurate in particular cases where the heating due to nonlinearities becomes predominant, which could for example occur in water in the presence of strong shocks [5].
Controlability of HIFU lesions in the presence of boiling
For the clinically used 4s pulses, it must be noted that among hundreds of in vitro and in vivo HIFU pulses, a large proportion of which induced boiling, and no tadpoleshaped lesion was ever observed with the Echopulse. This is probably due to the particular geometry of the transducer which has a large aperture and a hole for the imaging probe. For comparison, the fnumber of our transducer is 0.68 while other articles report an fnumber close to unity or bigger: 1.06 [8], 1 [21], 1.79 [22], and 0.98 [23]. Thus, in our case, the ultrasonic wave comes from a wide solid angle to the side of the focal spot and the enhanced heating due to the bubble cloud activity is distributed all around the focus. This probably explains that the lesion does not suddenly grow towards the transducer.
Moreover, the color of the liver tissue is white when usual thermal damages occur but in the presence of boiling, it becomes dark brown [31]. This strongly burned zone never reaches the periphery of the lesion, which shows that boiling is confined to the lesion center (see Fig. 5).
Nevertheless, this does not mean that boiling has no influence on heat deposition, especially for longer pulses. As treatments are monitored using Bmode imaging, boiling is clearly visible as hyperechoic marks. It is therefore possible to determine when it occurs. For the 12s in vivo unitary pulses, we represent in Fig. 6 the size of the lesion versus the boiling onset time.
It clearly appears that boiling favors heat deposition as the lesions are bigger when boiling occurs early. One could argue that boiling is not the cause but a correlated consequence of variations in the power reaching the focus. However, as shown in Fig. 7, this tendency is also observed on single rabbits for neighboring pulses delivered at the same depth. Therefore, slight local defocusing effects are more likely to be responsible for the differences in the boiling time and the in situ power can reasonably be considered as constant, which shows that boiling is responsible for this enhanced heating.
Thermal simulations
Calibration of A and σ _{defoc}
To calibrate these two parameters, we used the boiling onset times for the 4 and the 12s pulses. The median experimental values were, respectively, 3 and 8.1 s. These experimental data were found to be adequate because they allow for a simple calibration process.
Indeed, the 4s pulses were delivered at fix positions while the focus was moved by the robotic head during the 12s pulses. Therefore, increasing A increases both boiling onset times but increasing σ _{defoc} mainly impacts the 4s pulses as the energy is already spread by the focus movement during the 12s pulses.
As a result, these two parameters were set to A=0.37 and σ _{defoc}=290 μm, which resulted in boiling onset times of, respectively, 3 and 7.9 s for the 4 and the 12s pulses. The value of A is higher than the value reported by Goss et al. [28] (0.31) but lower than what is reported by Pohlhammer et al. [27] ([0.45,0.77]). The defocusing induced by σ _{defoc} is low, which is consistent with the fact that the experiments were performed in the liver, which is qualitatively homogeneous.
Calibration of R _{ SE }, η _{intercept}, W _{+}, and θ _{cone}
The best agreement was obtained for R _{ SE }=2.5 mm, η _{intercept}=0.31, W _{+}=10, and $\theta _{\text {cone}} = \frac {\pi }{18}$ based on qualitative visual evaluation of the 3D lesions (see Fig. 8).
The influence of η _{intercept} is logically very important. The value of 0.31 seems physically acceptable when compared to A=0.37: on the one hand, the bubble cloud is generating high frequencies which tend to increase the absorption compared to classical heat deposition but, on the other hand, a proportion of this energy is used to “burn” the tissue and to mechanically damage it.
The influence of R _{ SE } on the final lesion is not very important although it has a significant impact on the computational time. As the additional power due to boiling decreases linearly with the distance to the boiling zone, increasing R _{ SE } above this value simply leads to consider grid points for which the additional heating is negligible. Decreasing R _{ SE } leads to smaller values of the coefficient η _{intercept} and is physically less plausible.
Increasing W _{+} increases the frequency of the shifting between the two heating modes described in Fig. 3. When W _{+} is lower, the bubble cloud can become visually concave, which has not been observed on the experimental Bmode images. In addition, W _{+}=10 enables for properly reproducing the height of the 12s lesion.
θ _{cone} must be small in order to model a refocusing above the bubble cloud. It must be noted that, for classical spherical transducers without imaging probe, the bubble cloud is very unlikely to become concave. Therefore, the expression of W _{3} can be restricted to the case of the locally convex bubble cloud, thus eliminating the need for W _{+} and θ _{cone}.
All calibration data are reported in Table 2.
Evolution of the lesion size with boiling onset time
To evaluate if the influence of boiling was properly modeled, we performed a linear regression on the lesions sizes with respect to the boiling onset times (see Fig. 6). The boiling onset time is noted tm for time to mark as it was detected based on hyperechoic marks on the Bmode images. The resulting slopes are expressed in mm s^{−1} and, as boiling tends to increase the efficiency of the heat deposition, they are negative: the later the boiling, the smaller the lesion. When the modulus of the slope is big, it means that boiling has a strong influence on the lesion size along this dimension.
Boiling onset time estimation was subject to uncertainties of the order of 0.1 s, and the precision on the lesion sizes was of the order of 0.1 mm; therefore, we performed orthogonal linear regressions. A classical linear regression would have lead to errors of more than 50 % on the slopes, which highlights the importance of a proper analysis of the data in the context of uncertainties.
For the 12s pulse simulations, the slopes characterizing the evolution of the lesion sizes along x, y, and z with the boiling onset time were, respectively, −0.21, −0.27, and −0.63 mm s^{−1}. The slopes obtained by regression on the experimental data were, respectively, −0.27, −0.17, and −0.25 mm s^{−1}. These data are reported in Table 3 with 90 % confidence intervals.
Along x and y, the slopes computed based on the model are in good agreement with the experiments as they are within the 90 % confidence interval. The experimental dependence of the size along z on boiling onset time is found to be smaller in the measurements than what was predicted by the model. The value derived from the model could be reduced by setting a smaller value for the coefficient W _{+}, but this discrepancy probably highlights the limit of the equivalent model. In general, the influence of boiling appears more complex to model along the z axis.
For the 4s pulse simulations, the slopes computed based on regressions on the experimental sizes along x and y were, respectively, −0.5 and −0.6 mm s^{−1}. The slopes obtained based on the simulations were, respectively, −0.6 and −0.7 mm s^{−1} (see Table 4). The values derived from the simulations were therefore in good agreement with the slopes estimated by regression on the experimental data. As described in the “Methods” section, the values along z were not considered for those pulses due to the slicing method which biased the sizes along z.
Growth of the bubble cloud for the 12s pulses
With these parameters, the final area of the bubble cloud in the xz plane for the 12s pulses is 5 mm^{2}, which equals the median experimental value evaluated based on a manual segmentation of the Bmode images. As for the lesion sizes, we also performed a linear regression on the bubble cloud area with respect to the boiling onset time. The slope was negative as, logically, the later boiling occurs, the smaller is the bubble cloud area. The estimated slope was −1.1 mm^{2} s^{−1} (see Fig. 9). By comparison, the slope obtained from simulated data was −1.2 mm^{2} s^{−1}. This further validates the fact that our equivalent model is relevant to simulate the influence of boiling on the heat deposition pattern.
8s pulses and groups of 4s pulses
Despite the good agreement between simulated and experimental lesions for the 4 and 12s pulses, we denoted a poor agreement for the 8s pulses. As shown on Fig. 10, the simulated lesion are significantly smaller than the segmented lesions. These discrepancies are likely to be due to the fact that the segmented lesions do not come from unitary pulses for the 8 and 12 s. In the case of the 8s pulses, the spacing between two adjacent lesions was 4 mm instead of 5 mm for the 12s pulses (see Table 5). This resulted in a greater overlap between neighboring lesions. Consequently, although they could be distinguished one from another, each pulse may have benefited from additional heating due to the presence of the bubble cloud created by the previous pulse. This hypothesis has been qualitatively confirmed by the retrospective visualization of captured movies acquired by the Bmode imaging system. In addition, the median interpulse cooling time was 48 s for the 12s pulses and only 39 s in median for the 8s pulses which could reinforce the influence of the adjacent pulses.
The same phenomenon occurs for groups of 4s pulses if we consider groups of 19 pulses with an important overlap. Even when these groups of 19 pulses were simulated including relevant cooling time, the height of the simulated lesion was 6 mm, which is smaller than the 9 mm experimentally measured. The thermal buildup, which was accounted for in the simulations, is therefore not the only factor which increases the size of adjacent lesions.
Validity domain
The last observations highlight the main limit of the proposed model. It is efficient for simulating unitary pulses but, in the context of whole treatments, it does not account for the additional heating due to remaining bubble clouds created by the previous pulses.
It must also be pointed out that the model has been specifically designed for the studied case. In the field of HIFU thermal ablation, the transducer geometry and treatment parameters can greatly vary depending on the clinical application: therefore, some phenomena which are not modeled here can be of significant importance for other devices. For instance, it has been shown that increasing the fnumber of the transducer increases the influence of thermal lensing [32]: as our transducer had a low fnumber of 0.68, it was considered to be negligible in our case.
The modifications of the tissue properties with the thermal dose were not modeled either, except for perfusion. Therefore, the lesion resulting from several pulses delivered at the same location would probably not be properly modeled.
More importantly, this model has been developed based on pulses delivered on liver, while the Echopulse is designed to treat breast fibroadenomas and thyroid nodules. Validating the simulations in several different animal tissues could however be relevant to evaluate the potential differences in lesion creation.
Conclusions
In this study, important issues regarding the modeling of HIFU thermal ablations are discussed. A complete numerical model of HIFU pulses was designed and carefully validated based on experimental data. Its calibration was thoroughly discussed, and we described a simple method to estimate the ratio between absorption and attenuation and the defocusing related to the medium heterogeneity.
An equivalent model, adapted from [1], was proposed to simulate the modified heat deposition pattern in the presence of boiling. First, we showed that boiling did not make the lesion sizes unpredictable with our transducer. Then, in order to fit the experimental lesion shape, it was found that the influence of boiling on the heat deposition was complex and could be modeled as an alternation between two regimes. The influence of the remaining bubbles clouds in the context of adjacent pulses was also highlighted. Thus, beyond the predictive capability of the equivalent model, which is by nature limited, this article illustrates the use of numerical simulation as a tool to confront our understanding of the lesion creation process to real data.
Moreover, the equivalent model was validated against cheap and usually unexploited experimental data including the evolution of bubble clouds areas evaluated on Bmode images and the slope characterizing the evolution of the lesion sizes with the boiling onset time. Therefore, this study is also intended to encourage the creative use of experimental data to stimulate the design of complete models of HIFU treatments, from the acoustic beam to the induced lesions.
Future works should include modeling the additional heating due to previously created bubble clouds in the context of whole treatments. Our approach could also benefit from more fundamental experiments aiming at measuring the acoustic field scattered by the bubble cloud.
Ethics, consent, and permissions
The herein reported in vivo experiments have been approved by the ComEth ANSES/ENVA/UPEC ethics committee with the reference number 2015030416468883.
References
 1
Meaney PM, Cahill MD, Ter Haar GR. The intensity dependence of lesion position shift during focused ultrasound surgery. Ultrasound Med Biol. 2000; 26:441–50.
 2
Chavrier F, Chapelon JY, Gelet A, Cathignol D. Modeling of highintensity focused ultrasoundinduced lesions in the presence of cavitation bubbles. J Acoust Soc Am. 2000; 108:432–40.
 3
Zabolotskaya EA, Khokhlov RV. Quasiplane waves in the nonlinear acoustic of confined beams. Sov Phys Acoust. 1969; 15:35–40.
 4
Westervelt PJ. Parametric acoustic array. J Acoust Soc Am. 1963; 35:535–7.
 5
Hamilton MF, Blackstock DT. Nonlinear acoustics. New York: Acoustical Society of America; 2008.
 6
Filonenko EA, Khokhlova VA. Effect of acoustic nonlinearity on heating of biological tissue by highintensity focused ultrasound. Acous Phys. 2001; 47:468–75.
 7
Hynynen K. The threshold for thermally significant cavitation in dog’s thigh muscle in vivo. Ultrasound Med Biol. 1991; 17:157–69.
 8
Khokhlova VA, Bailey MR, Reed JA, Cunitz BW, Kaczkowski PJ, Crum LA. Effects of nonlinear propagation, cavitation, and boiling in lesion formation by high intensity focused ultrasound in a gel phantom. J Acoust Soc Am. 2006; 119:1834–848.
 9
Canney MS, Bailey MR, Crum LA, Khokhlova VA, Sapozhnikov OA. Acoustic characterization of high intensity focused ultrasound fields: a combined measurement and modeling approach. J Acoust Soc Am. 2008; 124:2406–420.
 10
Khokhlova VA, Souchon R, Tavakkoli J, Sapozhnikov OA, Cathignol D. Numerical modeling of finiteamplitude sound beams: shock formation in the near field of a cw plane piston source. J Acoust Soc Am. 1998; 110:95–108.
 11
Tavakkoli J, Cathignol D, Souchon R, Sapozhnikov OA. Modeling of pulsed finiteamplitude focused sound beams in time domain. J Acoust Soc Am; 104:2061–072.
 12
Liebler M, Dreyer T, Riedlinger RE. Modeling of interaction between therapeutic ultrasound propagation and cavitation bubbles. Ultrasonics. 2006; 44:319–24.
 13
Wang K, Teoh E, Jaros J, Treeby BE. Modelling nonlinear ultrasound propagation in absorbing media using the kWave toolbox: experimental validation. In: IEEE International Ultrasonics Symposium: 7–10 October 2012. Dresden: 2012. p. 523–6.
 14
Bhowmik A, Singh R, Repaka R, Mishra SC. Conventional and newly developed bioheat transport models in vascularized tissues: a review. J Theor Biol. 2013; 38:107–25.
 15
Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948; 1:93–122.
 16
Wissler EH. Pennes’ 1948 paper revisited. J Appl Physiol. 1998; 85:35–41.
 17
Nelson DA. Invited editorial on Pennes’ 1948 paper revisited. J Appl Physiol. 1998; 8523:2–3.
 18
Sapareto SA, Dewey WC. Thermal dose determination in cancer therapy. Int J Radiat Oncol Biol Phys. 1984; 10:787–800.
 19
Shaw A, ter Haar GR, Haller J, Wilkens V. Towards a dosimetric framework for therapeutic ultrasound. Int J Hyperth. 2015; 6736:1–11.
 20
Yarmolenko PS, Moon EJ, Landon C, Manzoor A, Hochman DW, Viglianti BL, Dewhirst MW. Thresholds for thermal damage to normal tissues: an update. Int J Hyperth. 2011; 27:320–43.
 21
Wojcik G, Mould JJ, Abboud N, Ostromogilsky M, Vaughan D. Nonlinear modeling of therapeutic ultrasound In: IEEE, editor. IEEE Ultrasonics Symposium Proceedings. Seattle, WA: IEEE: 1995. p. 1617–622.
 22
Watkin NA, Ter Haar GR, Rivens I. The intensity dependence of the site of maximal energy deposition in focused ultrasound surgery. Ultrasound Med Biol. 1996; 22:483–91.
 23
Jensen CR, Ritchie RW, Gyöngy M, Collin JRT, Leslie T, Coussios CC. Spatiotemporal monitoring of highintensity focused ultrasound therapy with passive acoustic mapping. Radiology. 2012; 262:252–61.
 24
Kovatcheva R, Guglielmina JN, Abehsera M, Boulanger L, Laurent N, Poncelet E. Ultrasoundguided highintensity focused ultrasound treatment of breast fibroadenomaa multicenter experience. J Ther Ultrasound. 2015;1. http://jtultrasound.biomedcentral.com/about.
 25
Esnault O, Franc B, Leenhardt L, Rouxel A, Ménégaux F, Lacoste F. Highintensity focused ultrasound (Hifu) treatment for thyroid nodules: experimental and first clinical studies. In: 6th International Symposium on Therapeutic Ultrasound. Melville, NY: AIP: 2007. p. 98–103.
 26
Treeby BE, Cox BT. kWave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. J Biomed Opt. 2010; 15:021314–112.
 27
Pohlhammer JD, Edwards CA, O’Brien WD Jr.Phase insensitive ultrasonic attenuation coefficient determination of fresh bovine liver over an extended frequency range. Med Phys. 1981; 8:692–4.
 28
Goss SA, Frizzell LA, Dunn F. Ultrasonic absorption and attenuation in mammalian tissues. Ultrasound Med Biol. 1979; 5:181–6.
 29
Lepock JR, Frey HE, Ritchie KP. Protein denaturation in intact hepatocytes and isolated cellular organelles during heat shock. Eur J Cell Biol. 1993; 122:1267–276.
 30
Ramachandran T, Sreenivasan K, Sivakumar R. Water vaporization from heated tissue: an in vitro study by differential scanning calorimetry. Lasers Surg Med. 1996; 19413415:413–5.
 31
Mast TD, Makin IRS, Faidi W, Runk MM, Barthe PG, Slayton MH. Bulk ablation of soft tissue with intense ultrasound: modeling and experiments. J Acoust Soc Am. 2005; 118:2715–724.
 32
Connor CW, Hynynen K. Bioacoustic thermal lensing and nonlinear propagation in focused ultrasound surgery using large focal spots: a parametric study. Phys Med Biol. 2002; 47:1911–928.
Acknowledgements
We would like to deeply thank Jeremie Anquez, Quentin Jacob, and Nesrine Barnat for providing us with the experimental data and Björn Gerold for his comments about the manuscript.
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
AG conducted the study, designed the model, analyzed the data, and wrote the manuscript. SY supervised the physical aspects of the study, reviewed the manuscript, and was responsible for the final approval. VL supervised the data analysis, reviewed the manuscript, and was responsible for the final approval. PL supervised the numerical modeling, reviewed the manuscript, and was responsible for the final approval. All authors read and approved the final manuscript.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 HIFU
 HITU
 Modeling
 Thermal model
 Nonlinear acoustics
 Numerical simulation
 Bioheat
 Boiling