Non-Fourier based thermal-mechanical tissue damage prediction for thermal ablation

ABSTRACT

Prediction of tissue damage under thermal loads plays important role for thermal ablation planning. A new methodology is presented in this paper by combing non-Fourier bio-heat transfer, constitutive elastic mechanics as well as non-rigid motion of dynamics to predict and analyze thermal distribution, thermal-induced mechanical deformation and thermal-mechanical damage of soft tissues under thermal loads. Simulations and comparison analysis demonstrate that the proposed methodology based on the non-Fourier bio-heat transfer can account for the thermal-induced mechanical behaviors of soft tissues and predict tissue thermal damage more accurately than classical Fourier bio-heat transfer based model.

KEYWORDS:

• bio-heat transfer
• non-Fourier
• soft tissue
• thermal damage
• thermo-mechanical

Introduction

Thermal ablation is a minimally invasive thermal treatment for curing unrespectable tumors located inside human organs such as the liver and stomach. Currently, the thermal ablation process is controlled using temperature elevation to form tissue coagulation zones.¹ This method is inaccurate, as an optimized temperature field does not necessarily imply optimized tissue damage. It also ignores the effect of thermal-mechanical deformation during the thermal therapy. In fact, even a small variation in thermal-induced mechanical deformation can lead to various effects such as protein denaturation, altered production of hormones, and suppressed immune response.² Therefore, it is an absolute necessity to study tissue response under thermal and mechanical loads to provide an effective guidance for precise control of thermal ablation tasks.

Several authors have studied thermal-mechanical behaviors of soft tissues under thermal therapy. Xu et al.³ treated human skin as a kind of composite material, laminate, where heat-induced mechanical response (thermal stress) in different layers of skin was calculated using traditional composite material mechanics. However, this method is only limited to layered soft tissues and how thermal-mechanical response is incorporated into existed damage model to qualify thermal-mechanical damage was not discussed. In order to identify the influence of thermally induced mechanical deformation on thermal damage prediction, Shen et al.^4,5 developed a tissue damage model using Fourier bio-heat transfer equation. It shows that thermally induced mechanical deformation decreases the activation energy for protein denaturation, making soft tissues more easily to be damaged. Recently, the author's group developed a model to predict thermally induced mechanical deformation and thermal damage of soft tissues by combining the Fourier bio-heat transfer equation with the theory of linear thermo-elasticity.^6,7 In general, most of the existing methods are mainly dominated by the classical Fourier bio-heat transfer equation,⁸ which assumes infinite speed of heat propagation in soft tissues. In fact, a change temperature in a local area does not instantaneously cause a perturbation on the temperature distribution in the rest area of the medium. The heat transfer in a tissue like non-homogenous medium needs a characteristic time to propagate to the nearest point. This non-Fourier thermal behavior has been observed in biological tissues via physical experiments.^9-11

This paper gives a new methodology for modeling of thermal-mechanical behaviors and associated damage of soft tissues during thermal ablation, where the modeling process combines non-Fourier bio-heat transfer, continuum mechanics, as well as non-rigid motion of dynamics to predict and analyze temperature distribution and thermal-induced mechanical deformations of soft tissues. An improved tissue damage model is developed by accounting for thermal-induced mechanical damage. The finite difference scheme and numerical algorithms are utilized to construct the proposed thermal-mechanical tissue damage model. Simulations and comparison analysis with the classical Fourier bio-heat transfer have been performed to evaluate the performance of the proposed methodology.

Material model and numerical modeling

Non-Fourier bio-heat transfer equation

The classical Fourier based bio-heat transfer equation considering coupled thermo-elastic effect can be described as,¹²(2) $\rho C\frac{\partial T}{\partial t}+\beta T_0\frac{\partial \left(\nabla u\right)}{\partial t}=k{\nabla }^{2}T+C_b{\omega }_b{\rho }_b\left(T_b-T\right)+Q_m+Q_{ext}$(2) where ρ is the density soft tissues, C the specific heat, β the coefficient of thermal expansion, u the displacement vector, k the thermal conductivity, T the temperature at time t, ${\displaystyle {\omega }_b}$ the blood perfusion, ${\displaystyle {\rho }_b}$ the density of blood, C_b the specific heat of blood, T_b the temperature of blood vessel, T₀ the initial temperature, ${\displaystyle Q_{ext}}$ the external heat source, and Q_m the generated metabolic heat.

