Modeling of soft tissue thermal damage based on GPU acceleration

Abstract

Hyperthermia treatments require precise control of thermal energy to form the coagulation zones which sufficiently cover the tumor without affecting surrounding healthy tissues. This has led modeling of soft tissue thermal damage to become important in hyperthermia treatments to completely eradicate tumors without inducing tissue damage to surrounding healthy tissues. This paper presents a methodology based on GPU acceleration for modeling and analysis of bio-heat conduction and associated thermal-induced tissue damage for prediction of soft tissue damage in thermal ablation, which is a typical hyperthermia therapy. The proposed methodology combines the Arrhenius Burn integration with Pennes’ bio-heat transfer for prediction of temperature field and thermal damage in soft tissues. The problem domain is spatially discretized on 3-D linear tetrahedral meshes by the Galerkin finite element method and temporally discretized by the explicit forward finite difference method. To address the expensive computation load involved in the finite element method, GPU acceleration is implemented using the High-Level Shader Language and achieved via a sequential execution of compute shaders in the GPU rendering pipeline. Simulations on a cube-shape specimen and comparison analysis with standalone CPU execution were conducted, demonstrating the proposed GPU-accelerated finite element method can effectively predict the temperature distribution and associated thermal damage in real time. Results show that the peak temperature is achieved at the heat source point and the variation of temperature is mainly dominated in its direct neighbourhood. It is also found that by the continuous application of point-source heat energy, the tissue at the heat source point is quickly necrotized in a matter of seconds, while the entire neighbouring tissues are fully necrotized in several minutes. Further, the proposed GPU acceleration significantly improves the computational performance for soft tissue thermal damage prediction, leading to a maximum reduction of 55.3 times in computation time comparing to standalone CPU execution.

Keywords:

• Thermal damage
• Arrhenius Burn integration
• Pennes’ bio-heat equation
• finite element method
• Graphics Processing Unit
• thermal ablation

1. Introduction

Hyperthermia treatments use energy sources such as radiofrequency, laser, focused ultrasound and microwaves to generate thermal injury to tumors. These therapies require precise control of thermal energy to form the coagulation zones which sufficiently cover the tumor without affecting surrounding healthy tissues. However, with the current clinical practice, it is difficult to achieve the precise control of thermal energy [1]. Consequently, complications such as burns to neighboring organs or the bile duct are often occurred due to excessive thermal energy, and in reverse tumor survivals are also frequently occurred due to insufficient thermal energy delivery.

Modeling of soft tissue thermal damage is of crucial importance for thermal dosage control and monitoring in hyperthermia procedures and developing associated surgical simulations [2, 3]. Radiofrequency thermal ablation (RFA) is a minimally invasive surgical hyperthermia procedure used to eradicate tumors by generating heat in the tumor tissues while maintaining the functionality of surrounding healthy organs [4]. The heat generated in the RFA induces cellular necrosis to the soft tissue, leading to irreversible tissue damage to burn the tumor to death. Without inducing unintended thermal damage to the surrounding healthy tissues, it is important to predict the extent of thermal-induced tissue damage from temperature field of soft tissues for effective clinical outcomes.

The tissue thermal damage mechanism involved in the hyperthermia therapy can be characterized by the Arrhenius Burn integration model [5], in which the heat transfer in biological soft tissues is described by the bio-heat conduction process. Various techniques were studied to model the bio-heat conduction process in 3-D space. A simple discretization technique is the finite difference method (FDM) [6–8], approximating spatial derivatives using difference equations; however, this method can be applied to regular geometric meshes only. The finite element method (FEM) is a popular discretization technique for irregular geometric meshes describing objects of irregular shapes such as human organs and tumors. Various finite element models were developed, such as the ones for catheter-based thermal ablation with ultrasonic heating sources [5, 9], RFA [10], and microwave ablation [11]. In general, the existing studies on FEM modeling of thermal ablation are mainly focused on a static analysis using commercial FEM analysis packages, without consideration of real-time computational performance required in clinical practice [12, 13]. The literature has reported that the computation time for simulating a thermal ablation procedure of 15 min is around 3–6 h [14, 15].

Computational acceleration using Graphics Processing Unit (GPU), which is a highly parallel and powerful processor, provides a solution to facilitate the computational process of FEM. GPU is much superior to the Central Processing Unit (CPU) for parallel computation in terms of large datasets. However, the utilization of GPU to facilitate FEM modeling of thermal ablation is still very limited. Borsic et al. [16] studied the simulation of thermal ablation using GPU; however, this work is based on the simplifications of physical laws and provides simple results only. Chen et al. [17] also studied a GPU-acceleration simulation of thermal therapy via non-contact microwave imaging; however, this GPU-acceleration technique is combined with FDM for regular meshes rather than with FEM for irregular meshes.

