Numerical analysis of thermal response of tissues subjected to high intensity focused ultrasound

Abstract

The present work is concerned with the numerical investigation of the thermal response of tissue-mimicking biological phantom(s) subjected to high intensity focused ultrasound (HIFU). Simulations have been performed on the 3-dimensional physical domain for two-layered as well as multi-layered medium consisting of water and liver tissue. Local pressure distribution within the body of the phantom has been calculated by solving the complete full-wave nonlinear form of Westervelt equation. The solution of the pressure wave equation has been coupled with Pennes bioheat transfer equation to determine the full field temperature distribution. Results in the form of pressure fields, temperature distributions and the corresponding thermal dosage in the targeted region of the tissue domain have been presented. Magnitudes of the maximum pressure (and hence the resultant temperature levels) in the focal region as obtained using the nonlinear form of Westervelt equation are found to be significantly higher than that determined based on the linear form of the equation. Compared to water, wherein the acoustic intensity is maximum, the addition of sub-layers of skin, fat, and muscle into water resulted in the reduction of the peak intensity and also shifted the intensity profiles along the direction of propagation of the acoustic waves. However, addition of liver tissue into water led to the shifting of intensity profile in the opposite direction i.e., towards the transducer. The results further reveal that due to the dependence of attenuation coefficient on the source frequency, the temperature at the focal region increases with an increase in the transducer frequency. The work is further extended from single lesion to multiple lesion generation through controlled movement of the transducer and the resultant transient full field temperature distribution has been presented. The concerned observations highlight the need of optimizing the thermal energy for each lesion, the inter spatial distance between different lesions and the delay time so as to ensure minimal thermal damage to the surrounding healthy cells as well as to reduce the total treatment duration.

Keywords:

• High intensity focused ultrasound (HIFU)
• nonlinear Westervelt equation
• Bio-heat transfer
• thermal ablation
• multiple lesion generation

1. Introduction

High intensity focused ultrasound (HIFU) holds its importance in therapeutic applications because of its potential to achieve selective destruction of the targeted region in a calculated and controlled manner. This technique is one of the minimal invasive/non-invasive methods, which generates heat in the body when the targeted delivery of ultrasound waves causes a sudden change in temperature (more than 56 °C) of the local tissue environment [1]. The central idea behind HIFU-based thermal therapy approach is the conversion of the acoustic energy into its thermal form. In the localized region, wherein the high-intensity ultrasound waves get focused, the absorption of the acoustic energy leads to an increase in the temperature and depending upon the operating parameters, temperature levels required for complete necrosis of the embedded abnormal cells may be achieved. In this context, one of the major challenges is the selective destruction of the abnormal cells (malignant and/or benign) without any unwanted thermal damage to the surrounding healthy/normal cells. Thus, it becomes important to develop a detailed understanding of the whole field temperature distribution within the body of the HIFU- subjected tissue sample.

Generally, HIFU transducers are excited at a single frequency that lies in the range of 0.5–3.5 MHz depending on their usage in clinical application [2]. If the sound wave propagates linearly through the tissue volume, the achievable heating rates depend on the factors such as incident ultrasound intensity, acoustical properties of the medium (e.g., attenuation or absorption coefficient) etc. [3]. Under operating conditions wherein HIFU transducers produce relatively higher levels of intensities and also due to the dependence of absorption coefficient on the working frequencies, the propagation of the acoustic waves in the medium becomes a nonlinear phenomenon. The nonlinear effects lead to the generation of higher harmonics, which in turn, lead to higher attenuation and heat deposition in the tissue, primarily in the focal region [4,5]. The importance of such nonlinear effects as one of the possible means to enhance the ultrasound-based thermal effects in a controlled manner has been realized by a range of researchers in the past [6]. In this context, various attempts have been made to couple the solution of the acoustic model with an appropriate thermal model so as to determine the resultant temperature distributions. However, the literature shows that such studies have primarily been limited to the linear form of acoustic waves propagation model [7–10]. The studies concerned with the development of analytical and/or numerical models for understanding the phenomenon of nonlinear propagation of ultrasound waves in the tissue medium coupled with an appropriate thermal model are still limited. In addition to nonlinear effects, the loss of energy depends on the nature of the physical process, which includes viscous losses, thermal conduction and various forms of molecular relaxation processes [11]. In this context, the Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation, which is a nonlinear parabolic equation, has often been employed to predict the pressure distribution in soft tissue medium obtained as a result of directional sound beams [12]. The potential advantages of the KZK equation-based approach include significantly reduced simulation time for obtaining the desired solution with relatively good accuracies due to the retarded time co-ordinate [13]. As a result, the KZK equation has found its applicability in wide range of applications in recent times. However, the application of KZK equation is limited by the constraint that it does not take into account the effects of reflection and scattering and also it can only be applied for transducers with small aperture angles [14]. Moreover, this equation cannot be applied to more realistic physical models that take into account the multi-layered nature of the biological tissues samples [15].

With this background, in attempt to couple the nonlinear effects of acoustic wave propagation on the resultant thermal field in the tissue domain (single as well as multi-layered), the present work employs the full wave nonlinear form of Westervelt equation to determine the acoustic pressure distribution in the multi-layered tissue medium. The approach becomes important as the governing equation takes into account the effects of diffraction, absorption as well as nonlinear propagation. In order to control the reflection from the opposite boundary of the domain, the concept of Perfectly Matched Layer (PML) has been implemented at the boundary surface. The results of Westervelt equation (in the form of pressure distribution) have been coupled with the Pennes bio-heat transfer equation to obtain the full field temperature distribution calculated in the tissue medium. It is also to be mentioned that a majority of studies reported in the literature are primarily limited to the generation of single lesion only and the associated thermal effects. In this direction, as an advancement, the present work is extended to the generation of multiple lesions by controlled movement of the transducer.

2. Physical and numerical models

2.1. Propagation of acoustic wave into the tissue medium

