Influence of temperature-dependent acoustic and thermal parameters and nonlinear harmonics on the prediction of thermal lesion under HIFU ablation

1.   Introduction

High intensity focused ultrasound (HIFU) therapy is a non-invasive local treatment technology emerging in recent years. Its principle is to focus ultrasound in the target area and form local high energy density. The tissue temperature can reach above 65 ℃ in an instant and produce coagulation necrosis, so as to achieve the purpose of tumor treatment ^[1,2]. HIFU treatment of tumor is based on the thermal mechanism of tissue lesion, so real-time monitoring of tissue thermal lesion is very important.

Real time and accurate monitoring of thermal lesion caused by HIFU is the key to improving the ablation efficacy and reduce complications. In addition, thermal lesion monitoring provides a quantitative tool for the evaluation after HIFU treatment. The existing thermal lesion monitoring methods of HIFU are mainly non-contact, such as magnetic resonance imaging (MRI) monitoring, ultrasonic imaging monitoring, elastic imaging monitoring and simulation predictions. MRI can not only visualize the thermal lesion after HIFU ablation, but also has high spatial resolution. However, it has some disadvantages such as low temporal resolution of anatomical images, non-transplantable and poor compatibility with HIFU devices ^[3]. Ultrasonic imaging has the advantages of low cost, real-time imaging and simple integration with HIFU system. The integration of ultrasonic imaging probe and HIFU transducer can evaluate the acoustic status in the HIFU radiation propagation path in real time, and monitor the thermal lesion state of tissue. However, the image quality of ultrasound scanner is reduced due to the local inhomogeneity of sound velocity in tissue. The definition of the lesion edge and the high contrast between the lesion area and the surrounding tissue are still the difficulties in ultrasonic imaging monitoring ^[4]. Elastic imaging can be used to observe and characterize the formation of thermal lesion, because the stiffness of protein denaturation and soft tissue coagulation is significantly higher than that of normal tissue, but its wide clinical application is still limited due to its low time efficiency ^[5]. Thermal lesion simulation has the advantages of fast calculation speed ^[6,7] and low cost, and it is easy to be combined with HIFU treatment in practice. At present, many literatures have reported the prediction of thermal lesion caused by HIFU irradiation by simulation ^[8,9]. The simulation results can provide guidance for doctors in the treatment of HIFU, and timely adjust the increase and decrease of heat dose, which is conducive to improve the safety of HIFU in clinical treatment. The spherical beam equation ^[10] (SBE) and Pennes bio-heat transfer equation ^[11] were used to calculate the acoustic pressure field and temperature field, and the thermal lesion area of the tissue can be acquired by calculating the equivalent heat dose with cumulative equivalent minutes (Cumulative Equivalent Minutes, CEM) ^[12], then the thermal lesion area can be obtained. However, acoustic parameters (such as sound velocity, density, B/A, attenuation coefficient) and thermal parameters (such as specific heat capacity and thermal conductivity) related to biological tissue are used as the prior information of prediction in the above equation. The traditional methods usually use constant acoustic and thermal parameters at room temperature or body temperature to predict thermal lesion ^[13,14,15].

Choi et al. ^[13] measured the attenuation coefficient, sound velocity and B/A of porcine liver tissue in vitro at different temperature, and Guntur et al. ^[14] obtained the density, specific heat capacity and thermal conductivity of ex-vivo porcine liver tissue at different temperature. The target physical properties of HIFU heating can effectively change the temperature distribution in biological tissue, that is, the geometry of lesion. Guntur et al. ^[16] theoretically studied the effect of temperature related thermal parameters on the thermal distribution in ex-vivo porcine liver tissue. The thermal lesion predicted by using the constant thermal parameters at 25 ℃ was obviously overestimated, compared with the prediction made by using the thermal parameters measured at temperature which led to thermal denaturation. However, in the actual HIFU treatment, in addition to the change of thermal parameters with the temperature of biological tissue, the acoustic parameters will also change with the change of the temperature of biological tissue ^[17,18]. Therefore, it is necessary to study the effects of temperature-dependent acoustic and thermal parameters on the thermal lesion of biological tissue during HIFU treatment.