This paper presents a GPU-accelerated methodology for modeling of soft tissue thermal damage in thermal ablation, where the soft tissue thermal damage is determined according to the Arrhenius damage model. The Pennes’ bio-heat equation, which is employed to determine the temperature field in soft tissues, is spatially discretized using FEM and temporally discretized using FDM. The implementation of the proposed methodology is achieved via the HLSL (High-Level Shader Language) and associated Direct3D 11. Simulation results demonstrate that the proposed GPU methodology can produce an efficient prediction of soft tissue thermal damage.

2. Methodology

2.1. Tissue thermal damage and bio-heat transfer

When sufficient heat is attained to soft tissues, thermal denaturation is occurred. This irreversible damage to soft tissues can be quantified by the Arrhenius integration [5], given by Equation (1)(1) $$\phi ={\int }_0^{\mathrm{t}}\mu \hspace{0.17em}{\text{exp}}^{-E_a/RT}dt$$(1) . (1) $$\phi ={\int }_0^{\mathrm{t}}\mu \hspace{0.17em}{\text{exp}}^{-E_a/RT}dt$$(1) where $\phi$ denotes the thermal damage and has no dimension, $t$ is the exposure time at a given absolute temperature $T$, $\mu$ denotes the material parameter, which is also called the damage rate factor, $E_a$ is the energy of activation, and $R$ is the universal gas constant.

The thermal damage of soft tissues will take place at normal body temperatures if we employ Equation (1)(1) $$\phi ={\int }_0^{\mathrm{t}}\mu \hspace{0.17em}{\text{exp}}^{-E_a/RT}dt$$(1) straightaway for any temperature; therefore, it is common to assume that soft tissue thermal damage is occurred when tissue temperature is above 43 °C [5].

Soft tissue temperatures in the heating process of thermal ablation can be determined based on the Pennes’ bio-heat equation, which expresses the propagation of heat energy in soft tissues owing to blood flow, metabolic heat generation, and regional heat sources [18]. The Pennes’ bio-heat equation is given by (2) $$m_{ti}{C}_{ti}\frac{\partial T_{ti}}{\partial t}=\nabla ·\left(k_{ti}\nabla T_{ti}\right)+W_{bl}{C}_{bl}\left(T_{art}-T_{ti}\right)+Q_m+Q_r$$(2) where $\nabla ·$ denotes the divergence operator,$\mathrm{ }\nabla$ the gradient operator, and $m_{ti}$ the tissue density, $C_{ti}$ the specific heat capacity, $T_{ti}$ the temperature of soft tissues, $t$ the time, $k_{ti}$ the thermal conductivity of soft tissues, $W_{bl}$ the perfusion rate, $C_{bl}$ the specific heat capacity of blood, $T_{art}$ the arterial blood temperature, $Q_m$ the metabolic heat generation rate, and $Q_r$ the heat rate of regional heat sources [6].

2.2. Numerical formulation

To numerically determine soft tissue temperatures for computation of thermal damage, Equation (2)(2) $$m_{ti}{C}_{ti}\frac{\partial T_{ti}}{\partial t}=\nabla ·\left(k_{ti}\nabla T_{ti}\right)+W_{bl}{C}_{bl}\left(T_{art}-T_{ti}\right)+Q_m+Q_r$$(2) is spatially discretized using FEM and temporally using FDM. Applying a trial function $\left[H\right]$ for approximation of the field values (e.g., the temperature $T=\left[H\right]\left\{T\right\}$), Equation (2)(2) $$m_{ti}{C}_{ti}\frac{\partial T_{ti}}{\partial t}=\nabla ·\left(k_{ti}\nabla T_{ti}\right)+W_{bl}{C}_{bl}\left(T_{art}-T_{ti}\right)+Q_m+Q_r$$(2) can be written into (3) $$\left[A\right]\left\{\frac{\partial T_{ti}}{\partial t}\right\}=\left[K\right]\left\{T_{ti}\right\}+\left\{Q\right\}$$(3) where (4) $$\left[A\right]={\int }_{\mathrm{\Omega }}^{}{m}_{ti}{C}_{ti}{\left[H\right]}^T\left[H\right]d\mathrm{\Omega }$$(4) (5) $$\left[K\right]=-{\int }_{\mathrm{\Omega }}^{}{k}_{ti}{\left[B\right]}^T\left[B\right]+W_{bl}{C}_{bl}{\left[H\right]}^T\left[H\right]d\mathrm{\Omega }$$(5) (6) $$\left\{Q\right\}={\int }_{\mathrm{\Gamma }}^{}{\left\{q\right\}}^T\left\{n\right\}{\left[H\right]}^{T}d\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}^{}{W}_{bl}{C}_{bl}{T}_{art}{\left[H\right]}^T+Q_m{\left[H\right]}^T+Q_r{\left[H\right]}^{T}d\mathrm{\Omega }$$(6) (7) $$\left[B\right]=\left[\begin{array}{ccc}\frac{\partial H_1}{\partial x}& \frac{\partial H_2}{\partial x}& \dots \\ \frac{\partial H_1}{\partial y}& \frac{\partial H_2}{\partial y}& \dots \\ \frac{\partial H_1}{\partial z}& \frac{\partial H_2}{\partial z}& \dots \end{array}\right];\mathrm{ }{\left\{q\right\}}^T=\left[-k_{ti,x}\frac{\partial T_{ti}}{\partial x},\mathrm{ }-k_{ti,y}\frac{\partial T_{ti}}{\partial y},\mathrm{ }-k_{ti,z}\frac{\partial T_{ti}}{\partial z}\right];\mathrm{ }\left\{n\right\}=\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]$$(7) where $\left[A\right]$ is the thermal mass matrix; $\left[K\right]$ the thermal stiffness matrix; $\left\{Q\right\}$ the effective thermal load vector; $\left[B\right]$ the temperature gradient matrix; $\left\{q\right\}$ the heat flux; $\left\{n\right\}$ the unit vector in the outward direction normal to the boundary surface.

