Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/taylor-francis/modeling-of-soft-tissue-thermal-damage-based-on-gpu-acceleration-LTLUQWuB3K
References
References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.
- Publisher
- Taylor & Francis
- Copyright
- © 2019 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
- ISSN
- 2469-9322
- DOI
- 10.1080/24699322.2018.1557891
- Publisher site
- See Article on Publisher Site
Abstract
COMPUTER ASSISTED SURGERY 2019, VOL. 24, NO. S1, 5–12 https://doi.org/10.1080/24699322.2018.1557891 RESEARCH ARTICLE a a a b c d Jinao Zhang , Jeremy Hills , Yongmin Zhong , Bijan Shirinzadeh , Julian Smith and Chengfan Gu a b School of Engineering, RMIT University, Bundoora, Australia; Robotics and Mechatronics Research Laboratory, Department of Mechanical and Aerospace Engineering, Monash University, Clayton, Australia; Department of Surgery, School of Clinical Sciences at Monash Health, Monash University, Clayton, Australia; City of Whittlesea, Mill Park, Australia KEYWORDS ABSTRACT Thermal damage; Arrhenius Hyperthermia treatments require precise control of thermal energy to form the coagulation zones Burn integration; Pennes’ which sufficiently cover the tumor without affecting surrounding healthy tissues. This has led bio-heat equation; finite modeling of soft tissue thermal damage to become important in hyperthermia treatments to element method; Graphics completely eradicate tumors without inducing tissue damage to surrounding healthy tissues. This Processing Unit; paper presents a methodology based on GPU acceleration for modeling and analysis of bio-heat thermal ablation 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 expli- cit 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 effect- ively 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 tem- perature 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 per- formance for soft tissue thermal damage prediction, leading to a maximum reduction of 55.3 times in computation time comparing to standalone CPU execution. 1. Introduction in hyperthermia procedures and developing associated surgical simulations [2, 3]. Radiofrequency thermal Hyperthermia treatments use energy sources such as ablation (RFA) is a minimally invasive surgical hyper- radiofrequency, laser, focused ultrasound and micro- thermia procedure used to eradicate tumors by gener- waves to generate thermal injury to tumors. These ating heat in the tumor tissues while maintaining the therapies require precise control of thermal energy to functionality of surrounding healthy organs [4]. The form the coagulation zones which sufficiently cover heat generated in the RFA induces cellular necrosis to the tumor without affecting surrounding healthy tis- the soft tissue, leading to irreversible tissue damage to sues. However, with the current clinical practice, it is burn the tumor to death. Without inducing unin- difficult to achieve the precise control of thermal tended thermal damage to the surrounding healthy energy [1]. Consequently, complications such as burns tissues, it is important to predict the extent of ther- to neighboring organs or the bile duct are often mal-induced tissue damage from temperature field of occurred due to excessive thermal energy, and in soft tissues for effective clinical outcomes. reverse tumor survivals are also frequently occurred The tissue thermal damage mechanism involved in due to insufficient thermal energy delivery. the hyperthermia therapy can be characterized by the Modeling of soft tissue thermal damage is of crucial Arrhenius Burn integration model [5], in which the heat importance for thermal dosage control and monitoring transfer in biological soft tissues is described by the CONTACT Yongmin Zhong yongmin.zhong@rmit.edu.au School of Engineering, RMIT University, Bundoora, VIC 3082, Australia 2019 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 6 J. ZHANG ET AL. bio-heat conduction process. Various techniques were computational process of FEM. GPU is much superior to studied to model the bio-heat conduction process in the Central Processing Unit (CPU) for parallel computa- 3-D space. A simple discretization technique is the finite tion in terms of large datasets. However, the utilization difference method (FDM) [6–8], approximating spatial of GPU to facilitate FEM modeling of thermal ablation is derivatives using difference equations; however, this still very limited. Borsic et al. [16] studied the simulation method can be applied to regular geometric meshes of thermal ablation using GPU; however, this work is only. The finite element method (FEM) is a popular dis- based on the simplifications of physical laws and pro- cretization technique for irregular geometric meshes vides simple results only. Chen et al. [17] also studied a describing objects of irregular shapes such as human GPU-acceleration simulation of thermal therapy via non- organs and tumors. Various finite element models were contact microwave imaging; however, this GPU-acceler- developed, such as the ones for catheter-based thermal ation technique is combined with FDM for regular ablation with ultrasonic heating sources [5, 9], RFA [10], meshes rather than with FEM for irregular meshes. and microwave ablation [11]. In general, the existing This paper presents a GPU-accelerated methodology studies on FEM modeling of thermal ablation are mainly for modeling of soft tissue thermal damage in thermal focused on a static analysis using commercial FEM ana- ablation, where the soft tissue thermal damage is deter- lysis packages, without consideration of real-time com- mined according to the Arrhenius damage model. The putational performance required in clinical practice [12, Pennes’ bio-heat equation, which is employed to deter- 13]. The literature has reported that the computation mine the temperature field in soft tissues, is spatially dis- time for simulating a thermal ablation procedure of cretized using FEM and temporally discretized using 15 min is around 3–6h [14, 15]. FDM. The implementation of the proposed methodology Computational acceleration using Graphics Processing Unit (GPU), which is a highly parallel and is achieved via the HLSL (High-Level Shader Language) powerful processor, provides a solution to facilitate the and associated Direct3D 11. Simulation results 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). 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. COMPUTER ASSISTED SURGERY 7 @T demonstrate that the proposed GPU methodology can ti ðÞ ðÞ m C ¼r k rT þ W C T T þ Q þ Q ti ti ti ti bl bl art ti m r @t produce an efficient prediction of soft tissue ther- (2) mal damage. where r denotes the divergence operator, r the gra- dient operator, and m the tissue density, C the spe- ti ti 2. Methodology cific heat capacity, T the temperature of soft tissues, t ti 2.1. Tissue thermal damage and bio-heat transfer the time, k the thermal conductivity of soft tissues, ti W the perfusion rate, C the specific heat capacity of bl bl When sufficient heat is attained to soft tissues, thermal blood, T the arterial blood temperature, Q the art m denaturation is occurred. This irreversible damage to metabolic heat generation rate, and Q the heat rate soft tissues can be quantified by the Arrhenius inte- of regional heat sources [6]. gration [5], given by Equation (1). E =RT u ¼ l exp dt (1) 2.2. Numerical formulation To numerically determine soft tissue temperatures for where u denotes the thermal damage and has no computation of thermal damage, Equation (2) is spa- dimension, t is the exposure time at a given absolute tially discretized using FEM and temporally using FDM. temperature T, l denotes the material parameter, which Applying a trial function ½H for approximation of the is also called the damage rate factor, E is the energy of field values (e.g., the temperature T ¼ H fTg), Equation activation, and R is the universal gas constant. (2) can be written into The thermal damage of soft tissues will take place at normal body temperatures if we employ @T ti ½ ½ fg A ¼ Kfg T þ Q (3) ti Equation (1) straightaway for any temperature; there- @t fore, it is common to assume that soft tissue ther- where mal damage is occurred when tissue temperature is above 43 C[5]. T ½ ½ ½ A ¼ m C H H dX (4) ti ti Soft tissue temperatures in the heating process of thermal ablation can be determined based on the T T ½ ½ ½ ½ ½ K ¼ k B B þ W C H H dX (5) ti bl bl Pennes’ bio-heat equation, which expresses the propa- gation of heat energy in soft tissues owing to blood ð ð T T T flow, metabolic heat generation, and regional heat fg ½ ½ ½ Q ¼ q fg n H dC þ W C T H þ Q H fg bl bl art m C X sources [18]. The Pennes’ bio-heat equation is given by þ Q H dX (6) Figure 2. Distribution of temperature (T) for the mid-plane (xy) at 30 s, where the temperature at the point source is 54 C. 8 J. ZHANG ET AL. 2 3 shape, consists of 20 20 20 cubes. The side length @H @H 1 2 ::: 6 7 of each cube is 0.1 mm. Each cube is further divided @x @x 6 7 6 7 @H @H into six linear tetrahedrons. The temperature at the 6 1 2 7 ½ B ¼ ::: ; 6 7 6 @y @y 7 point source is 54 C, and it is exerted to the soft tis- 6 7 4 5 @H @H sue in the duration of 300 s. The time increment Dt 1 2 ::: @z @z is 0.002 s. 2 3 Figure 2 presents the temperature distribution at @T @T @T ti ti ti 6 7 30 s after the point-source temperature is applied. The q ¼k ; k ; k ; n ¼ n fg 4 5 fg ti;x ti;y ti;z y @x @y @z temperature distribution reaches a steady state at which the variation of temperature is very small. It is (7) noted that the variation of temperature is mainly ½ ½ where A is the thermal mass matrix; K the thermal stiff- dominated in the direct neighborhood of the point ness matrix; fQg the effective thermal load vector; B the source, and the peak temperature (54 C) is achieved temperature gradient matrix; fqg the heat flux; fg n the at the point source. This is mainly because the heat unit vector in the outward direction normal to the bound- flow input is set at the source node only while the ary surface. heat flow inputs at other nodes are zero. As men- @T ti The temporal derivative in Equation (3) may be tioned in Section 2.1, it is supposed that thermal dam- @t approximated by a first-order explicit forward time age is occurred at tissue temperatures greater than integration scheme [19–22], yielding 43 C; therefore, over the simulated time period, ther- mal damage will be occurred at the nodes in the @T TðÞ t þ Dt T ðÞ t ti ti ti ffi (8) neighborhood of the heat source if the temperatures @t Dt are greater than the threshold value 43 C. where Dt is the step size of time integration, T ðtÞ and ti The computed thermal damage is quantitively ðÞ T t þ Dt are the temperatures at time t and t þ Dt, ti assessed by comparing the data to a threshold value respectively. u ¼ 1:0 which denotes the beginning of the irrevers- ible tissue damage [23]. The duration of applied heat- 3. GPU implementation ing (maintaining constant temperature) is defined as the time t to cellular necrosis, i.e., reaching a thermal The GPU implementation of the proposed methodology damage of u ¼ 4:6[9]. is developed in Cþþ and HLSL with Microsoft’s Direct3D Owing to the large drop of temperature in the dir- 11 APIs, where numerical computations are parallelized ect neighborhood of the heat source node, the tem- and accelerated by utilizing compute shaders. The shad- perature difference between the heat source node ers are the programmable components of graphics ren- and its immediately neighboring nodes is important. dering pipeline and easy to manipulate and modify, Considering Equation (1), the rate of thermal damage suitable for general-purpose GPU programming. quickly increases with the increase of temperature. As The shaders with the embedded HLSL are carried out illustrated in Figure 3, the tissue at the heat source at various stages in the GPU rendering process. The shad- node necrotizes quickly within 20 s, reaching a thermal ers are set and executed in a sequential order shown in damage of u ¼ 4:6, while directly neighboring tissues Table 1. Figure 1 illustrates the GPU-accelerated compu- commence necrotising since 60 s and have yet necro- tation framework of the proposed methodology. tized completely until 300 s. This reveals that while the lower and safer treatment temperature of 54 C can 4. Results and Discussion quickly necrotize the tissue at the heat source node in a matter of seconds, it takes several minutes of con- Simulation analyses are conducted to evaluate the perform- tinuous application for neighboring tissues to necro- ance of the proposed methodology. It first presents results tize. Figure 4 presents the time needed to attain of thermal damage, followed by performance evaluation cellular necrosis (u ¼ 4:6) with reference to constant on the GPU-accelerated implementation. temperature applied during thermal treatment. It can be seen that the time to necrotize the tissue is drastic- 4.1. Thermal damage ally decreased with the increase of temperature. It Thermal damage is simulated using the developed requires more than 1000 s to necrotize the tissue at methodology by applying a point-source temperature the operation temperature 44 C, 200 s at 48 C and to a soft tissue. The test object, which is in cubic 25 s at 52 C. When the temperature is above 54 C COMPUTER ASSISTED SURGERY 9 the tissue will be necrotized almost instantly. Figure 5 GPU. The GPU simulation was executed on an NVIDIA presents the final distribution of thermal damage for GTX 770M, whereas the CPU was Intel(R) Core(TM) i7- the mid-plane of the cubic tissue, where the tissue in 4700MQ. Figure 6 shows the average GPU solution the direct neighborhood of the point-source is com- time for a single iteration. It can be seen that the GPU solution time is linearly increased with the increase of pletely necrotized. It also reflects that several minutes the node number. The horizontal red line indicates of continuous application of the point-source heat the mark of 40 ms. The intersection point (indicated energy are required to fully necrotize the entire neigh- boring tissues. 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 4.2. GPU performance can be solved in a short time to achieve real-time Simulations were conducted at different mesh resolu- simulation [24–27]. After the mark of 22,000 nodes, tions while maintaining the same simulation parame- greater node numbers will produce solution times ters to obtain computation times for iterations on the greater than 40 ms. Therefore, the mesh resolution for Figure 3. Transient thermal damage u 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. 10 J. ZHANG ET AL. 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 u distribution at the mid-plane (xy)at 1200s. Figure 6. GPU solution time versus the number of nodes: a mesh of 22,000 nodes can be simulated at 40 ms (intersecting lines). COMPUTER ASSISTED SURGERY 11 [8] Li X, Zhong Y, Jazar R, et al. Thermal-mechanical simulation of bio-heat conduction and associated tis- deformation modelling of soft tissues for thermal sue thermal damage in real-time is approximately ablation. Biomed Mater Eng. 2014;24:2299–2310. 22,000 nodes in the proposed GPU implementation. [9] Prakash P, Diederich CJ. Considerations for theoretical The GPU implementation provided a performance modelling of thermal ablation with catheter-based improvement of 53 times compared to the CPU coun- ultrasonic sources: implications for treatment plan- ning, monitoring and control. Int J Hyperthermia. terpart at this node counts, and it achieved a max- 2012;28:69–86. imum improvement of 55.3 times using the [10] Wang Z, Aarya I, Gueorguieva M, et al. Image-based designated GPU and CPU configuration. 3D modeling and validation of radiofrequency inter- stitial tumor ablation using a tissue-mimicking breast phantom. Int J Cars. 2012;7:941–948. 5. Conclusion [11] Rattanadecho P, Keangin P. Numerical study of heat transfer and blood flow in two-layered porous liver In this paper, a GPU-accelerated methodology is pre- tissue during microwave ablation process using single sented for efficient prediction of soft tissue thermal and double slot antenna. Int J Heat Mass Transf. damage in RFA. The proposed methodology employs 2013;58:457–470. the Arrhenius Burn model and Pennes’ bio-heat equa- [12] Watanabe H, Kobayashi Y, Hashizume M, et al. tion for computation of soft tissue thermal damage. It Modeling the temperature dependence of thermo- physical properties: study on the effect of tempera- applies FEM and FDM to spatially and temporally dis- ture dependence for RFA. Conf Proc IEEE Eng Med cretize the Pennes’ bio-heat equation, which is solved Biol Soc. 2009;2009:5100–5105. with the acceleration of GPU hardware. Simulation [13] Mariappan P, Weir P, Flanagan R, et al. GPU-based results demonstrate that the proposed methodology RFA simulation for minimally invasive cancer treat- ment of liver tumours. Int J Cars. 2017;12:59–68. can efficiently predict soft tissue thermal damage. [14] Kroger € T, Altrogge I, Preusser T, et al. Numerical simu- Future research work includes Compute Unified lation of radio frequency ablation with state depend- Device Architecture (CUDA) implementation of the ent material parameters in three space dimensions. proposed methodology and consideration of thermal- Berlin, Heidelberg: Springer Berlin Heidelberg; 2006. mechanical response of soft tissues in the heating pro- [15] Rieder C, Kroger € T, Schumann C, et al. GPU-based cess of RFA. real-time approximation of the ablation zone for radiofrequency ablation. IEEE Trans Vis Comput Graph. 2011;17:1812–1821. Disclosure statement [16] Borsic A, Hoffer E, Attardo EA. GPU-accelerated real time simulation of radio frequency ablation thermal No potential conflict of interest was reported by the authors. dose. 2014 40th Annual Northeast Bioengineering Conference (Nebec); Boston, MA: IEEE; 2014: p. 1–2. [17] Chen G, Stang J, Haynes M, et al. Real-time three- References dimensional microwave monitoring of interstitial ther- mal therapy. IEEE Trans Biomed Eng. 2018;65: [1] Kok H, Kotte A, Crezee J. Planning, optimisation and 528–538. evaluation of hyperthermia treatments. Int J [18] Pennes HH. Analysis of tissue and arterial blood tem- Hyperthermia. 2017;33:593–607. peratures in the resting human forearm. J Appl [2] Zhang J, Zhong Y, Gu C. Deformable models for sur- Physiol. 1948;1:93–122. gical simulation: a survey. IEEE Rev Biomed Eng. 2018; [19] Zhang J, Zhong Y, Smith J, et al. Cellular neural net- 11:143–164. work modelling of soft tissue dynamics for surgical [3] Zhang J, Zhong Y, Smith J, et al. ChainMail based simulation. Technol Health Care. 2017;25:337–344. neural dynamics modeling of soft tissue deformation [20] Zhang JN, Zhong YM, Gu CF. Energy balance method for surgical simulation. Technol Health Care. 2017;25: for modelling of soft tissue deformation. Comput- 231–239. Aided Desi. 2017;93:15–25. [4] Lopresto V, Pinto R, Farina L, et al. Treatment plan- [21] Zhang J, Zhong Y, Smith J, et al. Energy propagation ning in microwave thermal ablation: clinical gaps and modeling of nonlinear soft tissue deformation for sur- recent research advances. Int J Hyperthermia. 2017; gical simulation. Simulation 2018;94:3–10. 33:83–100. [22] Zhang J, Zhong Y, Smith J, et al. Neural dynamics- [5] Prakash P. Theoretical modeling for hepatic micro- based Poisson propagation for deformable modelling. wave ablation. Open Biomed Eng J. 2010;4:27–38. Neural Comput Appl. 2019;31:1091–1101. [6] Shen W, Zhang J, Yang F. Modeling and numerical [23] Xu F, Lu TJ, Seffen KA. Biothermomechanics of skin simulation of bioheat transfer and biomechanics in tissues. J Mech Phys Solids. 2008;56:1852–1884. soft tissue. Math Comput Modell. 2005;41:1251–1265. [24] Zhang J, Zhong Y, Gu C. Ellipsoid bounding region- [7] Li X, Zhong Y, Subic A, et al. Prediction of tissue ther- based ChainMail algorithm for soft tissue deformation mal damage. Technol Health Care. 2016;24 Suppl 2: in surgical simulation. Int J Interact Des Manuf. 2018; S625–S629. 12:903–918. 12 J. ZHANG ET AL. [25] Zhang J, Zhong Y, Smith J, et al. A new ChainMail diffusion mechanics. Med Biol Eng Comput. 2018;12: approach for real-time soft tissue simulation. 2163–2176. Bioengineered 2016;7:246–252. [27] Zhang J, Shin J, Zhong Y, et al. Heat conduction-based [26] Zhang J, Zhong Y, Gu C. Soft tissue deformation methodology for nonlinear soft tissue deformation. Int modelling through neural dynamics-based reaction- J Interactive Des Manufacturing. 2018;13:147–161.
Journal
Computer Assisted Surgery – Taylor & Francis
Published: Oct 1, 2019
Keywords: Thermal damage; Arrhenius Burn integration; Pennes’ bio-heat equation; finite element method; Graphics Processing Unit; thermal ablation