Liu et al. found that the nonlinear parameter B/A cannot only reflect the change of liquid temperature, but also reflect the temperature change of tissue more sensitively than attenuation coefficient and sound velocity ^[19,20]. The nonlinear parameter B/A is related to the density and sound velocity of tissue, while the generation of harmonic is related to the nonlinear parameter B/A. The harmonic energy caused by ultrasonic nonlinearity is more easily absorbed by tissue, which can produce tissue thermal lesion. Harmonics can also be used to reflect the degree of tissue lesion, which can provide a reference for the prediction of tissue thermal lesion.

Based on the existing experimental data, we have used the data fitting method to obtain the temperature-dependent acoustic and thermal characteristic parameters of porcine liver in vitro, and we have also established a new acoustic-thermal coupled model by combining SBE acoustic equation and Pennes bio-heat transfer equation. The temperature-dependent acoustic and thermal characteristic parameters of porcine liver were used to predict thermal lesion, and the relationship between harmonic amplitude ratio (P₂/P₁) and thermal lesion was analyzed.

2.   Principles and methods

The focused ultrasound transducer has good axial symmetry. According to the study of Tjøtta et al., it was found that the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation was only applicable to the nonlinear acoustic field with the aperture less than 16° of the transducer ^[21,22]. Kamakura et al. established the SBE equation of large open angle ultrasonic transducer based on Westervelt equation by using ellipsoid coordinate system, which includes nonlinear, diffraction and absorption effects, and it can be used to calculate the acoustic field of large angle transducer ^[23]. Figure 1 shows a flat spherical coordinate system in which variables $(\sigma, \eta, \varphi)$ are related to rectangular coordinates:

where $- \infty < \sigma < \infty$, ${\rm{0}} \leqslant \eta \leqslant {\rm{1}}$, ${\rm{0}} \leqslant \phi \leqslant {\rm{2}}\pi$.

where $2b$ is the length between focal points in a flat spherical coordinate system. Combined with $\sigma$, $\eta$ and $\phi$, the scalar acoustic equation can be expressed as follows:

where p is the sound pressure, ${c_0}$ is the velocity of sound, ${\nabla ^2} = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2} + {\partial ^2}/\partial {z^2}$, and $t$ is the time. For Eq (3), it is assumed that the sound field is axisymmetric, so the equation does not contain the variable $\varphi$.

Figure 2 shows the structure of a concave sound source with a circular hole. $d$ is the geometric focal length, $a$ is the outer radius of the circular hole, and ${\alpha _0}$ is the half aperture angle. The relationship between ${\alpha _0}$ and $a$, $d$ is as follows:

The focused sound field is divided into two regions, $\sigma < {\sigma _0} < 0$ and ${\sigma _0} \leqslant \sigma$, where $\sigma = {\sigma _0}( = const < 0)$ represents a specific transition position. In Figure 2, on the left is the spherical wave approximation region, and the right is the plane wave approximation region. For the calculation of spherical wave and plane wave, please refer to reference ^[23]. The SBE equation is used to simulate the acoustic field of HIFU.

The Pennes bio-heat transfer equation is used to simulate the temperature field, and it can be expressed as ^[11]:

where ${C_{\rm{t}}}$ is the specific heat capacity, ${k_t}$ is the thermal conductivity of organism, $T$ is the medium temperature. The heat loss caused by blood flow is not considered in the equation, and ${Q_v}$ is the heat accumulation caused by sound field.

The time domain algorithm of finite difference operator is used to solve the heat conduction equation of organism to obtain the temperature distribution in the biological tissue ^[24]. The thermal lesion of tissue is related to the denaturation of protein, and the lesion area can be judged by the state of denaturation and ablation, which is related to the heat dose. The heat dose is used to characterize the thermal lesion result which is equivalent to the equivalent heating time at 43 ℃. When T ≥ 43 ℃, R = 0.5; T < 43 ℃, R = 0.25. It is usually used as the equivalent threshold for the formation of visible thermal lesion in biological tissue at 43 ℃ for 240 minutes.