Propagation of ultrasound waves into the tissue medium causes its particles to oscillate back and forth and results in a series of compression and rarefaction (alternating cycles of increased and reduced pressures respectively) of the medium [16]. As the soft biological tissues are thermo-viscous in nature, increase in pressure increases the local temperature of the tissue domain leading to an increase in the local speed of sound in the medium. As a result, the wave during the high-pressure phase of oscillation travels faster than the one during the lower pressure phase. Thus, the traveling waves, which initially carry energy at the fundamental frequency, slowly move to the higher harmonics. This results in the faster travel of the peaks of the wave than the troughs, and the wavefront gets as steepened as a sawtooth wave.

Nonlinear acoustics play an important role in distortion of the waveform during transmission and can be described by a parameter, B/A, which is basically the ratio of the second and the first terms in the Taylor series expansion of pressure as a function of density. The time and distance of propagation of the traveling wave also influence the extent of medium distortion [17]. The nonlinear effects lead to the generation of harmonics that result in high attenuation and heat deposition in the tissue, mainly in the focal region. Full wave Westervelt equation, which is based on the effects of diffraction, absorption, and nonlinear propagation, is given by [18]: (1) $${\nabla }^{2}p=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{\delta }{c_0^4}\frac{{\partial }^{3}p}{\partial t^3}-\frac{\beta }{{\rho }_0{c}_0^4}\frac{{\partial }^2{p}^2}{\partial t^2}$$(1)

In Equation (1)(1) $${\nabla }^{2}p=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{\delta }{c_0^4}\frac{{\partial }^{3}p}{\partial t^3}-\frac{\beta }{{\rho }_0{c}_0^4}\frac{{\partial }^2{p}^2}{\partial t^2}$$(1) above, p denotes the local sound pressure, $\beta (=1+B/2A)$ denotes the coefficient of nonlinearity and δ denotes the diffusivity of sound resulting from thermo-viscous fluid viscosity and heat conduction and can be expressed as [19]: (2) $$\delta =\frac{2c_0^3\alpha }{{\omega }^2}$$(2) where, α is the acoustic absorption coefficient and ω is the source angular frequency. For most tissues, the absorption coefficient is related to the ultrasound frequency via a power law of the form: (3) $$\alpha \left(f\right)=a{f}^b$$(3)

Here f is the ultrasound source frequency, a and b are the tissue-specific constants [20]. From the above power law, it is seen that the absorption coefficient increases with an increase in the ultrasound excitation frequency. Higher absorption coefficient implies higher heat deposition. On the other hand, higher attenuation results in the reduced penetration depth. The transmission of ultrasound energy within two tissue layers depends particularly on the density, speed of sound and the thickness of each layer. In Equation (1)(1) $${\nabla }^{2}p=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{\delta }{c_0^4}\frac{{\partial }^{3}p}{\partial t^3}-\frac{\beta }{{\rho }_0{c}_0^4}\frac{{\partial }^2{p}^2}{\partial t^2}$$(1) , the first two terms in RHS of the equation describe the linear lossless wave which gets transmitted at a small-signal sound speed. The third term represents attenuation due to thermal conductivity and fluid viscosity. The last term describes the acoustic nonlinearity, which explains the distortion of a wave due to the nonlinear effects and affects the thermal and mechanical changes inside the tissue medium [21].

2.2. Thermal model for temperature distribution

In order to determine the temperature distribution in the tissue domain, solution of Equation (1)(1) $${\nabla }^{2}p=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{\delta }{c_0^4}\frac{{\partial }^{3}p}{\partial t^3}-\frac{\beta }{{\rho }_0{c}_0^4}\frac{{\partial }^2{p}^2}{\partial t^2}$$(1) , which mathematically models the phenomenon of acoustic wave propagation in the physical domain, needs to be coupled with an appropriate bio-heat transfer model. In the present work, the solution of the full wave Westervelt equation (Equation (1)(1) $${\nabla }^{2}p=\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}-\frac{\delta }{c_0^4}\frac{{\partial }^{3}p}{\partial t^3}-\frac{\beta }{{\rho }_0{c}_0^4}\frac{{\partial }^2{p}^2}{\partial t^2}$$(1) ) has been coupled with the Pennes bio-heat transfer equation [22] (Equation (4)(4) $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k{\nabla }^{2}T\hspace{0.17em}-\hspace{0.17em}{\omega }_b{C}_b(T\hspace{0.17em}-\hspace{0.17em}{T}_{art})+Q$$(4) , given below): (4) $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k{\nabla }^{2}T\hspace{0.17em}-\hspace{0.17em}{\omega }_b{C}_b(T\hspace{0.17em}-\hspace{0.17em}{T}_{art})+Q$$(4)

In Equation (4)(4) $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k{\nabla }^{2}T\hspace{0.17em}-\hspace{0.17em}{\omega }_b{C}_b(T\hspace{0.17em}-\hspace{0.17em}{T}_{art})+Q$$(4) , ρ_t, C_t, and k are density, specific heat and the thermal conductivity of the tissue, respectively. ρ_b and C_b are respectively the density and specific heat of blood and ${\omega }_{b\mathrm{ }}$is the blood perfusion rate. T and T_art respectively represent the tissue and arterial blood temperatures. In Equation (4)(4) $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k{\nabla }^{2}T\hspace{0.17em}-\hspace{0.17em}{\omega }_b{C}_b(T\hspace{0.17em}-\hspace{0.17em}{T}_{art})+Q$$(4) above, Q represents the volumetric heat generation due to the absorption of acoustic intensity in the tissue domain and can be calculated by [23,24]: (5) $$Q(r,z)=\frac{2\alpha }{{\rho }_0 c_0{\omega }^2}\frac{1}{T}{\underset{0}{\overset{T}{\int }}\left(\frac{\partial p}{\partial t}\right)}^{2}dt$$(5)