The temporal derivative $\frac{\partial T_{ti}}{\partial t}$ in Equation (3)(3) $$\left[A\right]\left\{\frac{\partial T_{ti}}{\partial t}\right\}=\left[K\right]\left\{T_{ti}\right\}+\left\{Q\right\}$$(3) may be approximated by a first-order explicit forward time integration scheme [19–22], yielding (8) $$\frac{\partial T_{ti}}{\partial t}\cong \frac{T_{ti}\left(t+\Delta t\right)-T_{ti}\left(t\right)}{\Delta t}$$(8) where $\Delta t$ is the step size of time integration, $T_{ti}\left(t\right)$ and $T_{ti}\left(t+\Delta t\right)$ are the temperatures at time $t$ and $t+\Delta t$, respectively.

3. GPU implementation

The GPU implementation of the proposed methodology is developed in C++ and HLSL with Microsoft’s Direct3D 11 APIs, where numerical computations are parallelized and accelerated by utilizing compute shaders. The shaders are the programmable components of graphics rendering pipeline and easy to manipulate and modify, suitable for general-purpose GPU programming.

The shaders with the embedded HLSL are carried out at various stages in the GPU rendering process. The shaders are set and executed in a sequential order shown in Table 1. Figure 1 illustrates the GPU-accelerated computation framework of the proposed methodology.

Figure 1. GPU-accelerated computation framework: the mesh data is loaded along with simulation parameters at the host (CPU), and the numerical computation is performed on the designed device (GPU) via a sequential execution of five compute shaders in the GPU rendering pipeline.

Table 1. Execution of compute shaders.

Phase	Shader	Function
1	Shader#1	Clear values of each element in the global coefficient matrix.
2	Shader#2	Calculate the matrices and vectors of each element.
3	Shader#3	Sum element inputs to form the global coefficient matrix and load vector for numerical iterations.
4	Shader#4	Perform the Gauss-Seidel iteration on the linear system of equations.
5	Shader#5	Store the computed temperatures and thermal damage for the subsequent copy from the device (GPU) to host (CPU).

4. Results and Discussion

Simulation analyses are conducted to evaluate the performance of the proposed methodology. It first presents results of thermal damage, followed by performance evaluation on the GPU-accelerated implementation.

4.1. Thermal damage

Thermal damage is simulated using the developed methodology by applying a point-source temperature to a soft tissue. The test object, which is in cubic shape, consists of 20 × 20 × 20 cubes. The side length of each cube is 0.1 mm. Each cube is further divided into six linear tetrahedrons. The temperature at the point source is 54 °C, and it is exerted to the soft tissue in the duration of 300 s. The time increment $\Delta t$ is 0.002 s.

Figure 2 presents the temperature distribution at 30 s after the point-source temperature is applied. The temperature distribution reaches a steady state at which the variation of temperature is very small. It is noted that the variation of temperature is mainly dominated in the direct neighborhood of the point source, and the peak temperature (54 °C) is achieved at the point source. This is mainly because the heat flow input is set at the source node only while the heat flow inputs at other nodes are zero. As mentioned in Section 2.1, it is supposed that thermal damage is occurred at tissue temperatures greater than 43 °C; therefore, over the simulated time period, thermal damage will be occurred at the nodes in the neighborhood of the heat source if the temperatures are greater than the threshold value 43 °C.