The external radius a = 2.5 cm, inner radius a₁ = 1 cm, focal length d = 6 cm, center frequency f = 1 MHz. The HIFU irradiation signal is sinusoidal pulse cluster signal, the duty cycle is D = 10%, and the distance of ultrasonic propagation in water is d₁ = 5 cm. Assuming that porcine liver is a homogeneous medium, the thickness of porcine liver for simulation is h = 2 cm, width w = 4 cm, and the geometric model of HIFU irradiated porcine liver is shown in Figure 3. The radiation sound pressure on the surface of the transducer can be adjusted, and the radial and axial step sizes are 8.9 × 10⁻³ mm and 6.8 × 10⁻³ mm, respectively. The time steps of sound field and temperature field are 4.2 × 10⁻⁸ s and 2 × 10⁻⁴ s, respectively, and the nonlinear harmonic order is 128. The presented method was evaluated on an Intel(R) Core (TM) i7-9750H CPU @ 2.60 GHz with 8.0 GB RAM PC using Windows 10 operating system, we simulated a 7.9 s HIFU ablation, and the total CPU computation time was 12625.7 s.

3.   Acoustic and thermal characteristic parameters of biological tissue related to temperature

There is a strong correlation between the acoustic and thermal parameters of porcine liver tissue and temperature. Choi and Guntur et al. measured the attenuation coefficient, sound velocity, B/A, density, specific heat capacity and thermal conductivity of porcine liver in vitro ^[13,14]. Because the treatment time of HIFU is short and the thermal lesion of porcine liver irradiated by HIFU is mainly studied, the effect of blood perfusion rate is ignored here. Among them, the current literature only gives the nonlinear parameter B/A value not exceeding 75 ℃, while the nonlinear parameter B/A above 75 ℃ change very little with the increase of temperature ^[13,25]. In order to simplify the problem analysis, the nonlinear parameter B/A above 75 ℃ is still replaced by the nonlinear parameter B/A corresponding to 75 ℃. Combining the known acoustic and thermal experimental data of porcine liver and using the cftool toolbox of MATLAB, the attenuation coefficient, B/A, density, sound velocity, specific heat capacity and thermal conductivity of porcine liver can be obtained with temperature-dependent expressions ^[13,14].

The corresponding graph of the above polynomial is shown in Figure 4. The acoustic and thermal parameters corresponding to different temperature are used in the following calculation according to the above expressions. Constant acoustic and thermal parameters of porcine liver and water at 30 ℃ are shown in Table 1.

4.   Thermal lesion prediction based on temperature-dependent acoustic and thermal characteristic parameters of biological tissue

4.1.   Acoustic and thermal coupled model

During HIFU irradiation, the temperature of porcine liver tissue increased rapidly, and its acoustic and thermal characteristic parameters also changed ^[26,27,28]. The acoustic and thermal characteristic parameters of porcine liver tissue need to be updated in real time according to the temperature, and the updated acoustic and thermal characteristic parameters are fed back to the calculation of sound field and temperature field, and the change of porcine liver temperature in turn affects the acoustic and thermal characteristic parameters. The sound field and temperature field are coupled by heat accumulation Q_v. In the simulation, the SBE equation is used to calculate the sound field and the Pennes equation is used to calculate the temperature field. The relationship between the acoustic and thermal characteristic parameters and the temperature of porcine liver is used to update the acoustic and thermal characteristic parameters. For biological tissue, the nonlinear coefficient β is related to the nonlinear parameter B/A, and B/A is closely related to the acoustic and thermal characteristic parameters of the tissue. Therefore, the harmonic amplitude is related to B/A, the density and sound velocity of the tissue, while the acoustic and thermal characteristic parameters such as density and sound velocity are related to temperature.

The relationship between the nonlinear harmonic amplitude ratio (P₂/P₁) at the focal point and the thermal lesion of porcine liver was calculated and analyzed dynamically through the coupling relationship between sound field and temperature field. Figure 5 shows the calculation of P₂/P₁ and the prediction of thermal lesion using temperature-dependent acoustic and thermal characteristic parameters under the coupling of sound field and temperature field.

In the process of numerical simulation, the acoustic and thermal parameters are updated every 0.1 s, the sound field is solved by the finite difference frequency domain algorithm, and the temperature field is solved by the finite difference operator separation time domain algorithm.

4.2.   Results and discussion