For the ablation of the diseased tissue, heat is generated in the body through the absorption of ultrasound energy. The extent of thermal damage of the tissue depends on the duration of exposure and the final temperature levels reached in the process [25]. Sapareto and Dewey proposed an empirical relationship for the calculation of the thermal dose from temperature-time history. This equation uses the numerical integration of temperature with the time for determining the equivalent thermal dose at the reference temperature [26,27]. (6) $$TD\left(t\right)=\underset{0}{\overset{t}{\int }}{R}^{(T_{ref}-T(t^{\prime }\left)\right)}d{t}^{\prime }$$(6) where, T_ref is the reference temperature. In general, the reference temperature is chosen to be equal to 43 °C as this is the standard temperature that is commonly used in various experimental thermal data. The value of R (iso-dose constant) has been chosen based on the animal experiments performed in the literature [26]. (7) $$R=\{\begin{array}{c}0.25,\hspace{1em}T\left(t\right)\hspace{0.17em}<\hspace{0.17em}43{ }^{\mathrm{°}}C\\ 0.5,\hspace{1em} T\left(t\right)\ge 43{ }^{\mathrm{°}}C\end{array}$$(7)

For nearly complete necrosis of the biological tissues, TD needs to be in the range of 25 min to 240 min [10]. Ultimately, the thermal dose calculation is required to predict the dimensions of lesions generated in the HIFU therapy.

3. Numerical methodology

3.1. Computational domain

Numerical simulations have been performed on a 3 D Cartesian domain (symmetric about the z-axis) with a volume of 10 cm³. Schematic drawing of the physical domain considered in the present work has been shown in Figure 1(a). As shown, the direction of propagation of the acoustic waves (originating from the transducer) has been set along the z-direction. The space between tissue and transducer is filled by water medium. Based on the number of layers present in the tissue sample, simulations have been performed on two-layered medium (consisting of water and liver tissues) and multi-layered medium (wherein layers of water, skin, fat, muscle, and liver tissue have been considered). These two domains, with appropriate dimensions, have schematically been shown in Figure 1(b.i, b.ii). In order to control the reflection of acoustic waves from the interior of the boundary surface, a perfectly matched layer (PML) of 3 mm thickness has been added in the domain for calculating the acoustic pressure.

Figure 1. (a) Schematic of the HIFU-subjected tissue phantom geometry for simulations. (b) Schematics of ZX plane (axis-symmetric along z-axis) around y = 0: (b.i) two-layered medium, (b.ii) multi-layered medium and (b.iii) multi-layered with inhomogeneity. (c) Boundary conditions for pressure (c.i) and temperature (c.ii). (All dimensions are in cm).

In addition to homogeneous samples, simulations have also been performed on non-homogeneous tissue phantoms wherein an inhomogeneity has been considered to be embedded inside an otherwise homogeneous multi-layered tissue phantom. This configuration has been schematically depicted in Figure 1(b.iii). As shown, an optical inhomogeneity of volume 1 cm³, has been considered to be present in the layer of the liver tissues. In physical terms, the optical inhomogeneity numerically simulates the presence of abnormal cells (malignant and/or benign) present in an otherwise homogenous domain. Thus, the thermo-physical properties of the embedded inhomogeneity have been kept identical to that of the cancerous cells, as documented in the literature [10].

3.2. Initial and boundary conditions

The present section discusses the initial and boundary conditions that the computational domain has been subjected to for solving the governing equations presented in the previous section. In any given physical domain considered, the acoustic pressure distribution has been determined through the solution of the Westervelt equation formulated in 3 D Cartesian space coordinate. The input signal from a focused circular transducer working on a single frequency has been modeled using planar pressure source approximation [28, 29]. Pressure produced by the transducer (placed at z = 0) has been calculated using the following expression: (8) $$p(r,z=0,t)=\{\begin{array}{c}\frac{p_0}{\sqrt{1+r^2/d^2}}F\left(t+\frac{d}{c_0}\left(\sqrt{1+r^2/d^2}-1\right)\right)\hspace{1em}for\hspace{0.5em}r\le d \hspace{0.17em}\text{tan}\hspace{0.17em}\alpha \\ \hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} for\hspace{0.5em}r\hspace{0.17em}>\hspace{0.17em}d \hspace{0.17em}\text{tan}\hspace{0.17em}\alpha \end{array}$$(8) where p₀ is the source pressure amplitude and F(t) is taken as sin(ωt). As mentioned, for Cartesian space coordinate, $r=\sqrt{x^2+y^2}$ and $\text{tan}\hspace{0.17em}\alpha =a/\sqrt{1-\frac{a^2}{d^2}}$, where α denotes the aperture angle, a is the radius of the transducer and d denotes the focal length. Table 1 summarizes the parameters used for the source transducer used in the present work. As presented in the table, the transducer frequency has been set as 1.0 MHz. It is worth mentioning here that frequencies in the range of 0.8–1.6 MHz are commonly employed for liver cancer applications i.e., for the treatment of abnormal cells embedded relatively deep in the biological samples [3]. Furthermore, the source pressure oscillates sinusoidally in the range of +0.5 MPa to −0.5 MPa (peak to peak). The choice of the range of source pressure employed has been made to ensure that the maximum temperatures achieved in the focal region do not exceed 100 °C and hence to avoid any possibility of cavitation effects.

Table 1. Transducer Parameters.

Symbol	Property	Value
$a$	Radius of transducer (Aperture = $2a$)	3.5 cm
$d$	Focal length (Transducer radius of curvature)	7.0 cm
$f_0$	Transducer frequency	1.0 MHz

Figure 1(c) shows the schematic representation of the computational domain depicting the initial and boundary conditions used for numerically solving the Westervelt equation. As shown, initially, throughout the body of the physical domain, the pressure and its derivative with respect to time have been set at zero (i.e.,$\mathrm{ }p\left(x,y,z,t\right)=0,\mathrm{ }\partial p/\partial t=0$). At the transducer surface, i.e., at z = 0, the pressure distribution$\mathrm{ }p\left(r,z=0,t\right)\mathrm{ }$has been obtained using Equation (8)(8) $$p(r,z=0,t)=\{\begin{array}{c}\frac{p_0}{\sqrt{1+r^2/d^2}}F\left(t+\frac{d}{c_0}\left(\sqrt{1+r^2/d^2}-1\right)\right)\hspace{1em}for\hspace{0.5em}r\le d \hspace{0.17em}\text{tan}\hspace{0.17em}\alpha \\ \hspace{1em}0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} for\hspace{0.5em}r\hspace{0.17em}>\hspace{0.17em}d \hspace{0.17em}\text{tan}\hspace{0.17em}\alpha \end{array}$$(8) . The other boundary conditions have been shown in Figure 1(c.i). As the computational domain is considered axis-symmetric around the z-axis, the simulations have only been performed for one quarter of the total size of the computational domain (i.e.,$\mathrm{ }0\le x\le 5\mathrm{ }\mathrm{c}\mathrm{m}$,$\mathrm{ }0\le y\le 5\mathrm{ }\mathrm{c}\mathrm{m}$ and $0\le z\le 10\mathrm{ }\mathrm{c}\mathrm{m}$) to reduce the computational time.

In order to control the reflection from the opposite boundary of the domain, the concept of Perfectly Matched Layer (PML) has been implemented at the boundary surface. This condition of PML is utilized at the end of the domain, thus $p\left(x,y,z=10\mathrm{ }\mathrm{c}\mathrm{m},t\right)\mathrm{ }$has been taken to be equal to zero, as the layer ultimately decays the solution to almost zero. It is to be mentioned here that the PML is a layer of artificial absorbing material that is placed adjacent to the edges of the grid to change the oscillating solution into an exponentially decaying solution without any reflection [30]. The following steps have been followed for deriving the wave propagation into the PML:

1. Fourier transform of the 2^nd order linear wave equation.

2. PML transformation of $\frac{\partial }{\partial z}\to \frac{1}{1+i\frac{{\sigma }_z}{\omega }}\frac{\partial }{\partial z}$

3. Introduce new auxiliary variables$\mathrm{ }{D}_1$ and${\mathrm{ }D}_2$.

4. Truncate the domain.

5. Inverse transform to convert the Fourier domain back into the time domain with PML absorbing boundaries in the z-direction.

In equation form, the PML can be expressed as [31] (9) $$\frac{1}{c_0^2}\frac{{\partial }^{2}p}{\partial t^2}=\frac{{\partial }^{2}p}{\partial x^2}+\frac{{\partial }^{2}p}{\partial y^2}+\frac{\partial D_2}{\partial t}$$(9) (10) $$\frac{\partial D_1}{\partial t}+{\sigma }_z{D}_1=\frac{\partial p}{\partial z}$$(10) (11) $$\frac{\partial D_2}{\partial t}+{\sigma }_z{D}_2=\frac{{\partial }^2{D}_1}{\partial z\partial t}$$(11) Where, $D_1$ and $D_2$ are the auxiliary variables. The conductivity profile of the PML [31] is given by (12) $${\sigma }_z\left(z\right)=-{\left(\frac{z}{Z_p}\right)}^m\frac{c_0(m+1)\hspace{0.17em}\text{ln}\left(R_0\right)}{2Z_p}$$(12) where,${\mathrm{ }Z}_p$ is the PML thickness, m is the order and $R_0$ is the reflection of the PML at normal incidence. These parameters are optimized and taken for 100 matching layers and thus m = 3 and $R_0={10}^{-7}\mathrm{ }$have been considered in the present simulations [31]. The number of matching layers depends on the requirement of the extent of the pulse absorption.

Similarly, the initial temperature of the computational domain has been taken as the normal body temperature, i.e., 37 °C and other boundary conditions have been shown in Figure 1(c.ii). Acoustic and the thermophysical properties for various computational domains considered in the present study have been summarized in Table 2 [14].

Table 2. Acoustic and thermo-physical properties.

Symbol	Properties	Unit	Water	Skin	Fat	Muscle	Liver tissue	Blood
$c_0$	Speed of sound	m/s	1520	1540	1430	1560	1570	1540
${\rho }_0$	Density of material	kg/m^3	1000	1100	910	1050	1055	1060
${\alpha }_0$	Attenuation coefficient	Np/m	0.026	16	7	9	8.1	1.5
$\beta$	Nonlinearity coeff.		3.5	4.9	6.5	4.5	4.5	4.0
k	Thermal conductivity	W/m °C	0.6	0.28	0.24	0.51	0.51	0.53
Cp	Heat capacity	J/kg °C	4200	3500	3800	3700	3600	3770
ωb	Volumetric Blood perfusion rate	kg/m^3s	0	0	0	0.5	10	–

3.3. Grid independence

Governing equations (defined in 3 D Cartesian coordinate system) along with the initial and boundary conditions have been discretized using explicit finite difference scheme in time domain [32,33]. A detailed grid independence study has been performed before converging to an optimum grid size employed for the detailed thermal investigation of HIFU-subjected biological tissue phantoms. Results of the grid independence study have briefly been discussed here. Figures 2(a) and 2(b) show the variations of pressure with respect to the axial direction (z) for different grid sizes Δz and Δx respectively expressed in terms of the sound wavelength (λ = c/f, where f is the driving frequency). In the present study, an optimized grid size of Δz = 0.1λ and Δx = Δy = 0.2λ has been considered. Also, the discretization step size in time for solving the Westervelt equation (for pressure field) has been set as Δt = 0.01/f = 10⁻⁸s while for temperature field distribution, the time step of Δt = 0.001 s has been considered. The difference in the two time scales employed is to be attributed to the multi-physics numerical modeling wherein the acoustic wave takes about microseconds (50–60 µs in this study) to reach the targeted area. On the other hand, the temperature changes can only be observed in seconds (limited by the thermal diffusivity of the medium considered). For PML layer, all the grid sizes remain same, except Δz that has been taken as 0.1λ/5, so that a total number of 100 layers can be introduced within PML thickness and hence can have better control on reflection. The complete code has been written in C++ programming language on a computer with the following specifications: RAM 16.0 GB, 64-bit operating system processor Intel (R) Core (TM) i7-2600K CPU @ 3.40 GHz.

Figure 2. Pressure variation with respect to axial direction (along the direction of propagation of the ultrasound waves) for different grid sizes in space (a) for different dz (b) for different dx, to optimise the discretized step size in space.

4. Results and discussion

4.1. Verification of computational method

The numerical methodology developed in the present work has been verified against the known available solutions available in the literature. For the verification studies, the linear acoustic pressure field obtained through the numerical solution of the developed numerical model (based on 3D Westervelt equation) has been compared with the results obtained by the analytical solution of Rayleigh integral for planer source approximation [28]. Figure 3 illustrates the time (a) and space variation (b) of acoustic pressure in the water medium for the source pressure of p_o = 0.5 MPa and transducer frequency f = 1 MHz. It can be clearly seen that both the pressure variations show a good agreement between the present numerical simulations and the analytical solution. As expected, the peak of acoustic pressure is highest at the focal spot due to the focused nature of the HIFU (z = 7.0 cm, Figure 3(b)).

Figure 3. Comparison of the numerical solution of Westervelt equation as obtained in the present work with that of the analytical Rayleigh integral solution given by Hamilton [28]; (a) time and (b) space variation of acoustic pressure in water medium.

4.2. Thermal response of the HIFU- subjected two-layered homogeneous medium

The present section discusses the results of the numerical simulations performed for understanding the thermal response of two-layered homogeneous medium (physical domain without any presence of inhomogeneity) that has been subjected to high intensity focused ultrasound. The homogeneous two-layered tissue phantom comprises of the layers of water and liver tissue. The schematic representation of the physical domain has earlier been shown in Figure 1(b.i) (Domain 1). Here, the effects of nonlinearity and attenuation have also been taken into account. Figure 4 shows the axial distribution (z-direction) of peak compressional and rarefactional pressures as obtained through the solution of linear and nonlinear forms of Westervelt equation. It is to be seen that due to the nonlinear nature of propagation of sound waves, the magnitude of nonlinear peak positive pressure (P₊=16.08 MPa) is significantly higher than the corresponding peak negative pressure (P- = 10.36 MPa) in the focal region (z = 7 cm). On the other hand, linear Westervelt equation predicts identical values for the compression and rarefaction pressure distributions (P = 13.03 MPa). These results show that the nonlinear effects are concentrated primarily over a small region spread on either side of the focal spot (focal area) of the soft tissue.

Figure 4. Nonlinear peak positive and peak negative pressure distribution compared with linear pressure distribution (a) with respect to time and (b) along the axial distance, (c) volumetric heat generation and (d) temperature distribution along the axial distance.

A comparison of the profiles shown in Figure 4(a,b) brings out the clear differences in the linear and nonlinear forms of the solutions of full wave Westervelt equation. The observed differences are predominately concentrated in the close vicinity of the focal region and the profiles nearly overlap in the regions away from the focal spot. The peak (+ve) pressure predicted based on the nonlinear Westervelt equation is considerably higher in the focal region in comparison with that predicted through the linear model. The possible implications of these observations on the efficiency of the HIFU based non-invasive destruction of abnormal cells may be realized in the form of significant differences in the resultant temperature levels predicted by the two approaches. The strong nonlinear effects enhance the heat-affected zone. Thus, for the same specification of the source used, while relatively higher peak pressures can be realized with the nonlinear model, there is also a possibility of experiencing some other phenomena such as boiling, cavitation etc. [34].

Volumetric heat generation has been calculated through temporal averaging of nonlinear propagation of wave for one time period and was kept steady for the temperature calculation until the system (transducer) is kept on. Following this approach, the pressure distribution calculated in Domain 1 has been employed to calculate the variation of volumetric heat generation due to the attenuation of the propagating wave. Figure 4(c) illustrates the axial variations of the linear and nonlinear volumetric heat generation in the tissue domain along the center line of the transducer. It is to be noted that the attenuation coefficient for water is significantly low, hence, the magnitude of volumetric heat generation in the region i.e., for $0<z<2\text{\hspace{0.17em}cm }$is almost zero (Q = 0). It is worth mentioning here that in clinical applications, this is one of the primary reasons behind the usage of ultrasound gels, which have their acoustic impedance almost close to that of water, as this layer drastically eliminates the possibility of any reflection of the ultrasound waves that would otherwise occur in the case of air-tissue interface. Figure 4(c) shows the gradual build-up of the net volumetric heat generation term beyond the water layer and peaks in the focal region of the transducer (z = 7 cm). Further, as expected, the nonlinear form of the Westervelt equation predicts relatively higher amount of heat generation in the focal region, in comparison with that of the linear profile.

The volumetric heat generated due to the absorption of the ultrasound wave intensity in the medium leads to a rise in temperature of the physical domain, a phenomenon that forms the basis of the HIFU-based therapeutic techniques. With the volumetric heat distribution known, Equation (4)(4) $${\rho }_t{C}_t\frac{\partial T}{\partial t}=k{\nabla }^{2}T\hspace{0.17em}-\hspace{0.17em}{\omega }_b{C}_b(T\hspace{0.17em}-\hspace{0.17em}{T}_{art})+Q$$(4) has been solved to determine the temperature distribution within the body of the physical domain. Temperature variation has been predicted at the focus of the simulation geometry for a total sonication period of 0.2 s, (Figure 4(d)). As can be seen, the peak of temperature is realized at the focal spot and its value is ≈84 °C, which is practically sufficient for complete necrosis of the abnormal cells. As expected, the nonlinear theory predicts higher temperature peak at the focal spot compared to the linear model.

4.3. Effect of various layers on the focal intensity

Acoustic wave propagates into the target volume through various tissue layers with different acoustic and thermal properties. Thus, to simulate the more realistic model of the tissue sample, multilayer modeling (Domain 2, Figure 1(b.ii)) has been carried out. To investigate the layers and liver tissue effect, HIFU intensity maps for various combinations (water, water + tissue, water + layer + water and water + layer + tissue) have been shown in Figure 5. Since water is a lossless medium, the acoustic energy is maximum if solely water is considered. Introducing a layer consisting of sub-layers of skin (2 mm), fat (1 cm), muscle (1.5 cm) in water reduces the peak acoustic energy and also slightly shifts the profile in the direction of the wave propagation (away from the transducer) (Table 3). However, the addition of liver tissue into the water medium (water + tissue) shifts the profile in the direction opposite to that of the propagating wave i.e., towards the transducer with further reduction in the peak. Furthermore, considering Domain 2, which consists of water, layers of fat, and liver tissue, shift in the profile is in the direction of wave propagation, which is relatively higher in comparison with the shift that is realized for only water and tissue domain (Domain 1). This concludes that the overall effect of the layer is the shift of the intensity profile in the direction of wave propagation with relatively reduced peak.

Figure 5. Intensity plot along the axial distance for different domains. The layer consists of skin (0.2 cm) + fat (1 cm) + muscle (1.5 cm). In all the cases, 2 cm of water is placed between transducer and the medium.

Table 3. Quantitative analysis of various medium intensity.

Medium	zfocalplane (cm)	In (W/cm^2)
Water	6.9616	15148.7
W + Tissue	6.7488	5886.36
W + L+Water	6.992	8383.46
W + L+Tissue	6.8552	5387.41

To check the individual effect of fat and muscle layers on the acoustic intensity, these two layers, with different thicknesses, have been considered separately. If a particular layer is sandwiched between the lossless medium i.e., water and the homogenous liver tissue, any change in the intensity profile will primarily depend on the properties of the sandwiched layer. Figure 6 illustrates the individual effects in lossless and tissue medium. In all the cases, the gap between the transducer and tissue medium is filled with water. As seen in Figure 6(a), a single layer of fat or muscle with different thicknesses, when combined with water, results into an opposite shift (with respect to the focal point) of the intensity profiles. However, in both cases, the peak of the acoustic energy reduces with an increase in the layer thickness due to higher absorption. The observed trend may be explained on the basis of the differences in the speed of propagation of the acoustic waves in the layers of fat and muscle compared to water. The speed of propagation of acoustic waves in fat is lower compared to that in water. Thus, with an increase in the thickness of the fat layer, the point of maximum acoustic intensity shifts away from the transducer. However, relatively higher propagation speed in muscle layer leads to a shift of profile towards the transducer. Furthermore, the differences in the absorption coefficient of these two layers also lead to the corresponding changes in the peak intensity values.

Figure 6. Intensity plots along the axial distance to show the effect of single layer (fat/muscle) sandwiched between (a) two layers of water and (b) water and liver tissue. For better clarity, the layer configurations have schematically been shown in both the figures.

Figure 6(b) shows the profiles of intensities when the single layer of fat/muscle has been sandwiched between water and liver tissue (this combination is different from that presented in Figure 6(a) wherein the single layer was combined only with water). Compared to the profiles shown in Figure 6(a), a significant reduction in the magnitude of peak intensity for any given layer is to be seen. In addition, Figure 6(b) also shows that the profiles corresponding to the fat layer show the similar nature of shift (away from the transducer) as was seen in the case of the combination shown in Figure 6(a). However, the corresponding shift in the profiles of the acoustic waves in muscles is ≈0.99 times of that in the liver tissue.

4.4. Effect of source parameters: samples with embedded inhomogeneity

Discussions presented in the previous sections pertaining to the cases wherein source pressure and the transducer frequency were fixed at 0.5 MPa and 1 MHz respectively for Domains 1 and 2. In order to quantify the influence of variations in these parameters on the resultant temperature field, simulations have been conducted for different source pressures, p_o (0.2, 0.5 and 0.7 MPa) and frequencies, f (0.8, 1.0 and 1.2 MHz) for Domain 3 (Figure 1(b.iii))¹. With reference to Figure 1(b.iii), the embedded inhomogeneity numerically simulates the local presence of cancerous cells in the body. The acoustic and thermal properties of the embedded inhomogeneity have been considered to be as${\mathrm{ }c}_0\mathrm{\hspace{0.17em}}$= 1550 m/s,${\mathrm{}\rho }_0$= 1000 kg/m³,${\mathrm{ }\alpha }_0$= 8.1 Np/m,$\mathrm{ }\beta \mathrm{ }$= 4.5, k = 0.55 W/m°C and C_p = 3800 J/kg°C. For all the cases considered here, the heating duration has been kept constant at 0.2 s. Results of the study have been presented in the form of two-dimensional temperature contours in Figure 7.

Figure 7. Temperature contour plots in the focal plane (z-x plane) for varying fundamental frequencies (0.8, 1 and 1.2 MHz) and of varying source pressures (0.2, 0.5 and 0.7 MPa).

It is to be seen that for any given transducer frequency, an increase in the strength of the source pressure results in an increase in the heat deposited area due to higher acoustic energy loss. On the other hand, irrespective of the source pressure magnitude, as the transducer frequency is increased, the heat affected area narrows down and the acoustic energy gets more concentrated in the focal region. This observation can be attributed to the dependence of the attenuation coefficient on the source frequency, which increases with an increase in the frequency. As a result, more acoustic energy gets dissipated, resulting in higher temperature rise. This effect is quite prominent in the focal region due to the generation of secondary harmonics, as can be clearly seen from Figure 8, which shows the two-dimensional temperature contours (in the x-y plane) at the focus for different frequencies (source pressure has been kept constant at 0.5 MPa). The temperature contours presented in the figure show that the acoustic energy loss near the focus is more, but reduces with an increase in the frequency as one moves away from the focal region.

Figure 8. Temperature contour plots in the x-y plane for the focus point (z = 7 cm) for different fundamental frequencies (a) 0.8 MHz, (b) 1 MHz and (c) 1.2 MHz. For each case, the source pressure has been taken 0.5 MPa and heating duration is 0.2 s.

For more direct comparison, Figure 9 shows the temperature change with respect to the axial direction for the different transducer frequencies considered in the simulations. As the source frequency increases, the peak magnitude of temperature also increases due to the deposition of more number of acoustic wave cycles per second. Furthermore, as the diffusivity of acoustic energy is inversely proportional to frequency, the extent of the lateral spread of temperature profiles (on either side of the focal spot) reduces with increasing frequency, resulting into relatively narrower profiles at higher values of source frequency.

Figure 9. Temperature variation along the axial direction for different frequencies.

4.5. Generation of multiple lesions and energy modulation

Laboratory-based experiments, as well as clinical applications of HIFU technique, have shown that, compared to the overall size of the tumor affected region of the biological sample, the lesion size for a given position of single focused HIFU transducer is significantly small. (Here the lesion is defined as the spatial stretch of the biological tissue sample over which the local tissue temperature exceeds the threshold value required for complete destruction of the tumorous cells as the tissue sample gets subjected to HIFU waves.) In order to circumvent such limitations and to cover the entire spread of the embedded tumor, lesions must be placed systematically side-by-side. This arrangement can be achieved by moving the transducer mechanically (keeping the sample static) [35,36], using multi-element phased array transducer [37] or with the use of multiple transducers using multidirectional heating scheme [38]. The acoustic intensity of each lesion can be adjusted depending upon the size of the targeted tumor and its location. However, the main concern in each of these methods is to minimize the heating of healthy intervening tissues and the total treatment time. In this context, an overall planning and methodology for the generation of multiple lesions becomes important.

4.5.1. Planning methodology

For uniform and complete tumor ablation, it is required to generate similar size of each lesion to cover the entire tumor-affected region [39]. However, due to the thermal diffusion effects, the lesion size increases even after the completion of the heating duration, which makes the initial spot undertreated while the later one gets over treated [36]. To overcome this limitation, pre-controlled unequal distribution of acoustic energy is required to produce the modulated size of the lesion. There are two parameters which can be altered to generate the transient and modulated temperature in the tissue domain for achieving controlled lesion size:

1. Acoustic power,

2. Heating duration.

In both the cases, a modulation factor is required which can be adjusted according to the lesion position and gives controlled lesion size. The modulation factor, according to the usage and treatment requirement, can be defined as: (13) $$m=1-\frac{1}{N-(i-1)}$$(13)

Here, m is the modulation factor for the i^th lesion (i = 1,2,…N-1), N represents the total number of lesions (including the fundamental lesion) required to fill the entire tumor-affected region in any given direction. It is to be noted that i = 0 stands for the fundamental lesion for which the value of the modulated factor (m) would be equal to one since no modulation is required in the original acoustic intensity provided by the transducer. The advantage of using the above scheme is that the last lesion can always have half of the acoustic power or duration of heating used for the fundamental lesion (i = 0). It is to be mentioned that in the present study, it is assumed that the properties of the medium do not significantly change with temperature. Thus, the modulation factor can be directly multiplied by the volumetric power generation or the heating duration to provide a different level of temperature distribution for each lesion.

In addition to the total number of lesions, another factor that plays an important role in the planning methodology is the inter-spatial distance (χ) of different lesions. Since the domain size remains fixed, according to the application, the number of lesions or the gap between each lesion can be varied. However, to place more number of lesions in a confined space, the gap between each lesion should be optimized. The extent of this gap depends on the radius of the transducer and on the total number of lesions required to be generated. In place of the transducer radius, the boundary of the tumor can be considered so that the lesion remains inside the tumor domain. (14) $$\chi =\frac{a}{N}$$(14) where, a is the transducer radius and N represents the total number of lesions. It can be clearly seen from the above relation that the inter-spatial gap between each lesion is constant for all the lesions independent of its size. (To make it more generalized and dependent on each lesion position and size, χ can also be considered as a function of modulation factor (m) which is dependent on the i^th lesion.)

4.5.2. Modulation in volumetric heat generation and heating time

As mentioned earlier, the modulation factor for each lesion can be directly multiplied with the respective volumetric power distribution and can be considered as the modulated acoustic power generated by the moving transducer for a particular lesion. It is to be noted that the simulations for multiple lesions generation have been performed for the two-layered (water and tissue) system i.e., for Domain 1 (shown in Figure 1(b.i)). Images shown in Figure 10(a–d) explain the generation of total five lesions, (N = 5), in the tissue domain from fundamental (i = 0) to lesion five (i = 4). Q, T, and TD, respectively, represent the volumetric power, temperature, and thermal dosage. Images shown in the first column of Figure 10 show the volumetric power generated for a given position of the transducer as it is translated along the vertical direction. As the transducer is moved along the vertical direction for creating more number of lesions, the acoustic energy deposited in the focal region tends to decrease. Ultimately, for the last lesion, it reaches a value that is half of the energy of the fundamental lesion. The second column of the figure illustrates the transient temperature distribution in the overall domain due to the generation of each lesion. After the creation of one lesion, once the transducer is moved to create the other lesion, the previously generated lesion starts to cool down. For computational modeling, the initial conditions get updated each time after the generation of every lesion. The third column of Figure 10 explains the thermal dosage, which demonstrates the actual lesion size for each lesion generated. It can be clearly seen that, as the deposited energy decreases, the size of the lesion that is achieved for a given position of the transducer (along with the vertical direction) also reduces.

Figure 10. Generation of multiple lesions: (a) 1 lesion (b) 2 lesions (c) 3 lesions and (d) 5 lesions. First, second and the third column respectively represent the volumetric power, corresponding transient temperature and thermal dosage for each lesion. In the contours shown, L_i (i = 1, 2, ., 5) represents the number of the lesion created.

Similar to the modulation of the acoustic power, the heating duration can also be modified or controlled depending on the location of any given lesion. The procedure for the calculation of modulation factor and inter spatial gap between two lesions remains the same. However, the modulation factor m of the specific lesion gets multiplied with the heating duration (0.2 s) employed for the creation of the fundamental lesion (L₁), while the acoustic power remains same for all the lesions. It further implies that every step movement of the transducer provides the same absorbed power density but its total energy, which is the product of the absorbed power density and heating duration, gets changed. Figure 11 shows the comparison between modulation of energy and modulation of time. Figure 11(a) illustrates the temperature variation (at focus z = 7 cm) with respect to time. From the figure, it can be clearly seen that for m = 1 (i.e., no modulation), the temperature peaks of all the lesions reach to the same magnitude, whereas with modulation, the temperature peaks get reduced.

Figure 11. Temperature for m = 1, $m\times \mathit{\text{heating}}\mathit{\text{ time}}$ and $m\times \mathit{\text{volumetric}}\mathit{\text{ power}}$ for different lesions (a) with respect to time, (b) along the axial direction.

It should be noted here that though the modulation affects the volumetric power (or heating time), the total energy remains the same. However, the modulation of volumetric power gives slightly lower peak compared to the same modulation if it is incorporated in the heating time. This can be clearly seen in Figure 11(b), where the variations of temperature values with respect to axial direction have been plotted. The reduction in the peak values of temperature at the focal spot (z = 7 cm) can be attributed to the fact that modulation in heating time gives equal heating time to each lesion, which is higher than the corresponding heating time used for the case of modulated heating duration. Thus, before the maximum values of temperature get realized, the thermal energy gets diffused, resulting in a lower peak. Thus, higher lag in heating period causes the reduction in the peak values of temperature because of the dominance of thermal dissipation effects. Also, the total treatment time can be expected to be lesser for the case of modulated heating time.

4.5.3. Delay time period

The delay time period includes the transducer movement time as well as the cooling off period. It depends on the treatment application, and here, it is set to be equal to 0.1 s (half of the total heating time). Figure 12 introduces the effect of the delay period, in which (a) shows the modulation of volumetric power and (b) for the modulation of heating time. Since here the lesions generated are independent of each other, the delay period only shifts the peak position. Also, it can be clearly seen that modulated heating time period provides more cooling period compared to the modulated volumetric power. For complete necrosis of the targeted tumor, any extent of the gap between the lesions is not recommended. However, if the inter-spatial gap becomes very small, the lesions start linking with each other and may get overexposed. Under these conditions, sufficient cooling period is required before generating the next lesion [8]. Thus, the delay time period for such a situation is expected to play a significant role.

Figure 12. Delay time period in (a) modulation of volumetric power generation, and (b) modulation of heating time.

4.5.4. Effect of total number of lesions and inter spatial gap

Figure 13 illustrates the effect of increase in the number of lesions on temperature distribution in the tissue domain. Increase in number of lesions (from 5 to 20) results in the reduction in the inter-spatial distance between any given two lesions and beyond a certain limit, the lesions start merging with each other making it difficult to predict the exact number of lesions created in the region of interest. Merging of different lesions significantly affects the resultant temperature distribution field around each lesion. In addition to the temperature realized at the focal spot, this phenomenon of lesions merging with each other also results in an increase in the temperature of the surrounding region, which under ideal conditions should remain un-affected in order to prevent any possible thermal damage to the surrounding healthy cells. Thus, for ensuring selective destruction of abnormal cells spread over a finite region in the tissue domain, the clinical application of HIFU requires optimization of the inter-spatial distances between the relative positions of different lesions created by traversing the transducer.

Figure 13. Effect of total number of lesions on the resultant temperature distributions realised in the x-z (left column) and x-y planes (right column).

As mentioned earlier, if the inter spatial gap (χ) between two lesions is very small, the lesions start merging into each other, which leads to the preheating of the next lesion area even before its generation. This can be clearly seen from Figure 14(a) where temperature variations with respect to time for various total number of lesions have been reported for lesion 2 (L₂) and 5 (L₅). As the total number of lesions increase for any given domain, the temperature peak also increases for the lesions under consideration (lesion 2 and lesion 5 as mentioned in Figure 14(a)). This is expected due to the fact that as the inter-spatial gap between the lesions reduces, each lesion temperature gets affected by the dissipated energy of the previously generated lesion. Figure 14(b) and (c) show the enlarged views of preheating time of lesion 2 (L₂) and lesion 5 (L₅) for varying total number of lesions. Following the similar trend, preheating temperature is also gets higher as the total number of lesions increases.

Figure 14. Temperature at the focus point for (a) lesion 2 and lesion 5 for total number of lesions 5, 10 and 20, enlarged version (b) for lesion 2 (L₂) (c) for lesion 5 (L₅) shows the preheating time.

5. Conclusions

For therapeutic applications, it is important to develop an in-depth understanding of the tissue heating and the corresponding thermal response of the tissue. The present work is an attempt to delve upon some of these aspects with the use of multi-physics modeling scheme and reports the thermal analysis of biological tissue phantoms that are subjected to high intensity focused ultrasound waves. Full form of nonlinear Westervelt equation has been numerically solved in a three-dimensional computational domain for determining the pressure field. The pressure field has been coupled with Penne's bio-heat transfer equation to obtain the resultant whole field temperature distributions for various process parameters. Simulations were reported for two-layer as well as multi-layered tissue phantoms. Results were presented in the form of pressure field, temperature distributions and the corresponding thermal dosage in the targeted region of the tissue phantoms. Simulation results revealed significant differences between the acoustic pressure fields obtained through the linear and nonlinear form of Westervelt equations. Consequently, the magnitude of maximum temperature achieved at the focal plane of the transducer, as predicted using nonlinear Westervelt equation, was found to be significantly higher than that realized using the linear form of the equation. The observed trend was attributed to the generation of the higher harmonics associated with the nonlinear propagation of the acoustic waves in the tissue medium. Presence of the fat and the muscle layers in multi-layered tissue phantoms resulted into a clear shift in the focal spot of the pressure waves, which in turn changed the position of the heat affected zone with respect to the original focal plane of the transducer. In the context of multiple lesion generation, the importance of energy modulation for optimizing the inter-spatial distance between two consecutive lesions was highlighted by controlled movement of the transducer (while keeping the tissue phantom static).

Disclosure statement

No potential conflict of interest was reported by the authors.

Notes

1 Following the reports available in the literature (for instance, [3]), frequencies in the range 0.8 to 1.6 MHz are generally employed for liver cancer surgery and hence hold significance for clinical applications. This has been one of the primary motivations for selecting the range of frequencies employed in the reported numerical simulations.