Applying the concept of finite heat propagation velocity,³ a modified unsteady heat transfer equation can be obtained(2) $q\left(r,t+{\tau }_q\right)=-k\nabla \text{T}\left(r,t\right)$(2) where ${\displaystyle {\tau }_q}$ is the characteristic time of biological soft tissues.

Using first-order Taylor expansion, the above equation becomes(3) $q\left(r,t+{\tau }_q\right)=q\left(r,t\right)+{\tau }_q\frac{\partial q\left(r,t\right)}{\partial t}$(3)

Substituting (3(3) $q\left(r,t+{\tau }_q\right)=q\left(r,t\right)+{\tau }_q\frac{\partial q\left(r,t\right)}{\partial t}$(3) ) into (1(2) $\rho C\frac{\partial T}{\partial t}+\beta T_0\frac{\partial \left(\nabla u\right)}{\partial t}=k{\nabla }^{2}T+C_b{\omega }_b{\rho }_b\left(T_b-T\right)+Q_m+Q_{ext}$(2) ), the non-Fourier bio-heat transfer equation can be obtained,(4) ${\tau }_p\rho C\frac{{\partial }^{2}T}{\partial t^2}+\left({\tau }_p{C}_b{\omega }_b{\rho }_b+\rho C\right)\frac{\partial T}{\partial t}+{\tau }_p\beta T_0\frac{{\partial }^2\left(\nabla u\right)}{\partial t^2}+\beta T_0\frac{\partial \left(\nabla u\right)}{\partial t}=k{\nabla }^{2}T+{\tau }_p\frac{\partial Q_m}{\partial t}+{\tau }_p\frac{\partial Q_{ext}}{\partial t}+C_b{\omega }_b{\rho }_b\left(T_b-T\right)+Q_m+Q_{ext}$(4)

Linear elastodynamics

Biologically, soft tissues are complex in terms of material compositions and structural formations. While soft tissue structure shows time-dependent, nonlinear and anisotropic behaviors, in terms of small deformation caused by thermal load, soft tissues can be investigated by linear thermo-elastic models to a high precision.

From the constitutive elastic material law under thermal loads, the stress tensor is related to the strain tensor and temperature change θ, which equals ${\displaystyle T-T_0}$(5) ${\sigma }_{ij}=\mu \left(u_{i,j}+u_{j,i}\right)+\lambda u_{k,k}{\delta }_{ij}-\beta \theta {\delta }_{ij}$(5) where ${\displaystyle u_{ij}}$ is the displacement components, and ${\displaystyle {\delta }_{ij}}$ is the Kronecker's symbol defined as${\delta }_{ij}=\{\begin{array}{l}1 for i=j ,\\ 0 for i\ne j ,\end{array}$

Rewriting the governing equation of non-rigid mechanics of motion as(6) $\rho {\ddot{u}}_{\hspace{0.17em}i}={\sigma }_{ij,j}+F_i$(6)

Rearranging the motion equation as,(7) $\rho {\ddot{u}}_{\hspace{0.17em}i}={\left(\mu \left(u_{i,j}+u_{j,i}\right)+\lambda u_{k,k}{\delta }_{ij}-\beta \theta {\delta }_{ij}\right)}_{,j}+F_i$(7) where λ and μ are the Lame constants, F is the exerted external force and,$2\mu =\frac{E}{1+\nu },\lambda =\frac{\nu E}{\left(1+\nu \right)\left(1-2\nu \right)},\beta =\frac{\alpha E}{1-2\nu }=\alpha \left(3\lambda +2\mu \right),\mu =G$

Constructing the proposed model on a regular cubic mesh straightforwardly using a finite difference scheme, gives the following global characteristic equation for all mesh points:(8) $M\underset{\hspace{0.17em}\_}{\ddot{u}}+C\underset{\hspace{0.17em}\_}{\dot{u}}+K\underset{\_}{u}=F, \underset{\_}{u}=\left\{\begin{array}{l}U\\ \theta \end{array}\right\}$(8)

