Simulation of highintensity focused ultrasound lesions in presence of boiling
 Anthony Grisey^{1, 2}Email author,
 Sylvain Yon^{1},
 Véronique Letort^{2} and
 Pauline Lafitte^{2}
https://doi.org/10.1186/s4034901600569
© Grisey et al. 2016
Received: 4 September 2015
Accepted: 17 March 2016
Published: 31 March 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.
Keywords
HIFU HITU Modeling Thermal model Nonlinear acoustics Numerical simulation Bioheat BoilingBackground
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.
Medium acoustic and thermal properties
Water  Superficial tissue (skin...)  Liver  

Sound speed [ m.s ^{−1}]  1447  1547  1614 
Density [ kg.m ^{−3}]  1000  1214  996 
Attenuation [ dB.cm ^{−1}.MHz^{−y }]  2.2×10^{−3}  0.62  0.37 
y  2  1.27  1.27 
Specific heat [ J.kg ^{−1}.K ^{−1}]  ×  ×  3700 
Thermal conductivity [ W.m ^{−1}.K ^{−1}]  ×  ×  0.55 
Perfusion [ s ^{−1}]  ×  ×  0.018 
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
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.
A comparison with the linear planewave approximation is provided in the “Results and discussion” section.
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 equation was solved with Matlab using finite differences with a firstorder Euler explicit scheme in time and a centered scheme in space.
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.
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.
During the simulations, when the temperature reached T _{boil}, the heat deposition profile Q was modified at each time step.

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:Let z _{shield} denote the z coordinate where the maximum is reached.$$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)} $$

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 }:R _{ SE } is a parameter of the model to be calibrated.$$\mathcal{H}(t) = \mathcal{B}(t) \oplus \mathcal{S}(R_{SE}) $$

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}:with η _{intercept} as the coefficient characterizing the efficiency of the power absorption, which has to be calibrated.$$P_{\mathcal{H}} = \eta_{\text{intercept}} \, r_{\text{shield}} \, P(z_{\text{shield}}) $$

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 :where the subscript._{foc} refers to the coordinates of the actual focus. Then,$${{\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}}} $$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}}}\)


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.
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
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.
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 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}.
Calibration data
4s pulses: unitary lesion sizes [ mm]  Boiling onset time [ s]  

x  y z  4s pulses  12s pulses  
Measurements  2.55  2.65 6  3  8.1  
Simulations  2.4  2.8 4.2  3  7.9 
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.
Validation data: slope of lesion size vs boiling onset time in mm/s for the 12s pulses. Slopes obtained by orthogonal linear regression on the experimental data are represented by the best fit and the 90 % confidence interval
Slope of lesion size vs boiling onset time [ mm/s]  

Along x  Along y  Along z  
Regression on experimental data  −0.27 [ −0.34; −0.15]  −0.17 [ −0.27; −0.1]  −0.25 [ −0.4; −0.03] 
Simulations  −0.21  −0.27  −0.63 
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.
Validation data: slope of lesion size vs boiling onset time in mm/s for the 4s pulses. Experimental data along z were not considered due to the bias introduced by the slicing method
Slope of lesion size vs boiling onset time [ mm/s]  

Along x  Along y  
Regression on  
experimental data  −0.5  −0.6 
Simulations  −0.6  −0.7 
Growth of the bubble cloud for the 12s pulses
8s pulses and groups of 4s pulses
Experimental conditions related to the influence of neighboring pulses
Median cooling  Spacing along  Spacing along  Overlap  

time (s)  x (mm)  y (mm)  (%)  
4s pulses  22  1.8  1.6  45 
8s pulses  39  4  3.5  19 
12s pulses  48  5  4.3  6.7 
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.
Declarations
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.
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.
Authors’ Affiliations
References
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 Zabolotskaya EA, Khokhlov RV. Quasiplane waves in the nonlinear acoustic of confined beams. Sov Phys Acoust. 1969; 15:35–40.Google Scholar
 Westervelt PJ. Parametric acoustic array. J Acoust Soc Am. 1963; 35:535–7.View ArticleGoogle Scholar
 Hamilton MF, Blackstock DT. Nonlinear acoustics. New York: Acoustical Society of America; 2008.Google Scholar
 Filonenko EA, Khokhlova VA. Effect of acoustic nonlinearity on heating of biological tissue by highintensity focused ultrasound. Acous Phys. 2001; 47:468–75.View ArticleGoogle Scholar
 Hynynen K. The threshold for thermally significant cavitation in dog’s thigh muscle in vivo. Ultrasound Med Biol. 1991; 17:157–69.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedPubMed CentralGoogle Scholar
 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.View ArticleGoogle Scholar
 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.Google Scholar
 Liebler M, Dreyer T, Riedlinger RE. Modeling of interaction between therapeutic ultrasound propagation and cavitation bubbles. Ultrasonics. 2006; 44:319–24.View ArticleGoogle Scholar
 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.Google Scholar
 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.View ArticleGoogle Scholar
 Pennes HH. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948; 1:93–122.PubMedGoogle Scholar
 Wissler EH. Pennes’ 1948 paper revisited. J Appl Physiol. 1998; 85:35–41.PubMedGoogle Scholar
 Nelson DA. Invited editorial on Pennes’ 1948 paper revisited. J Appl Physiol. 1998; 8523:2–3.Google Scholar
 Sapareto SA, Dewey WC. Thermal dose determination in cancer therapy. Int J Radiat Oncol Biol Phys. 1984; 10:787–800.View ArticlePubMedGoogle Scholar
 Shaw A, ter Haar GR, Haller J, Wilkens V. Towards a dosimetric framework for therapeutic ultrasound. Int J Hyperth. 2015; 6736:1–11.Google Scholar
 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.View ArticleGoogle Scholar
 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.Google Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.
 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.Google Scholar
 Treeby BE, Cox BT. kWave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields. J Biomed Opt. 2010; 15:021314–112.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 Goss SA, Frizzell LA, Dunn F. Ultrasonic absorption and attenuation in mammalian tissues. Ultrasound Med Biol. 1979; 5:181–6.View ArticlePubMedGoogle Scholar
 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.View ArticleGoogle Scholar
 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.View ArticleGoogle Scholar
 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.View ArticlePubMedGoogle Scholar
 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.View ArticlePubMedGoogle Scholar