Figure 6 shows the comparison of predicted thermal lesion using the acoustic and thermal characteristic parameters of ex-vivo porcine liver tissue at constant temperature (30 ℃) and temperature-dependent acoustic and thermal characteristic parameters under the conditions of P₀ of 0.23 MPa and 0.25 MPa and different irradiation time, respectively.

From Figure 6 (a, b), it can be found that under a certain sound pressure, the thermal lesion increased with the increase of irradiation time. At the same sound pressure and irradiation time, the thermal lesion predicted by temperature-dependent acoustic and thermal characteristic parameters was larger than that predicted by constant characteristic parameters.

Figure 7 shows the corresponding relationship between P₂/P₁ at the focus and the irradiation time when P₀ = 0.23 MPa and P₀ = 0.25 MPa, respectively.

It can be seen from Figure 7 (a, b) that the variation trend of P₂/P₁ with irradiation time was basically the same. The turning point of P₂/P₁ appeared in Figure 7 (a) and Figure 7 (b) at irradiation time t = 1.4 s and t = 1.1 s, respectively. P₂/P₁ increased first and then decreased with the increase of irradiation time.

Figure 8 shows the relationship between the calculated P₂/P₁ and the predicted thermal lesion using temperature-dependent acoustic and thermal characteristic parameters.

It can be seen from Figure 8 (a, b) that under different sound pressure conditions, at the beginning of irradiation, the P₂/P₁ at the focus position increased from small to large at first, and no visible thermal lesion appeared in porcine liver tissue. Then, the P₂/P₁ decreased gradually from large to small, but the thermal lesion of porcine liver tissue increased gradually.

Figures 7 and 8 show the two-dimensional (2D) relationship between P₂/P₁ and irradiation time, the lesion area and P₂/P₁, respectively. In order to better describe the relationship among thermal lesion area, irradiation time and P₂/P₁ of porcine liver, their three-dimensional (3D) relationship can be obtained by using interpolant function of cftool toolbox of MATLAB, as shown in Figure 9.

Figure 9 shows the relationship among thermal lesion area, P₂/P₁ and irradiation time. Although P₂/P₁ increased gradually with the increase of irradiation time when P₀ = 0.23 MPa and P₀ = 0.25 MPa, respectively, there was no visible thermal lesion in porcine liver tissue when the irradiation time was less than 1.4 s and 1.1 s, respectively, because the thermal dose did not reach the threshold of lesion. Then, P₂/P₁ decreased gradually when the irradiation time was longer than 1.4 s and 1.1 s, respectively, but the thermal dose in porcine liver reached the threshold of lesion. In addition, with the increase of irradiation time, more and more heat accumulated in porcine liver, so the thermal lesion of porcine liver tissue increased gradually.

5.   Conclusions

The effects of temperature-dependent acoustic and thermal characteristic parameters and nonlinear harmonics on thermal lesion prediction were studied, and a new acoustic-thermal coupled model based on nonlinear SBE acoustic equation and Pennes bio-heat transfer equation was proposed. The thermal lesion predicted by temperature-dependent acoustic and thermal characteristic parameters was compared with that predicted by constant temperature acoustic and thermal characteristic parameters, and it was found that the thermal lesion predicted by the temperature-dependent acoustic and thermal characteristic parameters was larger than that predicted by traditional method, and the harmonic amplitude ratio increased gradually at first, but there was no thermal lesion in porcine liver. Then the harmonic amplitude ratio decreased gradually, but the thermal lesion in porcine liver increased gradually, which may provide reference for clinical safety treatment of HIFU.

Future research work will focus on two aspects: First, experimental conditions will be created and experimental studies will be carried out to back up the simulation results. Second, besides thermal analysis, the effect of mechanical deformation on biological tissue lesion will be investigated.

Acknowledgments

This work was supported by the Hunan Provincial Natural Science Foundation of China under grant No. 2018JJ3557, the Scientific Research Fund of Hunan Provincial Education Department under grant No. 17B025, and the Key Program of Changsha Normal University under grant No. 2019XJZK12. The authors sincerely thank the anonymous reviewers for their helpful comments and suggestions.

Conflict of interest

The authors declare no conflict of interest in this paper.