To solve Eq. (8(8) $M\underset{\hspace{0.17em}\_}{\ddot{u}}+C\underset{\hspace{0.17em}\_}{\dot{u}}+K\underset{\_}{u}=F, \underset{\_}{u}=\left\{\begin{array}{l}U\\ \theta \end{array}\right\}$(8) ), it is necessary to specify the boundary conditions under concerning. Here, the boundary condition for solving Eq. (8(8) $M\underset{\hspace{0.17em}\_}{\ddot{u}}+C\underset{\hspace{0.17em}\_}{\dot{u}}+K\underset{\_}{u}=F, \underset{\_}{u}=\left\{\begin{array}{l}U\\ \theta \end{array}\right\}$(8) ) is realized by Dirichlet boundary condition, in which specified temperature and displacement are enforced to the related points on the solution domain boundaries at all times.

Thermal-mechanical damage prediction model

To predict thermal damage with consideration of thermal induced mechanical deformation effect, the relevant governing equation is expressed in Eq. (9(9) $\text{Ω}=\varsigma \underset{0}{\overset{t}{{\displaystyle \int }}}\text{exp}\left(-\frac{\text{Δ}E-E_{mechanical}}{R{T}_{\left(x,t\right)}}\right)dt$(9) ). Basically, it assumes that soft tissues subject to thermal stress cause tissue deformation. Thus, tissues are more easily being damaged under such a condition, because tissue protein denaturation may occur at relatively small activation energy.(9) $\text{Ω}=\varsigma \underset{0}{\overset{t}{{\displaystyle \int }}}\text{exp}\left(-\frac{\text{Δ}E-E_{mechanical}}{R{T}_{\left(x,t\right)}}\right)dt$(9) where ${\displaystyle E_{mechanical}}$ is the strain energy per mole, which is represented by Eq. (10). ${\displaystyle V_{mol}}$ stands for molar volume of target tissue.(10) $E_{mechanical}=\frac{1}{2}{V}_{mol}{\displaystyle \sum }_{i,j=1}^3{\sigma }_{ij}{ϵ}_{ij}$(10)

Results and discussion

Several simulations were conducted to test the performance of proposed thermal-mechanical model based on non-Fourier bio-heat transfer. Consider a cubic volumetric tissue model, which contains 1000 elements with 1331 nodes, with a needle inserted inside (see Fig. 1). The needle has 5 evenly spaced point sources at the front to generate heat energy and they are positioned in the Y direction at the center of the middle horizontal plane inside the cubic tissue. The needle generates thermal energy at the 5 points at ${\displaystyle 4\text{×}{10}^8\text{ W}}$ to heat up the tissue model. Table 1 shows the thermal and mechanical parameters used for the cubic tissue under thermal loads.^13-18 The boundary of the cubic model is fixed and its temperature is set to 310 K. The cubic model has the initial temperature of 310 K at the rest state.

Figure 1. Cubic tissue model with a needle inserted inside.

Table 1. Tissue parameters and constants.

Parameters	Value
Density $\rho$; $kg/m^3$	1060
Specific heat$C$; $J/kgK$	4192
Thermal conductivity $\text{ k}$; $\text{W}/\text{mK}$	0.613
Thermal expansion coefficient $\alpha$; $1/\text{°C}$	$1\times {10}^{-4}$
Young's modulus $E$; $Pa$	$1.1\times {10}^6$
Poisson's ratio $\nu$;	0.48
Blood perfusion rate ${\omega }_b$; $kg/\left(m^{3}s\right)$	0.5
Arterial temperature $T_b$; K	310
Metabolic heat generation rate $Q_m$; $W/m^3$	33800
Blood specific heat $C_b$; $J/kgK$	3600
Characteristic Time ${\tau }_q$; s	16
Activation energy $\Delta E$; $J/mol$	667000
Frequency factor $\varsigma$; $1/s$	$1.98\times {10}^{106}$
Gas constant $R$; $J/\left(molK\right)$	8.314

Temperature distributions

To study the thermal behavior of the cubic tissue, 3 observation points P₁ (0.005, 0.005, 0.005), P₂ (0.005, 0.005, 0.006) and P₃ (0.005, 0.005, 0.007) (blue points in Fig. 1), were chosen to observe the temperature variation inside the cubic tumor during the heating process. Figure 2 illustrates the temperature distributions at these 3 points within a heating period of 120 s for both models. As we can see, at the initial stage of heating (within about 40 s), the non-Fourier model involves oscillations in the temperature rise. After the initial stage of heating, with the increase of the heating time, the temperature rise for the non-Fourier model is larger than Fourier model, but the difference between both models gradually becomes smaller and finally the temperature curves of both models are converged to each other. This demonstrates that although the Fourier and non-Fourier models have similar accuracy for temperature prediction after a sufficient time of heating, the non-Fourier bio-heat model accounts for the transition state of bio-heat transfer, thus more suitable for the prediction of tissue's instantaneous thermal-mechanical behavior.