Figure 2. Distribution of temperature ($T$) for the mid-plane ($xy$) at 30 s, where the temperature at the point source is 54 °C.

The computed thermal damage is quantitively assessed by comparing the data to a threshold value $\phi =1.0$ which denotes the beginning of the irreversible tissue damage [23]. The duration of applied heating (maintaining constant temperature) is defined as the time $t$ to cellular necrosis, i.e., reaching a thermal damage of $\phi =4.6$ [9].

Owing to the large drop of temperature in the direct neighborhood of the heat source node, the temperature difference between the heat source node and its immediately neighboring nodes is important. Considering Equation (1)(1) $$\phi ={\int }_0^{\mathrm{t}}\mu \hspace{0.17em}{\text{exp}}^{-E_a/RT}dt$$(1) , the rate of thermal damage quickly increases with the increase of temperature. As illustrated in Figure 3, the tissue at the heat source node necrotizes quickly within 20 s, reaching a thermal damage of $\phi =4.6$, while directly neighboring tissues commence necrotising since 60 s and have yet necrotized completely until 300 s. This reveals that while the lower and safer treatment temperature of 54 °C can quickly necrotize the tissue at the heat source node in a matter of seconds, it takes several minutes of continuous application for neighboring tissues to necrotize. Figure 4 presents the time needed to attain cellular necrosis ($\phi =4.6$) with reference to constant temperature applied during thermal treatment. It can be seen that the time to necrotize the tissue is drastically decreased with the increase of temperature. It requires more than 1000 s to necrotize the tissue at the operation temperature 44 °C, 200 s at 48 °C and 25 s at 52 °C. When the temperature is above 54 °C the tissue will be necrotized almost instantly. Figure 5 presents the final distribution of thermal damage for the mid-plane of the cubic tissue, where the tissue in the direct neighborhood of the point-source is completely necrotized. It also reflects that several minutes of continuous application of the point-source heat energy are required to fully necrotize the entire neighboring tissues.

Figure 3. Transient thermal damage $\phi$ at the mid-plane ($xy$) at (a) 1 s; (b) 10 s; (c) 20 s; (d) 60 s; (e) 150 s; and (f) 300 s; the tissue at the heat source node necrotizes quickly within 20 s, while directly neighboring tissues have not necrotized completely until 300 s.

Figure 4. Time to cellular necrosis (seconds) versus tissue temperature (°C); the time to necrotize the tissue is drastically decreased with the increase of temperature, and the tissue will be necrotized almost instantly when the temperature is above 54 °C.

Figure 5. Thermal damage $\phi$ distribution at the mid-plane ($xy$) at 1200 s.

4.2. GPU performance

Simulations were conducted at different mesh resolutions while maintaining the same simulation parameters to obtain computation times for iterations on the GPU. The GPU simulation was executed on an NVIDIA GTX 770M, whereas the CPU was Intel(R) Core(TM) i7-4700MQ. Figure 6 shows the average GPU solution time for a single iteration. It can be seen that the GPU solution time is linearly increased with the increase of the node number. The horizontal red line indicates the mark of 40 ms. The intersection point (indicated by the two dash lines) indicates around 22,000 nodes, that is, about 120,000 tetrahedrons, for modeling of bio-heat conduction and associated thermal damage can be solved in a short time to achieve real-time simulation [24–27]. After the mark of 22,000 nodes, greater node numbers will produce solution times greater than 40 ms. Therefore, the mesh resolution for simulation of bio-heat conduction and associated tissue thermal damage in real-time is approximately 22,000 nodes in the proposed GPU implementation. The GPU implementation provided a performance improvement of 53 times compared to the CPU counterpart at this node counts, and it achieved a maximum improvement of 55.3 times using the designated GPU and CPU configuration.

Figure 6. GPU solution time versus the number of nodes: a mesh of 22,000 nodes can be simulated at 40 ms (intersecting lines).

5. Conclusion

In this paper, a GPU-accelerated methodology is presented for efficient prediction of soft tissue thermal damage in RFA. The proposed methodology employs the Arrhenius Burn model and Pennes’ bio-heat equation for computation of soft tissue thermal damage. It applies FEM and FDM to spatially and temporally discretize the Pennes’ bio-heat equation, which is solved with the acceleration of GPU hardware. Simulation results demonstrate that the proposed methodology can efficiently predict soft tissue thermal damage. Future research work includes Compute Unified Device Architecture (CUDA) implementation of the proposed methodology and consideration of thermal-mechanical response of soft tissues in the heating process of RFA.

Disclosure statement

No potential conflict of interest was reported by the authors.