To further study the thermal behavior of the cubic tumor, the plane at Y = 0.005 m, which is perpendicular to the heating line, was selected to observe the temperature variation during the heating process. Figure 3 shows the temperature distributions on this plane at the initial time of 10 s for both Fourier and non-Fourier models. It can be seen that the temperature distributions of both models have a similar behavior with the highest temperature at the center of the selected plane, which is also the center of the heating sources. However, the non-Fourier model has smaller temperature rise than the Fourier model. This is because the non-Fourier model involves a process of time relaxation, leading to the smaller temperature rise at the initial heating stage.

Figure 2. Temperature distributions at certain points for both Fourier and non-Fourier models.

Figure 3. Comparison of temperature distributions between the Fourier (left) and non-Fourier (right) based methods.

Thermal-induced mechanical deformations

Trials were also conducted to study the thermal-induced deformation behaviors of the cubic tumor. Figure 4 shows the strain distributions at the above plane (Y = 0.005 m) at the heating time of 10 s for both Fourier and non-Fourier based methods. It can be seen that following the temperature distributions shown in Fig. 3, the strain in the Y direction of the non-Fourier based method is also smaller than that of the Fourier based method. It should be noted that as the selected observation plane is perpendicular to the Y axis, the strains in the X and Z directions are of a similar behavior, but are different from that of the strain in the Y direction. The above similar phenomenon can also be observed for the stresses in the 3 axis directions shown in Fig. 5. Figure 6 shows the strain energy distributions for both models. It can be seen that the non-Fourier model has smaller strain energy than the Fourier model at the initial heating time of 10 s. This proves that the strain energy distribution also follows the temperature distributions shown in Fig. 3. From the above, it is evident that the temperature distribution has a significant effect on thermal-induced strain and stress as well as strain energy.

Figure 4. Comparison of strain distributions between the Fourier (top) and non-Fourier (bottom) based methods.

Figure 5. Comparison of stress distributions between the Fourier (top) and non-Fourier (bottom) based methods.

Figure 6. Comparison of strain energy distribution between the Fourier (left) and non-Fourier (right) based methods.

Damage distribution

Simulations were also conducted to test the performance of the proposed thermal-mechanical tissue damage model. Figure 7 shows the damage comparison between the Fourier and non-Fourier models at the 3 selected points P₁-P₃. As the non-Fourier based method involves the transition state at the initial heating stage, at the points (e.g. P1 and P2) nearer to the center of the heat sources (i.e. the center of the selected observation plane), it needs more time to reach the threshold of irreversible tissue thermal damage (Ω = 1)³ than the Fourier based method. This phenomenon basically follows the above analysis at initial heating time of 10 s for points P1 and P2, where both temperature and strain energy prediction values were larger for the Fourier model. However, at the points (e.g., P3) further away from the center of the heat sources, it needs less time to reach the tissue damage threshold than the Fourier based method. This can be attributed to larger temperature prediction values toward the time point where the damage threshold was reached (around 60 s) for the non-Fourier model at point P3, as shown in Fig. 2. Moreover, Fig. 8 shows the damage distributions at the above plane (Y = 0.005 m) at the end heating time of 120 s. It can be seen that the final tissue damage zone for the non-Fourier based method (in green) is clearly larger for that of the Fourier based method (in Bangladesh green).

Figure 7. Damage vs. time comparison.

Figure 8. Comparison of tissue damage.

Conclusion

This paper presents a new methodology to predict thermal-induced mechanical deformation and associated tissue damage of soft tissues under thermal load. This methodology establishes a thermal-mechanical model and tissue damage model based on non-Fourier bio-heat transfer and continuum mechanics. Simulations and comparison analysis demonstrate that the proposed methodology based on non-Fourier bio-heat transfer accounts for thermal-induced mechanical deformation of soft tissues more accurately, leading to more accurate prediction of tissue thermal damage. Future work will focus on the construction of the proposed methodology on irregular human organ models using finite element method.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.