Towards real-time finite-strain anisotropic thermo-visco-elastodynamic analysis of soft tissues for thermal ablative therapy
Jinao Zhang, Remi Jacob Lay, Stuart K. Roberts, Sunita Chauhan
[email protected]; [email protected] Department of Mechanical and Aerospace Engineering, Monash University, Wellington Road, Clayton, VIC 3800, Australia
Towards real-time finite-strain anisotropic thermo-visco-elastodynamic analysis of soft tissues for thermal ablative therapy
Jinao Zhang , Remi Jacob Lay , Stuart K. Roberts , Sunita Chauhan Department of Mechanical and Aerospace Engineering, Monash University, Clayton, Victoria, Australia Department of Gastroenterology, The Alfred Hospital, Melbourne, Victoria, Australia [Preprint]
Abstract.
Background and Objectives : Accurate and efficient prediction of soft tissue temperatures is essential to computer-assisted treatment systems for thermal ablation. It can be used to predict tissue temperatures and ablation volumes for personalised treatment planning and image-guided intervention. Numerically, it requires full nonlinear modelling of the coupled computational bioheat transfer and biomechanics, and efficient solution procedures; however, existing studies considered the bioheat analysis alone or the coupled linear analysis, without the fully coupled nonlinear analysis.
Methods : We present a coupled thermo-visco-hyperelastic finite element algorithm, based on finite-strain thermoelasticity and total Lagrangian explicit dynamics. It considers the coupled nonlinear analysis of (i) bioheat transfer under soft tissue deformations and (ii) soft tissue deformations due to thermal expansion/shrinkage. The presented method accounts for anisotropic, finite-strain, temperature-dependent, thermal, and viscoelastic behaviours of soft tissues, and it is implemented using GPU acceleration for real-time computation.
Results : The presented method can achieve thermo-visco-elastodynamic analysis of anisotropic soft tissues undergoing large deformations with high computational speeds in tetrahedral and hexahedral finite element meshes for surgical simulation of thermal ablation. We also demonstrate the translational benefits of the presented method for clinical applications using a simulation of thermal ablation in the liver.
Conclusion : The key advantage of the presented method is that it enables full nonlinear modelling of the anisotropic, finite-strain, temperature-dependent, thermal, and viscoelastic behaviours of soft tissues, instead of linear elastic, linear viscoelastic, and thermal-only modelling in the existing methods. It also provides high computational speeds for computer-assisted treatment systems towards enabling the operator to simulate thermal ablation accurately and visualise tissue temperatures and ablation zones immediately.
Keywords:
Thermal ablation; bioheat transfer; biomechanics; finite-strain thermo-visco-elastodynamics; surgical simulation. Introduction
Hyperthermia and thermal ablation treatments play an important role in the management of cancers following surgery, chemotherapy, and radiotherapy. Thermal ablation is a sub-category of hyperthermia in which local body tissues are heated until 50-100 โ over a few seconds or minutes to cause coagulation necrosis to damage and kill cancer cells [1]. The commonly used thermal ablative techniques include radiofrequency ablation (RFA), microwave ablation (MWA), laser ablation, and high-intensity focused ultrasound ablation (HIFU ablation). These techniques are typically minimally invasive (such as percutaneous RFA) or non-invasive (such as HIFU ablation), which are associated with fewer complications, shorter hospital stays, and less morbidity and mortality than those undergoing surgery. Percutaneous thermal tumour ablation under imaging guidance is an important and widely accepted treatment option, often given with curative intent particularly in patients with early-stage hepatocellular carcinoma (HCC, primary liver cancer). The two most commonly employed methods are RFA and MWA, where RFA is currently the most widely recommended first-line ablation technique for those not suitable for surgery [2]. Conventional RFA involves a monopolar electrode that induces local heating using alternating current. Tissue heating results in coagulative necrosis, but temperatures decrease both with distance from the electrode and with the presence of blood flow. The latter phenomenon results in the so-called โheat sink effectโ where adjacent large blood vessels cause tissue cooling, reducing the effectiveness of ablation [3]. MWA is a relatively recent and increasingly popular thermal ablative technique because of its ability to produce more rapid heating and higher maximum tissue temperatures. MWA has the advantages of producing wider and more predictable ablation volumes resulting in high complete ablation rates, and the ability to simultaneously treat multiple lesions [4] and potentially treat larger lesions more effectively [5]. Moreover, MWA requires less procedural time and is not subject to a heat-sink effect, which is known potential short-comings of RFA [6]. Rates of complications are low with both RFA (4.1%) and MWA (4.6%) [7, 8]. Still, despite the absence of convincing evidence, uptake of MWA is on the increase in many large academic centres in preference to RFA because of its ability to achieve a more rapid ablation, providing a workflow advantage particularly when performing multiple ablations. In the treatment of hepatic cancer, it was shown that the thermal ablation of solid, localised small ( โค ๐๐ ) tumours was very effective, achieved complete tumour necrosis and provided similar survival outcomes to those undergoing liver resection [9]. To assist the operator of thermal ablation with better patient safety and greater ablation accuracy, computer-integrated treatment system is becoming increasingly important to provide in silico analysis of patient-specific soft tissue temperatures and thermal damage [10]. The aim of achieving a safe and effective ablation is to induce a desired level of thermal injury to the target tumour while avoiding unintended thermal damage to nearby healthy tissues. To this end, our objectives are to employ computational bioheat transfer and biomechanics to enable the operator to simulate thermal ablation accurately and visualise tissue temperatures and ablation zones on patient-specific anatomical models immediately. The challenges in achieving these are two-fold: (i) the mathematical models must be able to describe soft tissue thermo-mechanical responses accurately, and (ii) the numerical solution procedures must be formulated efficiently for fast computation. The specific scenario we consider is towards real-time simulation of anisotropic thermo-visco-elastodynamic responses of soft tissues at finite-strains during thermal ablation. Due to body organ movements (such as respiratory motion) and tool-tissue interactions (such as percutaneous ablation), soft tissues are constantly in a state of motion consisted of dislocations and deformations; the former can be described by rigid-body displacements including translations and rotations, while the latter involves change in organ shapes that affect the distribution of thermal energy which is essential for computation of tissue temperatures. Furthermore, due to the variation of temperatures, the effects of thermal expansion [11] and shrinkage [12] induce thermal stresses, leading to thermally-induced deformations that also affect energy distributions. The thermally-induced deformations have also been evidenced in the studies of thermo-poroelastic modelling for fluid transport in tumours [13], thermal strain effects on low-density lipoprotein deposition in arterial walls [14], and thermo-diffusion effects [15]. Mathematically, these require full nonlinear modelling of the thermo-elastodynamic behaviour of soft tissues to account for (i) bioheat transfer under soft tissue deformations and (ii) soft tissue deformations due to thermal expansion/shrinkage. However, many of the existing studies were focused on thermal analysis alone without considering the mechanical behaviour [16-18]; therefore, these methods are incompatible and not suitable. Also, for such coupled nonlinear models, the need for efficient solution procedures for real-time computation on cost-effective and readily-available computing hardware for surgical simulation becomes more difficult [19]. Thermo-elastic analysis of soft tissues is not new to hyperthermia modelling; however, most of the existing methods are not suitable for surgical simulation. Key shortfalls of the existing methods stem from (i) expensive computational cost, (ii) solely steady-state analysis, and (iii) incompatibility with nonlinear thermo-elastic problems. Most of the reported methods [20, 21] were focused on sophisticated thermo-elastic models but failed to consider computational speed where solutions are often computationally expensive to obtain (e.g., Krรถger et al . [22] reported a simulation of an 8 minutes RFA in three-dimensional (3D) space took about 6 hours). Some studies were focused on solely the steady-state thermo-elastic analysis [23], which ignored the transient thermo-elastic response of soft tissues in dynamics. Finally, some works considered only homogeneous and isotropic material properties with linear thermo-elasticity [21, 23] or linear thermo-visco-elasticity [24]. As a result, these works did not describe anisotropic viscoelastic properties of soft tissues and are incompatible with the nonlinear problem of soft tissues undergoing large deformations where geometric and material nonlinearities are involved. To address the above issues, we develop a thermo-visco-hyperelastic total Lagrangian explicit dynamics finite element algorithm and utilise Graphics Processing Unit (GPU) acceleration for real-time finite-strain anisotropic thermo-visco-elastodynamic analysis of soft tissues. The following contributions are presented: (i) we develop a thermo-visco-hyperelastic model that is fully nonlinear to describe the nonlinear characteristics of bioheat transfer in deformed soft tissues with thermal expansion/shrinkage (Section 3), (ii) we develop efficient solution procedures for solving the coupled model in both tetrahedral and hexahedral computational grids (Section 4), (iii) we develop a GPU implementation for real-time computation (Section 5), and (iv) we demonstrate the translational benefits using a simulation of thermal ablation in the liver (Section 6.3). As the outcome, the presented method can achieve thermo-visco-elastic analysis of anisotropic soft tissues undergoing large deformations in tetrahedral and hexahedral finite element meshes with high computational speeds for surgical simulation of thermal ablation. The remainder of this work is organised as follows: Section 2 introduces the bioheat and biomechanics models, Sections 3, 4 and 5 present the proposed model, solution procedure, and GPU implementation, respectively. Algorithm verification, performance evaluation, and medical application are presented in Section 6. Discussions are presented in Section 7, and Section 8 concludes the present work. Bioheat and biomechanics models
Computational thermo-visco-elastodynamic analysis of soft tissues requires the definition of bioheat and biomechanics models. We use the Pennes model for bioheat transfer and the visco-hyperelastic model for biomechanics analysis.
Bioheat transfer model
Various bioheat transfer models were reported in literature to characterise heat transfer in biological soft tissues, including Pennes model [25], Wulff model [26], Klinger model [27], Chen and Holmes model [28], porous-media model [29], discrete vasculature model [30], and dual-phase-lag model [31] (see review papers [1, 32]), among those the Pennes model is widely used [33-35] and has been compared against experimental data [36, 37] to provide reliable tissue temperature predictions in the absence of large blood vessels (diameters larger than 0.5 ๐๐ ) [38]. Despite the development of other more sophisticated and rigorous bioheat models which may be seen as modifications from the classical Pennes model, they are usually obtained at the cost of increased numerical complexities and computational loads. Ge et al . [39] compared temperature profiles and histories between the Pennes model, dual-phase-lag model, and porous model, their results showed a lower temperature predicted by the porous model but nearly the same temperatures by the Pennes and the dual-phase-lag model. In addition, there exist controversies as to whether or not the dual-phase-lag conduction and any non-Fourier conduction is important for biological tissues [40] as some of the non-Fourier evidence has been called into question repeatedly [41, 42]. In light of this, we employ the Pennes model in the present work, which has proven to be remarkably effective for bioheat transfer analysis [43]; it is also computationally efficient and consumes less data storage space. The Pennes bioheat model accounts for the thermal effects of anisotropic heat conduction in solid tissues, isotropic blood perfusion, metabolic heat generation, and regional heat sources. The governing equation of the transient Pennes model in 3D is given by [25] ๐ ๐ก ๐ ๐ก ๐ ๐(๐ฑ) ๐ก ๐๐ก = โ โ ( ๐ ๐ก โ ๐(๐ฑ) ๐ก ) โ ๐ค ๐๐ก ๐ ๐๐ก ( ๐(๐ฑ) ๐ก โ ๐ ๐๐ก ) + ๐ ๐๐ก + ๐ ๐๐ก โ๐ฑ โ ๐บ (1) where ๐ is the tissue mass density, ๐ the specific heat capacity, ๐(๐ฑ) the temperature at a material point ๐ฑ(๐ฅ, ๐ฆ, ๐ง) in the continuum ๐บ , ๐ก the time, ๐ the thermal conductivity, ๐ค ๐ the blood perfusion rate, ๐ ๐ the blood specific heat capacity, ๐ ๐ the arterial blood temperature, ๐ ๐ the metabolic heat generation rate, ๐ ๐ the regional heat source rate, โ โ the divergence operator, and โ the gradient operator; a left superscript โข ๐ก denotes the system configuration in which a quantity occurs ( ๐ก denotes the current time system configuration); the material properties are, in general, temperature-dependent [44, 45] and can be directly incorporated for tissue temperature analysis. Soft tissue temperatures must be determined from the deformed configuration due to dynamic motions of soft tissues; however, the classical Pennes model is based on only the static non-moving state [46]. To address this, we apply a mapping function to transform the unknown current (deformed) configuration to the known reference (undeformed) configuration for bioheat transfer analysis, which will be demonstrated in Section 3.1. Visco-hyperelastic biomechanics model
Soft tissues undergo large deformations under external forces where geometric and material nonlinearities are involved, and they exhibit time-dependent viscous and elastic behaviours [47]. The widely used linear elastic and linear viscoelastic models are based on the assumptions of infinitesimally small deformations with the linear stress-strain relationship, but they are not compatible with large deformations of soft tissues. For fully compatible nonlinear modelling, we use the visco-hyperelastic biomechanics model based on finite elasticity for which nonlinear continuum mechanics is the fundamental basis. By referring to the reference system configuration, the variation of strain at a material point is expressed by the Green-Saint Venant strain tensor ๐ = ( ๐ โ ๐) where ๐ is the right Cauchy-Green tensor and ๐ is the identity matrix of the second rank. The corresponding strain energy density is expressed by ๐น which is typically a function of the strain invariants of ๐ . Accordingly, the energy conjugated stress measure is the second Piola-Kirchhoff stress tensor ๐ = ๐ ๐น( ๐ ) ๐ ๐ . Mechanical behaviours are accommodated in ๐น . For anisotropic soft tissues exhibiting directional-dependent behaviours, ๐น can be modified by employing unit vectors to describe local fibre directions [48]. For viscoelastic behaviours, ๐น can be modified by employing a time-dependent ๐นฬ such as ๐นฬ = โซ ๐(๐ก โ ๐ก โฒ ) ๐ ๐น ๐๐ก โฒ ๐๐ก โฒ๐ก0 in a convolution integral [49] where ๐(๐ก โ ๐ก โฒ ) = ๐ โ +โ ๐ ๐ ๐ (๐ก โฒ โ๐ก)/๐ ๐ ๐๐=1 is the relaxation function of a generalised Maxwell model [50], and ๐ โ , ๐ ๐ , and ๐ ๐ are positive constants. By referring to the reference configuration, the equation of motion is given by [51] โ โ ๐ + ๐ = ๐ ๐ฬ ๐ก (2) where ๐ is the first Piola-Kirchhoff stress tensor, ๐ the reference body force, and ๐ฬ the material acceleration field; a left subscript โข denotes the configuration with respect to which the quantity is measured (0 denotes the reference configuration). The above balance principle is valid for any arbitrary parts of a continuum, but the stress measure is expressed for elastodynamic analysis. For thermo-elastodynamic analysis, the constitutive expression for the stress is obtained based on the Helmholtz free energy which is a function of the finite-strain and temperature and by exploring the conservation of energy and the second law of thermodynamics [52]. To obtain this stress expression, we apply a multiplicative decomposition [52] to include the elastic and thermal deformations for the total stress, which will be demonstrated in Section 3.2. Coupling of bioheat and biomechanics models
A two-way interaction of bioheat and biomechanics models is identified: (i) due to soft tissue deformations, transient tissue temperatures must be computed based on the deformed configuration, and (ii) thermal expansion/shrinkage also induces soft tissue deformations. Both (i) and (ii) affect the prediction of thermal ablation outcomes and require the coupled thermo-elastodynamic modelling. Fig. 1 illustrates the notations used in the proposed coupled thermo-elastodynamic analysis. Using the Lagrangian description, soft tissues undergo deformations ๐บ โ ๐บ ๐ก described by the deformation gradient ๐ . An intermediate system configuration ๐บ ๐ is introduced, decomposing the total deformation gradient ๐ into ๐ ( ๐บ โ ๐บ ๐ ) and ๐ ๐๐ก ( ๐บ ๐ โ ๐บ ๐ก ) . We compute (i) bioheat transfer under deformations by transformation of the system configuration ๐บ ๐ก โ ๐บ (Section 3.1) and (ii) soft tissue deformations due to thermal expansion/shrinkage by multiplicative decomposition ๐ = ๐ ๐๐ก ๐ (Section 3.2). Fig. 1 Reference ๐บ , intermediate ๐บ ๐ , and current ๐บ ๐ก configurations, and the deformations described by ๐ ( ๐บ โ ๐บ ๐ก ) , ๐ ( ๐บ โ ๐บ ๐ ) , and ๐ ๐๐ก ( ๐บ ๐ โ ๐บ ๐ก ) . Mathematical modelling of bioheat transfer under soft tissue deformations
We employ a transformation of system configurations to achieve bioheat computation on deformed soft tissues [53]. The mapping function for deformations ๐บ โ ๐บ ๐ก can be defined by ๐ ๐ก โถ ๐บ ร โ โ ๐บ ๐ก , ๐ฑ ๐ก = ๐ ๐ก ( ๐ฑ ) = ๐( ๐ฑ , ๐ก) (3) where the mapping function ๐ is continuous in time, continuously differentiable, globally invertible, and orientation-preserving at all times. Using the current (deformed) configuration ๐บ ๐ก , the Pennes model can be integrated over the volume to get โซ ( ๐ ๐ก ๐ ๐ก ๐ ๐(๐ฑ) ๐ก ๐๐ก ) ๐บ ๐ก ๐ ๐บ ๐ก = โซ ( โ ๐ก โ ( ๐ ๐ก โ ๐ก ๐(๐ฑ) ๐ก ) โ ๐ค ๐๐ก ๐ ๐๐ก ( ๐(๐ฑ) ๐ก โ ๐ ๐๐ก ) + ๐ ๐๐ก + ๐ ๐๐ก ) ๐บ ๐ก ๐ ๐บ ๐ก (4) By performing a change of variables, it can be expressed in terms of the reference configuration ๐บ as โซ (โฆ ) ๐บ ๐ก ๐ ๐บ ๐ก = โซ (โฆ )det (โ๐( ๐ฑ , ๐ก)) ๐บ ๐ ๐บ = โซ (โฆ )det (๐ ๐ฑ ๐ก ๐ ๐ฑ ) ๐บ ๐ ๐บ = โซ (โฆ )det( ๐ ) ๐บ ๐ ๐บ (5) where โ is the differential operator. The divergence and gradient operators ( โ ๐ก โ and โ ๐ก ) are expressed in terms of ๐บ ๐ก which is unknown. By applying the chain rule, they can be expressed in terms of the known reference configuration ๐บ as โซ โ ๐ก โ ( ๐ ๐ก โ ๐ก ๐(๐ฑ) ๐ก ) ๐บ ๐ก ๐ ๐บ ๐ก = โซ ( โ ๐
0๐ก โ1 โ ( ๐ ๐ก โ ๐
0๐ก โ1
๐(๐ฑ) ๐ก )) det( ๐ ) ๐บ ๐ ๐บ (6) The above formulation allows for the computation of bioheat transfer on deformed configurations of soft tissues, compared to the classical Pennes model which was computed based on only the static non-moving state [46]. From the point of view of numerical efficiency, the above formulation is based on only ๐ ( ๐บ โ ๐บ ๐ก ) , ๐ ๐ก , and ๐(๐ฑ) ๐ก , whereas the other variables are defined over the reference configuration ๐บ that can be precomputed. Mathematically, it accounts for nonlinear characteristics of bioheat transfer at finite strains. Mathematical modelling of soft tissue deformations due to thermal expansion/shrinkage
We apply a multiplicative decomposition to build the stress expression for thermo-elastodynamic analysis. Since the deformations caused by thermal expansion and shrinkage are mathematically equivalent, we will use thermal expansion for the volume change of soft tissues in the remainder of the present work. Prior to the onset of thermal expansion, we consider the tissue is in a physiological state described by the stress-free intermediate configuration ( ๐บ ๐ in Fig. 1), which is a hypothetical configuration obtained by isothermal elastic destressing of the deformed configuration ๐บ ๐ก to zero stress [54]. With this intermediate configuration, the total deformation gradient ๐ for ๐บ โ ๐บ ๐ก can be decomposed into a thermal part ๐ ๐กโ๐๐0๐ and an elastic part ๐ ๐๐๐๐ ๐๐ก via a local multiplicative decomposition [52] expressed by ๐ = ๐ ๐๐๐๐ ๐๐ก ๐ ๐กโ๐๐0๐ (7) To meet the requirement for zero stress under zero strain, the constitutive expression must be written in terms of the elastic deformation gradient ๐ ๐๐๐๐ ๐๐ก , yielding ๐ ๐๐๐๐ ๐๐ก = ๐ ๐ ๐กโ๐๐0๐ โ1 (8) The second Piola-Kirchhoff stress ๐ ๐๐ก associated with the stress-free intermediate configuration ๐บ ๐ can then be expressed by ๐ ๐๐ก = ๐ ๐น( ๐ ๐๐ก ) ๐๐ก ๐ ๐ ๐๐ก = 2 ๐ ๐น( ๐ ๐๐ก ) ๐๐ก ๐ ๐ ๐๐ก = ๐ ๐๐ก ( ๐ ๐๐๐๐ ๐๐ก ) = ๐ ๐๐ก ( ๐ ๐ ๐กโ๐๐0๐ โ1 ) (9) By pulling back the stress to the reference configuration ๐บ , the total stress ๐ can be expressed by ๐ = det( ๐ ๐กโ๐๐0๐ ) ๐ ๐กโ๐๐0๐ โ1 ๐ ๐๐ก ( ๐ ๐ ๐กโ๐๐0๐ โ1 ) ๐ ๐กโ๐๐0๐ โ๐ (10) The above formulation is based on the deformation gradients ๐ and ๐ ๐กโ๐๐0๐ , which enables thermal expansion to be included entirely in the material model. The thermal deformation gradient ๐ ๐กโ๐๐0๐ can be specified depending on the type of material anisotropy. For an orthotropic thermal expansion with the principal axes parallel to unit vectors ๐ฆ , ๐ง and ๐ฆ ร ๐ง in ๐บ , ๐ ๐กโ๐๐0๐ can be specified by [54] ๐ ๐กโ๐๐0๐ = ๐ ๐ ๐ + (๐ ๐ โ ๐ ๐ ) ๐ฆ โ ๐ฆ + (๐ ๐ โ ๐ ๐ ) ๐ง โ ๐ง (11) where ๐ ๐ , ๐ ๐ and ๐ ๐ are the stretch ratios in the orthogonal directions ๐ฆ ร ๐ง , ๐ฆ and ๐ง , respectively; the above expression can be modified for transversely isotropic thermal expansion ๐ ๐กโ๐๐0๐ = ๐ ๐ ๐ + (๐ ๐ โ ๐ ๐ ) ๐ฆ โ ๐ฆ and isotropic thermal expansion ๐ ๐กโ๐๐0๐ = ๐ ๐ ๐ . The thermal stretch ratio ๐ is related to the thermal expansion coefficient ๐ผ by ๐ โ 1 +๐ผ( ๐(๐ฑ) ๐ก โ ๐(๐ฑ) ) [52]. Numerical solution procedures
We develop a thermo-visco-hyperelastic total Lagrangian explicit dynamics finite element algorithm for solutions of the coupled bioheat and biomechanics analysis.
Thermo-visco-hyperelastic total Lagrangian explicit dynamics finite element algorithm
To achieve fast and accurate numerical solutions, we develop solution procedures based on the (i) fast explicit dynamics finite element algorithm [55] for transient bioheat computation [56] under soft tissue deformations [53] and the (ii) total Lagrangian explicit dynamics finite element algorithm [57] for biomechanics computation of large-strain soft tissue deformations. The proposed formulation has not been evidenced in previous works. We employ the total Lagrangian formulation for computational biomechanics for its numerical efficiency without error accumulation. Displacement-based finite element formulations use the updated Lagrangian formulation or the total Lagrangian formulation to describe the frame of reference [58]. The updated Lagrangian obtains variable values at the current configuration based on the end of the previous time step via an incremental strain-description, which requires a re-calculation of spatial derivatives at each time step. In contrast, the total Lagrangian considers all variables referred to the reference configuration so that the strains lead to correct results after a load cycle without error accumulation [59], which is not the case in the incremental strain-description. More importantly, it allows for the spatial derivatives to be precomputed and stored [57]. It was shown that the total Lagrangian took 2.1 ๐๐ for a time step to obtain the solution for an ellipsoid indentation, which was about five times faster compared to 10.6 ๐๐ consumed by the updated Lagrangian [60]. We employ explicit finite element method (FEM) and explicit time integration for computational bioheat transfer and biomechanics for numerical efficiency and parallelisation. Explicit FEM [61] is a simplified form of FEM that is often used in surgical simulation. The internal and external loads and mass are lumped to the nodes, leading to block diagonal mass and damping matrices to enable computation to be performed at the element level without the need for assembling the stiffness matrix for the entire model [57]. The explicit FEM can be integrated in time either explicitly or implicitly for dynamics, and we use explicit time integration to allow for variable values in future states to be computed directly based on the current state without the need for stiffness matrix inversion which is required in the implicit integration. It was shown that the implicit integration consumed computation time that was at least one order of magnitude larger than that by the explicit counterpart [62]. The resultant explicit dynamics leads to an explicit formulation for unknown state variables that is well suited for parallel computation. The global system of equations can be split into independent equations for individual nodes, allowing each nodal computation to be assigned to a parallel task to perform calculations independently. Along with computationally efficient low-order finite elements, the above formulations lead to very efficient computation at run-time. Formulation for bioheat transfer under soft tissue deformations
Using the first-order explicit time integration and the lumped (diagonalised) mass approximation while ensuring the conservation of mass, the discretized matrix equation of the Pennes model can be written as ๐ ๐๐๐๐๐ก ( ๐(๐ฑ) ๐ก+โ๐ก โ ๐(๐ฑ) ๐ก โ๐ก ) = ๐ ๐ก ๐(๐ฑ) ๐ก โ ๐ ๐๐๐๐๐๐ก ๐(๐ฑ) ๐ก + ๐ ๐๐ก + ๐ ๐๐ก + ๐ ๐๐ก (12) where ๐ ๐๐๐๐ is the diagonalised thermal mass (mass and specific heat capacity) matrix, ๐(๐ฑ) the vector of nodal temperatures, โ๐ก the time step, ๐ the thermal stiffness (conduction) matrix, ๐ ๐๐๐๐๐ the diagonalised thermal stiffness (blood perfusion) matrix, ๐ ๐ , ๐ ๐ , and ๐ ๐ the vectors of heat flow of blood perfusion at arterial temperature, metabolic heat generation, and regional heat sources, respectively. The equation can be further rearranged into ๐(๐ฑ) ๐ก+โ๐ก = ๐(๐ฑ) ๐ก + โ๐ก ๐ ๐๐๐๐๐ก โ1 (โ ๐ ๐๐กโ๐๐0๐ก๐ โ ๐ ๐๐๐๐๐๐ก ๐(๐ฑ) ๐ก + ๐ ๐๐ก + ๐ ๐๐ก + ๐ ๐๐ก ) (13) where โ ๐ ๐๐กโ๐๐0๐ก๐ = ๐ ๐ก ๐(๐ฑ) ๐ก = ๐ ๐กโ๐๐0๐ก (14) where ๐ ๐๐กโ๐๐ is the vector of thermal load components due to heat conduction in an element ๐ , and ๐ ๐กโ๐๐ is the vector of global nodal thermal loads. ๐ ๐๐กโ๐๐0๐ก is computed by ๐ ๐๐กโ๐๐0๐ก = โซ โ ๐ก ๐ฑ๐ก ๐ ๐ ๐๐ก โ ๐ก ๐ฑ๐ก ๐ ๐ ๐ก๐ ๐๐ก ๐(๐ฑ) ๐ก (15) where โ๐ก ๐ฑ is the matrix of element shape function spatial derivatives, ๐ ๐ the element thermal conductivity matrix, and ๐ ๐ the element volume. By incorporating the proposed bioheat formulation under soft tissue deformations (Eq. (6)), the above equation at the current configuration ๐บ ๐ก can be written in terms of the reference configuration ๐บ as ๐ ๐๐กโ๐๐0๐ก = โซ (โ ๐ก ๐ฑ0 ๐
0๐ก โ1 ) ๐ ๐ ๐๐ก (โ ๐ก ๐ฑ0 ๐
0๐ก โ1 )det( ๐ ) ๐ ๐ ๐0 ๐(๐ฑ) ๐ก (16) where ๐ can be computed from the element nodal displacements matrix ๐ฎ ๐ by ๐ = ( ๐ฎ ๐๐ก ) ๐ โ ๐ก ๐ฑ0 + ๐ (17) For fast computation, the computationally efficient 3D low-order finite elements such as the eight-node reduced-integrated hexahedral elements and the four-node linear tetrahedral elements are used, leading to the following computations at the element level: (i) eight-node reduced-integrated hexahedral element: ๐ ๐๐กโ๐๐0๐ก = 8det( ๐ )(โ ๐ก ๐ฑ0 ๐
0๐ก โ1 ) ๐ ๐ ๐๐ก โ ๐ก ๐ฑ0 ๐
0๐ก โ1 det( ๐ ) ๐(๐ฑ) ๐ก (18) (ii) four-node linear tetrahedral element: ๐ ๐๐กโ๐๐0๐ก = ๐ (โ ๐ก ๐ฑ0 ๐
0๐ก โ1 ) ๐ ๐ ๐๐ก โ ๐ก ๐ฑ0 ๐
0๐ก โ1 det( ๐ ) ๐(๐ฑ) ๐ก (19) where ๐ is the element Jacobian matrix. At each time step, the temperature ๐ (๐) (๐ฑ) ๐ก+โ๐ก at a node ๐ can be obtained by ๐ (๐) (๐ฑ) ๐ก+โ๐ก = ๐ (๐) (๐ฑ) ๐ก + โ๐ก ๐ถ (๐)๐๐๐๐๐ก โ1 ( ๐ (๐)๐กโ๐๐0๐ก โ ๐พ ๐(๐)๐๐๐๐๐ก ๐ (๐) (๐ฑ) ๐ก + ๐ ๐(๐)๐ก + ๐ ๐(๐)๐ก + ๐ ๐(๐)๐ก ) (20) The above formulation allows for the computation of bioheat transfer with respect to the deformed configuration of soft tissues. Eq. (20) states an explicit formulation for the unknown temperature ๐ (๐) (๐ฑ) ๐ก+โ๐ก which can be obtained based on the variable values at ๐ก only; hence, the global system of equations can be split into independent equations for parallel computation. Furthermore, there is no need for iterations during any process of the algorithm. Temperature-dependent thermal properties and nonlinear thermal boundary conditions can also be directly incorporated. Formulation for soft tissue deformations due to thermal expansion
The discretized matrix equation of the mechanics of motion can be written as ๐ ๐๐๐๐๐ก ๐ ๐(๐ฑ) ๐ก ๐๐ก + ๐ ๐๐๐๐๐ก ๐ ๐(๐ฑ) ๐ก ๐๐ก + ๐(๐) ๐ก ๐(๐ฑ) ๐ก = ๐ ๐ก (21) where ๐ ๐๐๐๐ is the diagonalised mass matrix, ๐ ๐๐๐๐ = ๐พ๐ ๐๐๐๐ the mass-proportional damping (a special case of Rayleigh damping, ๐พ the damping coefficient), ๐(๐) the stiffness matrix,
๐(๐ฑ) the vector of nodal displacements, and ๐ the vector of externally applied forces. The displacements ๐(๐ฑ) ๐ก+โ๐ก at the next time step can be computed based on the explicit central-difference scheme [63], yielding
๐(๐ฑ) ๐ก+โ๐ก = ๐ ๐ก โ โ ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก๐ + 2 ๐ ๐๐๐๐๐ก โ๐ก ๐(๐ฑ) ๐ก + ( ๐ ๐๐๐๐๐ก ๐๐๐๐๐ก โ๐ก ) ๐(๐ฑ) ๐กโโ๐ก ๐ ๐๐๐๐๐ก ๐๐๐๐๐ก โ๐ก (22) where โ ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก๐ = ๐(๐) ๐ก ๐(๐ฑ) ๐ก = ๐ ๐กโ๐๐โ๐๐๐๐ 0๐ก (23) where ๐ ๐๐กโ๐๐โ๐๐๐๐ is the vector of nodal forces due to thermal and elastic stresses in an element ๐ , and ๐ ๐กโ๐๐โ๐๐๐๐ is the vector of global nodal forces. ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก can be computed by ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก = โซ ๐ ๐ โ ๐ก ๐ฑ0 ๐ ๐ ๐0 (24) where ๐ is the total stress to account for the effect of thermal expansion given by Eq. (10). The corresponding low-order finite element formulations are: (i) eight-node reduced-integrated hexahedral element: ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก = 8det( ๐ ) ๐ (det( ๐ ๐กโ๐๐0๐ ) ๐ ๐กโ๐๐0๐ โ1 ๐ ๐๐ก ( ๐ ๐ ๐กโ๐๐0๐ โ1 ) ๐ ๐กโ๐๐0๐ โ๐ ) โ ๐ก ๐ฑ0 (25) (ii) four-node linear tetrahedral element: ๐ ๐๐กโ๐๐โ๐๐๐๐ 0๐ก = ๐ ๐ (det( ๐ ๐กโ๐๐0๐ ) ๐ ๐กโ๐๐0๐ โ1 ๐ ๐๐ก ( ๐ ๐ ๐กโ๐๐0๐ โ1 ) ๐ ๐กโ๐๐0๐ โ๐ ) โ ๐ก ๐ฑ0 (26) At each time step, the displacement component ๐ข ๐ฅ(๐)๐ก+โ๐ก of ๐ฎ (๐)๐ก+โ๐ก ( ๐ข ๐ฅ(๐)๐ก+โ๐ก , ๐ข ๐ฆ(๐)๐ก+โ๐ก , ๐ข ๐ง(๐)๐ก+โ๐ก ) at a node ๐ can be obtained by ๐ข ๐ฅ(๐)๐ก+โ๐ก = ๐ ๐ฅ(๐)๐ก โ ๐ ๐ฅ(๐)๐กโ๐๐โ๐๐๐๐ 0๐ก + 2 ๐ ๐ฅ(๐)๐๐๐๐๐ก โ๐ก ๐ข ๐ฅ(๐)๐ก + ( ๐ท ๐ฅ(๐)๐๐๐๐๐ก ๐ฅ(๐)๐๐๐๐๐ก โ๐ก ) ๐ข ๐ฅ(๐)๐กโโ๐ก ๐ท ๐ฅ(๐)๐๐๐๐๐ก ๐ฅ(๐)๐๐๐๐๐ก โ๐ก (27) For viscoelastic modelling where the time-dependent ๐นฬ = โซ ๐(๐ก โ ๐ก โฒ ) ๐ ๐น ๐๐ก โฒ ๐๐ก โฒ๐ก0 is used (Section 2.2), the corresponding time-dependent stress is expressed by ๐ฬ = โซ ๐(๐ก โ ๐ก โฒ ) ๐ ๐ ๐๐ก โฒ ๐๐ก โฒ๐ก0 , and it can be computed by ๐ฬ = ๐ โ โ ๐
0๐ก ๐๐๐=1 [64] where ๐
0๐ก ๐ = โ๐ก๐ ๐ โ๐ก+๐ ๐ ๐ + ๐ ๐ โ๐ก+๐ ๐ ๐ . The corresponding thermo-visco-elastic nodal forces ๐ ๐๐กโ๐๐โ๐ฃ๐๐ ๐๐โ๐๐๐๐ 0๐ก can be obtained by substituting the time-dependent ๐ฬ for ๐ in Eq. (24), yielding ๐ ๐๐กโ๐๐โ๐ฃ๐๐ ๐๐โ๐๐๐๐ 0๐ก = โซ ๐ ๐ฬ โ ๐ก ๐ฑ0 ๐ ๐ ๐0 (28) which is used to replace ๐ ๐ฅ(๐)๐กโ๐๐โ๐๐๐๐ 0๐ก by ๐ ๐ฅ(๐)๐กโ๐๐โ๐ฃ๐๐ ๐๐โ๐๐๐๐ 0๐ก in Eq. (27) for viscoelastic modelling. The above formulation enables the time-dependent effect to be included entirely in the thermo-elastic stress for thermo-visco-elastic constitutive analysis. The constitutive material law, temperature-dependent mechanical properties and nonlinear mechanical boundary conditions can be directly incorporated. Fig. 2 illustrates an implementation of the proposed numerical solution procedures. Fig. 2. Key components of the proposed numerical solution procedures. GPU implementation
Having established the formulations for bioheat transfer under soft tissue deformations (Section 4.2) and soft tissue deformations due to thermal expansion (Section 4.3) for tetrahedral and hexahedral computational grids, the presented method is implemented using GPU parallel computation for real-time surgical simulation. The GPU implementation consists of a host (CPU) and a device (GPU) code for host-device interaction and device parallel computation. The host implementation is responsible for precomputing and storing the constant matrices and simulation parameters and for invoking device methods to interact with GPUs for memory allocation, texture binding, data copy from/to device, and kernel launching. As mentioned previously, many simulation parameters can be precomputed owing to the total Lagrangian formulation, such as the spatial derivatives, initial element volumes, initial Jacobian, and mass and damping matrices. The host code is written in the C++ programming language in Visual Studio 2017. The device implementation is responsible for computing new internal thermal loads, temperatures, internal forces and displacements for the temperature and displacement fields. The nodal loads, such as the thermal loads and internal forces, are computed by an โelementโ kernel for tetrahedron/hexahedron calculations. The nodal temperatures and displacements are computed by Eqs. (20) and (27), respectively, at each time step for time-stepping, and they are computed by a โnodeโ kernel for node calculation. A simulation time step is achieved by launching the โelementโ kernel across ๐ ๐ threads and the โnodeโ kernel across ๐ ๐ threads where ๐ ๐ and ๐ ๐ denote the numbers of finite elements and nodes in the organ models, respectively. The GPU-accelerated computation is achieved via a time-loop of such time-step computation. The device code is written using the NVIDIA CUDA programming API version 10.2. Numerical results
Numerical evaluations are conducted to assess the validity of the presented methodology for simulating large strain thermo-elastodynamics of soft tissues for surgical simulation, and assessments on CPU and GPU computational performance are presented. The translational benefits for a clinically relevant application are demonstrated using a simulation of thermal ablation in the liver.
Algorithm verification
Although living tissues need to be used for validation of the proposed method, experiments with living tissues raise ethical issues with technical difficulties. In the present work, a 3D numerical example with thermal and mechanical properties similar to those of soft tissues is used for algorithm verification, which offers a controlled environment for a quantitative assessment on numerical accuracy.
Fig. 3 illustrates the geometry and simulation settings of the employed model, and Table. 1 presents the values of the geometry, material properties and simulation time. The prescribed displacement ๐ข ๐ง = 0.02 ๐ in Fig. 3 corresponds to 40% extension compared to the height (0.05 ๐ ) of the model for large deformations. The right face of the model was assumed to be fixed in position during the simulation. Adiabatic boundary condition was applied to the exterior of the model. The established nonlinear procedures, โDynamic, Temp-disp, Explicitโ, from commercial finite element analysis package, ABAQUS/CAE 2018 (2017_11_08-04.21.41 127140), were used to produce reference solutions under the same conditions for comparison. To quantitatively evaluate the validity of the presented method, normalised relative errors ๐๐ ๐ธ ๐ = | ๐ ๐๐ด๐ต๐ด๐๐๐ โ๐ ๐๐๐๐๐๐๐ ๐๐ ๐ ๐๐๐ฅ๐ด๐ต๐ด๐๐๐ โ๐ ๐๐๐๐ด๐ต๐ด๐๐๐ | and ๐๐ ๐ธ ๐ข = | ๐ข ๐๐ด๐ต๐ด๐๐๐ โ๐ข ๐๐๐๐๐๐๐ ๐๐ ๐ข ๐๐๐ฅ๐ด๐ต๐ด๐๐๐ โ๐ข ๐๐๐๐ด๐ต๐ด๐๐๐ | and total error indicators ๐ ๐ = โ โ (๐ ๐๐ด๐ต๐ด๐๐๐ โ๐ ๐๐๐๐๐๐๐ ๐๐ ) โ (๐ ๐๐ด๐ต๐ด๐๐๐ ) and ๐ ๐ข = โ โ (๐ข ๐๐ด๐ต๐ด๐๐๐ โ๐ข ๐๐๐๐๐๐๐ ๐๐ ) โ (๐ข ๐๐ด๐ต๐ด๐๐๐ ) were used; ๐ข ๐ included the three components ๐ข ๐,๐ฅ , ๐ข ๐,๐ฆ and ๐ข ๐,๐ง . Fig. 3. Model geometry and simulation settings for algorithm verification. Table. 1 Simulation settings of the model in Fig. 3. Geometry Values Unit Description Nodes 4489 17956 degrees of freedom (DOFs, ๐ฅ, ๐ฆ, ๐ง, ๐ per node) Elements 24290 Four-node linear tetrahedrons Material properties Density ๐ [๐๐/๐ ] [21] Specific heat capacity ๐ [๐ฝ/(๐๐ โ โ)] [21] Thermal conductivity ๐ [๐/(๐ โ โ)] [21] Thermal expansion ๐ผ [1/โ] Isotropic neo-Hookean hyperelastic model
๐น = ๐2 (๐ผ ฬ โ 3) + ๐ 2 (๐ฝ โ 1) Shear modulus ๐ [๐๐] Corresponding to Youngโs modulus 3500 Pa and Poissonโs ratio 0.47 [53] Bulk modulus ๐ [๐๐] Simulation time 5 ๐ Time step 0.00008 ๐ (62500 steps) Note: the anisotropic thermal conductivity, anisotropic hyperelastic material models, and viscoelasticity have been previously verified in Refs. [55] and [64]. Fig. 4 presents a comparison between four cases of the numerical example for validation: (a) thermal analysis only, (b) mechanical analysis only, (c) thermo-mechanical analysis without and (d) with thermal expansion. Cases (a-c) can be achieved by setting ๐ข ๐ง , ๐ ๐ and ๐ ๐ , and ๐ผ to zero, respectively. The results show good agreement between the proposed method and ABAQUS for all four cases, demonstrating the validity of the presented method in conducting individual and full analysis of finite-strain thermo-elastodynamics of soft tissues. Compared to the temperature field in (a), the temperature distribution in (c) was affected by the deformations and was further affected by thermal expansion in (d) which had a broader distribution of temperatures and a lower maximum temperature at the heat source centre. Compared with the displacement fields in (b) and (c), the displacement field in (d) had considerable local expansions at the three heat source regions that compressed the cylindrical holes and yielded less volume deformations. A statistical summary is presented in Table. 2. (a) (b) (c) (d) T e m p e r a t u r e c o m p a r i s on D i s p l ace m e n t c o m p a r i s on Fig. 4. Comparisons of temperatures, displacements, and numerical errors between four cases of the numerical example: (a) thermal analysis, (b) mechanical analysis, (c) thermo-mechanical analysis without and (d) with thermal expansion; rows ๐ฅ๐ง plane, ๐ฆ๐ง plane, and ๐ฅ๐ฆ plane, respectively. Table. 2 A statistical summary of the comparisons in Fig. 4. (a) (b) (c) (d) Thermal Mechanical Thermo-mechanical without thermal expansion Thermo-mechanical with thermal expansion P A P A P A P A T e m p e r a t u r e c o m p a r i s on ( โ ) Min 37 37 37 37 37 37 37 37 Max 47.2238 47.2236 37 37 47.2085 47.2084 47.1112 47.1111 Median 37 37 37 37 37 37 37 37 Q1 37 37 37 37 37 37 37 37 Q3 37.4252 37.4252 37 37 37.4001 37.4002 37.7441 37.7437 ๐ ๐ D i s p l ace m e n t c o m p a r i s on ( ๐ ) Min 0 0 -0.0045 -0.0045 -0.0045 -0.0045 -0.0035 -0.0035 Max 0 0 0.0200 0.0200 0.0200 0.0200 0.0212 0.0212 Median 0 0 3.6013e-4 3.6046e-4 3.6013e-4 3.6046e-4 3.2390e-4 3.2365e-4 Q1 0 0 -1.2742e-4 -1.2665e-4 -1.2742e-4 -1.2665e-4 -1.4113e-4 -1.4104e-4 Q3 0 0 0.0037 0.0037 0.0037 0.0037 0.0038 0.0038 ๐ ๐ข
0 7.3887e-4 7.3887e-4 6.3711e-4 P: Proposed, A: ABAQUS
Fig. 5 presents the effect of thermal expansion on the model temperature and displacement fields. As shown in Fig. 5(a-b), higher temperature regions were observed around the outer expanded region of the three heat sources, and even higher temperatures were observed between the heat sources due to the concentration of heat by the expansion of heat source regions. On the other hand, lower temperatures were observed around the middle section of the expanded regions, where there was less heat compared to the case without thermal expansion. As shown in Fig. 5(c-d), the effect of thermal expansion induced noticeable deformations around the heat source regions, and even greater deformations were observed at the cylindrical hole regions. On the other hand, fewer deformations were observed between heat sources as the regions expanded in opposite directions, resulting in less and close to zero resultant displacements. (a) (c) (b) (d) Fig. 5. The effect of thermal expansion on the temperature and displacement fields; colour maps indicate the difference in temperature and displacement values with and without thermal expansion.
CPU and GPU computational performance
The presented method was evaluated on an Intel(R) Core(TM) i5-2500K CPU @ 3.30 GHz with 8.0 GB RAM PC using Windows 10 operating system. Fig. 6(a) presents CPU solution times for a single time step in four-node linear tetrahedral (T4) and eight-node reduced-integrated hexahedral (H8) meshes for temperature-independent thermo-mechanical (TherMechTI), thermo-mechanical with thermal expansion (TherMechExpanTI), and temperature-dependent thermo-mechanical with thermal expansion (TherMechExpanTD) at 11 different model sizes. It can be seen that the computation times of TherMechTI were the least for both T4 and H8, followed by those of TherMechExpanTI and TherMechExpanTD. Furthermore, it was noticed that the solution times increased almost linearly with the increase of model sizes, and linear interpolation and extrapolation could be used to determine the solution times at unknown model sizes. Fig. 6(b) presents the computational overheads in CPU computation with the averaged overheads given in Table. 3. The computation time of the hourglass control algorithm [60] for H8 was considered in all evaluations. (a) (b) Fig. 6. (a) CPU solution times per simulation time step and (b) computational overheads in T4 and H8 meshes. Table. 3 Averaged computational overheads in Fig. 6(b). T4 H8 Compared to TherMechTI TherMechExpanTI 18.35% 8.31% TherMechExpanTD 58.72% 37.25% Compared to TherMechExpanTI TherMechExpanTD 34.12% 26.71% The GPU implementation was executed on the same PC using an NVIDIA GeForce GTX 780 GPU with 2304 CUDA cores @ 863 MHz. The GPU computation was performed for the same model sizes and conditions as the CPU computation. Fig. 7(a-b) presents GPU solution times for a single time step in T4 and H8 for TherMechTI, TherMechExpanTI and TherMechExpanTD, and speed improvements (CPU/GPU time ratios). The GPU solution times exhibited a similar trend as the CPU times where they increased almost linearly with the increase of model sizes. Fig. 7(c-d) presents the computational overheads in GPU computation. A summary is presented in Table. 4. (a) (b) (c) (d) Fig. 7. GPU solution times per simulation time step and speed improvements (CPU/GPU) for (a) T4 and (b) H8 meshes, and (c-d) computational overheads. Table. 4 GPU speed improvements over CPU and averaged computational overheads in Fig. 7.
T4 H8 Maximum speed improvements Compared to CPU TherMechTI 17.94x 21.08x TherMechExpanTI 12.46x 16.92x TherMechExpanTD 16.38x 20.97x Computational overheads Compared to TherMechTI TherMechExpanTI 65.86% 31.37% TherMechExpanTD 68.62% 34.52% Compared to TherMechExpanTI TherMechExpanTD 1.67% 2.44%
Fig. 8 presents the ratios of computation time of T4 and H8 in CPU and GPU executions for comparing the computational performance of finite elements. Linear interpolation was used to determine H8 computation times at T4 model sizes due to different mesh densities. At the same model sizes, T4 consumed between 2.3 and 2.8 times more computation time than those of H8 in CPU computation, and between 1.6 and 3.5 times in GPU computation. (a) (b) Fig. 8. Ratios of computation time of T4 over H8 in the (a) CPU and (b) GPU computation.
Bioheat-based temperature prediction and biomechanics-based image registration for surgical simulation of thermal ablation
The presented methodology is applied to simulate and predict soft tissue temperatures and deformations for patient-specific treatment planning of thermal ablation. Fig. 9 illustrates the workflow employed in the present study. Based on the acquired patient-specific medical image dataset, organ geometries were extracted by image segmentation to create 3D surface models for finite element meshing. Soft tissue thermal and mechanical material properties and load and boundary conditions were assigned to generate computational bioheat and biomechanics models for the temperature and displacement fields of soft tissues. The temperature field was used to predict the ablation zones (60 โ isotherms were used due to their strong relationship with the visible boundaries of coagulated tissues based on experimental observations [65]). The displacement field was used to incorporate soft tissue deformations in the medical images by image registration. The 60 โ Table. 5 Simulation settings of the liver model for temperature-dependent anisotropic thermo-visco-hyperelastic analysis.
Geometry Values Unit Description Nodes 3268 13072 DOFs ( ๐ฅ, ๐ฆ, ๐ง, ๐ per node) Elements 18007 Four-node linear tetrahedrons Heat sources ร ๐ท = 0.01 [๐] Heat sources touch each other Material properties Density ๐ [๐๐/๐ ] Specific heat capacity ๐ โ [๐ฝ/(๐๐ โ โ)] Temperature-dependent [45] 4300 @ 90 โ [๐ฝ/(๐๐ โ โ)] Linear interpolation Thermal conductivity ๐ โ [๐/(๐ โ โ)] Temperature-dependent [45] 0.75 @ 90 โ [๐/(๐ โ โ)] Linear interpolation Thermal expansion ๐ผ [1/โ] [70] Transversely isotropic neo-Hookean hyperelastic model ๐น = ๐2 (๐ผ ฬ โ 3) + ๐ ๐ ฬ โ 1) + ๐ 2 (๐ฝ โ 1) Shear modulus ๐ [๐๐] Corresponding to Youngโs modulus 3500 Pa and Poissonโs ratio 0.47 [53] Bulk modulus ๐ [๐๐] Anisotropic coefficient ๐ ๐ ๐ ๐ = 2๐ [๐๐] [64] Unit vector ๐ [1 0 0] To indicate local fibre directions for anisotropic behaviours Prony terms ๐ ๐ ๐ ๐ [๐/๐ ] ๐ ๐ = ๐๐๐ด๐ where ๐๐ด๐ = 1.5 ร 6.104 ร 10 ๐/๐๐ [70] Metabolic heat generation ๐ ๐ [๐/๐ ] [70] Blood perfusion rate ๐ค ๐ [๐๐/(๐ โ ๐ )] [21] Blood specific heat capacity ๐ ๐ [๐ฝ/(๐๐ โ ๐พ)] [21] Arterial blood temperature ๐ ๐ [โ] Exterior of the domain: adiabatic boundary condition Simulation time 25 ๐ Time step 0.0002 ๐ (125000 steps) Fig. 10 presents the simulated temperatures, 60 โ isotherms, and displacement fields in the liver using TherMechTI, TherMechExpanTI and TherMechExpanTD simulations, and Table. 6 presents the results of maximum temperatures and maximum and minimum displacements. The temperature-independent simulation employed the specific heat and conductivity values at 37 โ . The TherMechExpanTD had the lowest maximum temperature, followed by those of TherMechExpanTI and TherMechTI. Due to higher tissue temperatures, TherMechExpanTI had a more pronounced thermal expansion than TherMechExpanTD. In all, TherMechExpanTD had the smallest predicted ablation volumes. Fig. 11 shows the combined results where the predicted ablation zones are displayed in the registered patient-specific Computed Tomography (CT) images considering soft tissue deformations for the operator of ablation to conduct treatment planning and predict treatment outcomes. Temperatures
60 โ isotherms Displacement magnitudes (a) (b) (c) Fig. 10. Simulation results of inducing three heat sources in the liver using (a) TherMechTI, (b) TherMechExpanTI and (c) TherMechExpanTD simulations (note: (a) TherMechTI does not involve thermal expansion; hence, N/A is used in โDisplacement magnitudesโ).
Table. 6 Temperatures and displacements in the liver in Fig. 10 to show the effects of thermal expansion and temperature-in/dependent material properties. ๐ ๐๐๐ฅ (โ) ๐ข ๐ฅ,๐๐๐ฅ (๐) ๐ข ๐ฆ,๐๐๐ฅ (๐) ๐ข ๐ง,๐๐๐ฅ (๐) ๐ข ๐ฅ,๐๐๐ (๐) ๐ข ๐ฆ,๐๐๐ (๐) ๐ข ๐ง,๐๐๐ (๐) TherMechTI 81.8196 0 0 0 0 0 0 TherMechExpanTI 81.7794 1.25e-5 1.65e-5 1.55e-5 -1.19e-5 -1.50e-5 -1.44e5 TherMechExpanTD 78.2727 1.17e-5 1.55e-5 1.46e-5 -1.11e-5 -1.42e-5 -1.36e-5 (a) (b) (c) (d) Fig. 11. The simulated 60 โ isotherms using TherMechExpanTD are visualised on the image registered patient-specific CT images from different views. Discussions
Importance of biomechanics modelling : Biomechanics modelling provides the basis for understanding of many thermo-mechanical behaviours during thermal ablation. In percutaneous RFA, the liver experiences deformations due to needle insertion. In MWA, soft tissues experience expansions close to the antennaโs active tip [11], and an overall tissue shrinkage at the end of the procedure [12]. In other non-invasive thermal ablation treatments such as HIFU ablation, target organs such as the liver and kidneys experience respiration-induced movements. These soft tissue behaviours require the coupled computational bioheat transfer and biomechanics modelling, instead of bioheat transfer alone which was used in many existing works. This will improve the understanding of parametric dependence of thermal doses and exposure outcomes while the target tissues are in dynamics.
Applicability to other hyperthermia treatments : Although the clinically relevant application is demonstrated by a simulation of thermal ablation for local hyperthermia, the principle of the coupled bioheat and biomechanics modelling is generic and can be straightforwardly adapted for regional and whole-body hyperthermia simulations, given the regional and whole-body finite element meshes, material properties, and load and boundary conditions.
Heat source based on modality : In the present work, we used spherical heat source regions with a uniform heat value (similar to Yuan [71]) to simplify the problem for the purpose of algorithm verification. However, the actual heat deposition is usually not uniform and should be determined based on the specific type of modality used for thermal ablation. For instance, it needs to solve the propagation of ultrasound in soft tissues for acoustic pressures for volumetric heat depositions in HIFU [72], Maxwellโs equation for the electromagnetic field in MWA [11] and RFA [73], and laser intensity described by Beer-Lambertโs law in laser ablation [23]. A review of heat source generated by different applications can be found in [74].
Thermal damage of tissues : The extent of thermal damage in tissues can be estimated based on (i) temperature thresholding, (ii) the cumulative equivalent minutes (CEM) at 43 โ , or (iii) an Arrhenius damage integral [46]. In the present work, we used the temperature thresholding based on 60 โ isotherms owing to their strong relationship with the visible boundaries of coagulated tissues based on experimental observations [65]. This is also supported by numerical estimations from the CEM43 and the Arrhenius integral, showing that the tissues will be necrotised almost instantly when the temperature is 60 โ . Anisotropy : The presented method accommodates anisotropic and temperature-dependent properties of soft tissues. Anisotropic thermal conductivity, anisotropic thermal expansion coefficient, and anisotropic hyperelasticity can be incorporated into the constitutive models for soft tissue temperature and deformation computation. Temperature-dependent specific heat capacity and thermal conductivity can also be done similarly.
Computational performance : The computational performance is mainly affected by two factors, the total number of time steps and the number of DOFs. These two factors are controlled by the total simulation time, time step size, and model discretization. A longer simulation time usually yields longer computation time. This can be reduced by using a larger time step size to decrease the total number of time steps; however, the maximum allowable time step size is limited by a critical value in the explicit time integration [75] which is related to the spatial discretization and material properties. With given material properties, a more refined discretization usually leads to a smaller critical time step, which is limited by the smallest element in the mesh, and also leads to an increase in number of DOFs, further increasing computational cost. Therefore, simulations need to consider the mesh discretization for desired efficiency. Some methods such as deep learning [76] and reduce order modelling [77] were reported to permit a larger time step size to be used for more efficient computation.
Finite element considerations : finally, considerations must also be given to the type of finite elements (T4 or H8) used for simulation. As shown in Fig. 8, T4 meshes are between 1.6 and 3.5 times computationally more expensive than H8 with the same number of DOFs; however, the generation of H8 meshes is typically more difficult due to the complicated irregular shape of body organs, requiring significant time and manual interventions to complete a single mesh [69]. In contrast, the key advantage of T4 meshes is that they can be generated automatically based on the surface geometry of patient-specific anatomy, and this is often a standard feature of many mesh processing packages, such as TetGen [68] mentioned previously. Conclusion
Simulation of thermal ablation requires the coupled computational bioheat transfer and biomechanics modelling for accurate prediction of soft tissue temperatures. We achieve this by presenting a thermo-visco-hyperelastic total Lagrangian explicit dynamics finite element algorithm. The essential advantage of the presented method is that it enables full nonlinear modelling of the anisotropic, finite-strain, temperature-dependent, thermal, and viscoelastic behaviours of soft tissues, instead of the linear elastic, linear viscoelastic, and thermal-only modelling in the existing works. To achieve simulations in real-time, the presented algorithm is developed to be well suited for GPU parallel computation. We demonstrate a clinically relevant scenario using a simulation of thermal ablation in the liver. Future works will be devoted to the application of the presented method in simulation of the planning, guidance and training of thermal ablation and other hyperthermia treatments and the development of computer-assisted treatment systems for personalised thermal ablation.
Acknowledgment
This work is funded by the National Health and Medical Research Council (NHMRC), Australia, Grant [1093314].
References [1]
A. Andreozzi, L. Brunese, M. Iasiello, C. Tucci, and G. P. Vanoli, โModeling Heat Transfer in Tumors: A Review of Thermal Therapies,โ
Ann Biomed Eng, vol. 47, no. 3, pp. 676-693, 2019. [2] D. J. Breen, and R. Lencioni, โImage-guided ablation of primary liver and renal tumours,โ
Nat Rev Clin Oncol, vol. 12, no. 3, pp. 175-86, 2015. [3] D. S. Lu, S. S. Raman, P. Limanond, D. Aziz, J. Economou, R. Busuttil, and J. Sayre, โInfluence of large peritumoral vessels on outcome of radiofrequency ablation of liver tumors,โ
J Vasc Interv Radiol, vol. 14, no. 10, pp. 1267-74, 2003. [4] C. Boutros, P. Somasundar, S. Garrean, A. Saied, and N. J. Espat, โMicrowave coagulation therapy for hepatic tumors: review of the literature and critical analysis,โ
Surg Oncol, vol. 19, no. 1, pp. e22-32, 2010. [5] J. Yu, P. Liang, X. Yu, F. Liu, L. Chen, and Y. Wang, โA comparison of microwave ablation and bipolar radiofrequency ablation both with an internally cooled probe: results in ex vivo and in vivo porcine livers,โ
European journal of radiology, vol. 79, no. 1, pp. 124-130, 2011. [6] O. Seror, โAblative therapies: advantages and disadvantages of radiofrequency, cryotherapy, microwave and electroporation methods, or how to choose the right method for an individual patient?,โ
Diagnostic and interventional imaging, vol. 96, no. 6, pp. 617-624, 2015. [7] N. V. Violi, R. Duran, B. Guiu, J.-P. Cercueil, C. Aubรฉ, A. Digklia, I. Pache, P. Deltenre, J.-F. Knebel, and A. Denys, โEfficacy of microwave ablation versus radiofrequency ablation for the treatment of hepatocellular carcinoma in patients with chronic liver disease: a randomised controlled phase 2 trial,โ
The Lancet Gastroenterology & Hepatology, vol. 3, no. 5, pp. 317-325, 2018. [8] L. C. Bertot, M. Sato, R. Tateishi, H. Yoshida, and K. Koike, โMortality and complication rates of percutaneous ablative techniques for the treatment of liver tumors: a systematic review,โ
European radiology, vol. 21, no. 12, pp. 2584-2596, 2011. [9] Y. Fang, W. Chen, X. Liang, D. Li, H. Lou, R. Chen, K. Wang, and H. Pan, โComparison of long-term effectiveness and complications of radiofrequency ablation with hepatectomy for small hepatocellular carcinoma,โ
J Gastroenterol Hepatol, vol. 29, no. 1, pp. 193-200, 2014. [10] M. Moche, H. Busse, J. J. Futterer, C. A. Hinestrosa, D. Seider, P. Brandmaier, M. Kolesnik, S. Jenniskens, R. Blanco Sequeiros, G. Komar, M. Pollari, M. Eibisberger, H. R. Portugaller, P. Voglreiter, R. Flanagan, P. Mariappan, and M. Reinhardt, โClinical evaluation of in silico planning and real-time simulation of hepatic radiofrequency ablation (ClinicIMPPACT Trial),โ
Eur Radiol, vol. 30, no. 2, pp. 934-942, 2020. [11] V. Lopresto, R. Pinto, L. Farina, and M. Cavagnaro, โTreatment planning in microwave thermal ablation: clinical gaps and recent research advances,โ
Int J Hyperthermia, vol. 33, no. 1, pp. 83-100, 2017. [12] L. Farina, Y. Nissenbaum, M. Cavagnaro, and S. N. Goldberg, โTissue shrinkage in microwave thermal ablation: comparison of three commercial devices,โ
International Journal of Hyperthermia, vol. 34, no. 4, pp. 382-391, 2018. [13] A. Andreozzi, M. Iasiello, and P. A. Netti, โA thermoporoelastic model for fluid transport in tumour tissues,โ
J R Soc Interface, vol. 16, no. 154, pp. 20190030, 2019. [14] S. Chung, and K. Vafai, โMechanobiology of low-density lipoprotein transport within an arterial wallโimpact of hyperthermia and coupling effects,โ
Journal of biomechanics, vol. 47, no. 1, pp. 137-147, 2014. [15] M. Iasiello, K. Vafai, A. Andreozzi, and N. Bianco, โLow-density lipoprotein transport through an arterial wall under hyperthermia and hypertension conditionsโAn analytical solution,โ
Journal of biomechanics, vol. 49, no. 2, pp. 193-204, 2016. [16] V. Suomi, J. Jaros, B. Treeby, and R. O. Cleveland, โFull Modeling of High-Intensity Focused Ultrasound and Thermal Heating in the Kidney Using Realistic Patient Models,โ
IEEE Trans Biomed Eng, vol. 65, no. 11, pp. 2660-2670, 2018. [17] K. Sumser, E. Neufeld, R. F. Verhaart, V. Fortunati, G. M. Verduijn, T. Drizdal, T. van Walsum, J. F. Veenland, and M. M. Paulides, โFeasibility and relevance of discrete vasculature modeling in routine hyperthermia treatment planning,โ
Int J Hyperthermia, vol. 36, no. 1, pp. 801-811, 2019. [18] J. Ma, X. Yang, Y. Sun, and J. Yang, โThermal damage in three-dimensional vivo bio-tissues induced by moving heat sources in laser therapy,โ
Sci Rep, vol. 9, no. 1, pp. 10987, 2019. [19] K. Miller, โComputational Biomechanics for Patient-Specific Applications,โ
Ann Biomed Eng, vol. 44, no. 1, pp. 1-2, 2016. [20] Y. Tong, โHigh precision solution for thermo-elastic equations using stable node-based smoothed finite element method,โ
Applied Mathematics and Computation, vol. 336, pp. 272-287, 2018. [21] W. Karaki, Rahul, C. A. Lopez, D. A. Borca-Tasciuc, and S. De, โA continuum thermomechanical model of in vivo electrosurgical heating of hydrated soft biological tissues,โ
Int J Heat Mass Transf, vol. 127, no. Pt A, pp. 961-974, 2018. [22] T. Krรถger, I. Altrogge, T. Preusser, P. L. Pereira, D. Schmidt, A. Weihusen, and H.-O. Peitgen, "Numerical Simulation of Radio Frequency Ablation with State Dependent Material Parameters in Three Space Dimensions,"
Medical Image Computing and Computer-Assisted Intervention โ MICCAI 2006. pp. 380-388. [23] P. Wongchadakul, P. Rattanadecho, and T. Wessapan, โImplementation of a thermomechanical model to simulate laser heating in shrinkage tissue (effects of wavelength, laser irradiation intensity, and irradiation beam area),โ
International Journal of Thermal Sciences, vol. 134, pp. 321-336, 2018. [24] X. Li, Q.-H. Qin, and X. Tian, โThermo-viscoelastic analysis of biological tissue during hyperthermia treatment,โ
Applied Mathematical Modelling, vol. 79, pp. 881-895, 2020. [25] H. H. Pennes, โAnalysis of tissue and arterial blood temperatures in the resting human forearm,โ
J Appl Physiol, vol. 1, no. 2, pp. 93-122, 1948. [26] W. Wulff, โThe energy conservation equation for living tissue,โ
IEEE Trans Biomed Eng , no. 6, pp. 494-495, 1974. [27] H. G. Klinger, โHeat transfer in perfused biological tissue. I. General theory,โ
Bull Math Biol, vol. 36, no. 4, pp. 403-15, 1974. [28] M. M. Chen, and K. R. Holmes, โMicrovascular contributions in tissue heat transfer,โ
Ann N Y Acad Sci, vol. 335, no. 1, pp. 137-50, 1980. [29] A. Nakayama, and F. Kuwahara, โA general bioheat transfer model based on the theory of porous media,โ
International Journal of Heat and Mass Transfer, vol. 51, no. 11-12, pp. 3190-3199, 2008. [30] A. Kotte, G. van Leeuwen, J. de Bree, J. van der Koijk, H. Crezee, and J. Lagendijk, โA description of discrete vessel segments in thermal modelling of tissues,โ
Phys Med Biol, vol. 41, no. 5, pp. 865-84, 1996. [31] D. Y. Tzou,
Macro-to microscale heat transfer: the lagging behavior : John Wiley & Sons, 2014. [32] A. Bhowmik, R. Singh, R. Repaka, and S. C. Mishra, โConventional and newly developed bioheat transport models in vascularized tissues: A review,โ
Journal of Thermal Biology, vol. 38, no. 3, pp. 107-125, 2013. [33] P. Gupta, and A. Srivastava, โNumerical analysis of thermal response of tissues subjected to high intensity focused ultrasound,โ
Int J Hyperthermia, vol. 35, no. 1, pp. 419-434, 2018. [34] J. Zhang, and S. Chauhan, โNeural network methodology for real-time modelling of bio-heat transfer during thermo-therapeutic applications,โ
Artif Intell Med, vol. 101, pp. 101728, 2019. [35] J. Zhang, J. Hills, Y. Zhong, B. Shirinzadeh, J. Smith, and C. Gu, โModeling of soft tissue thermal damage based on GPU acceleration,โ
Comput Assist Surg (Abingdon), vol. 24, no. sup1, pp. 5-12, 2019. [36] M. C. Kolios, A. E. Worthington, M. D. Sherar, and J. W. Hunt, โExperimental evaluation of two simple thermal models using transient temperature analysis,โ
Phys Med Biol, vol. 43, no. 11, pp. 3325, 1998. [37] E. H. Wissler, โPennes' 1948 paper revisited,โ
J Appl Physiol (1985), vol. 85, no. 1, pp. 35-41, 1998. [38] M. A. Solovchuk, M. Thiriet, and T. W. H. Sheu, โComputational study of acoustic streaming and heating during acoustic hemostasis,โ
Applied Thermal Engineering, vol. 124, pp. 1112-1122, 2017. [39] M. Ge, K. Chua, C. Shu, and W. Yang, โAnalytical and numerical study of tissue cryofreezing via the immersed boundary method,โ
International Journal of Heat and Mass Transfer, vol. 83, pp. 1-10, 2015. [40] Y. Zhang, โGeneralized dual-phase lag bioheat equations based on nonequilibrium heat transfer in living biological tissues,โ
International Journal of Heat and Mass Transfer, vol. 52, no. 21-22, pp. 4829-4834, 2009. [41] F. Xu, K. Seffen, and T. Lu, โNon-Fourier analysis of skin biothermomechanics,โ
International Journal of Heat and Mass Transfer, vol. 51, no. 9-10, pp. 2237-2259, 2008. [42] F. Xu, T. Lu, and K. Seffen, โBiothermomechanical behavior of skin tissue,โ
Acta Mechanica Sinica, vol. 24, no. 1, pp. 1-23, 2008. [43] G. C. Bourantas, M. Ghommem, G. C. Kagadis, K. Katsanos, V. C. Loukopoulos, V. N. Burganos, and G. C. Nikiforidis, โReal-time tumor ablation simulation based on the dynamic mode decomposition method,โ
Med Phys, vol. 41, no. 5, pp. 053301, 2014. [44] J. Zhang, J. Hills, Y. Zhong, B. Shirinzadeh, J. Smith, and C. Gu, โTemperature-Dependent Thermomechanical Modeling of Soft Tissue Deformation,โ
Journal of Mechanics in Medicine and Biology, vol. 18, no. 08, pp. 1840021, 2019. [45] S. R. Guntur, K. I. Lee, D. G. Paeng, A. J. Coleman, and M. J. Choi, โTemperature-dependent thermal properties of ex vivo liver undergoing thermal ablation,โ
Ultrasound Med Biol, vol. 39, no. 10, pp. 1771-84, 2013. [46] M. Schwenke, J. Georgii, and T. Preusser, โFast Numerical Simulation of Focused Ultrasound Treatments During Respiratory Motion With Discontinuous Motion Boundaries,โ
IEEE Trans Biomed Eng, vol. 64, no. 7, pp. 1455-1468, 2017. [47] Y.-C. Fung,
Biomechanics: Mechanical Properties of Living Tissues : Springer-Verlag, 1993. [48] J. A. Weiss, B. N. Maker, and S. Govindjee, โFinite element implementation of incompressible, transversely isotropic hyperelasticity,โ
Computer Methods in Applied Mechanics and Engineering, vol. 135, no. 1-2, pp. 107-128, 1996. [49] J. Zhang, Y. Zhong, and C. Gu, โDeformable Models for Surgical Simulation: A Survey,โ
IEEE Rev Biomed Eng, vol. 11, pp. 143-164, 2018. [50] G. A. Holzapfel, โOn Large Strain Viscoelasticity: Continuum Formulation and Finite Element Applications to Elastomeric Structures,โ
International Journal for Numerical Methods in Engineering, vol. 39, no. 22, pp. 3903-3926, 1996. [51] A. G. Holzapfel,
Nonlinear solid mechanics II , 2000. [52] L. Vujosevic, and V. Lubarda, โFinite-strain thermoelasticity based on multiplicative decomposition of deformation gradient,โ
THEORETICAL AND APPLIED MECHANICS, vol. 28, no. 29, pp. 379-399, 2002. [53] J. Zhang, and S. Chauhan, โFast computation of soft tissue thermal response under deformation based on fast explicit dynamics finite element algorithm for surgical simulation,โ
Comput Methods Programs Biomed, vol. 187, pp. 105244, 2020. [54] V. A. Lubarda, โConstitutive theories based on the multiplicative decomposition of deformation gradient: Thermoelasticity, elastoplasticity, and biomechanics,โ
Applied Mechanics Reviews, vol. 57, no. 2, pp. 95-108, 2004. [55] J. Zhang, and S. Chauhan, โFast explicit dynamics finite element algorithm for transient heat transfer,โ
International Journal of Thermal Sciences, vol. 139, pp. 160-175, 2019. [56] J. Zhang, and S. Chauhan, โReal-time computation of bio-heat transfer in the fast explicit dynamics finite element algorithm (FED-FEM) framework,โ
Numerical Heat Transfer, Part B: Fundamentals, vol. 75, no. 4, pp. 217-238, 2019. [57] K. Miller, G. Joldes, D. Lance, and A. Wittek, โTotal Lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation,โ
Communications in Numerical Methods in Engineering, vol. 23, no. 2, pp. 121-134, 2006. [58] K.-J. Bathe,
Finite element procedures , 2006. [59] G. Szekely, C. Brechbuhler, R. Hutter, A. Rhomberg, N. Ironmonger, and P. Schmid, โModelling of soft tissue deformation for laparoscopic surgery simulation,โ
Med Image Anal, vol. 4, no. 1, pp. 57-66, 2000. [60] G. R. Joldes, A. Wittek, and K. Miller, โAn efficient hourglass control implementation for the uniform strain hexahedron using the Total Lagrangian formulation,โ
Communications in Numerical Methods in Engineering, vol. 24, no. 11, pp. 1315-1323, 2007. [61] A. Nealen, M. Mรผller, R. Keiser, E. Boxerman, and M. Carlson, โPhysically Based Deformable Models in Computer Graphics,โ
Computer Graphics Forum, vol. 25, no. 4, pp. 809-836, 2006. [62] S. Cotin, H. Delingette, and N. Ayache, โA hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation,โ
The Visual Computer, vol. 16, no. 8, pp. 437-452, 2000. [63] Z. A. Taylor, M. Cheng, and S. Ourselin, โHigh-speed nonlinear finite element analysis for surgical simulation using graphics processing units,โ
IEEE Trans Med Imaging, vol. 27, no. 5, pp. 650-63, 2008. [64] Z. A. Taylor, O. Comas, M. Cheng, J. Passenger, D. J. Hawkes, D. Atkinson, and S. Ourselin, โOn modelling of anisotropic viscoelasticity for soft tissue simulation: numerical solution and GPU execution,โ
Med Image Anal, vol. 13, no. 2, pp. 234-44, 2009. [65] P. Prakash, โTheoretical modeling for hepatic microwave ablation,โ
The Open Biomedical Engineering Journal,
Magnetic resonance imaging, vol. 30, no. 9, pp. 1323-1341, 2012. [68] H. Si, โTetGen, a Delaunay-Based Quality Tetrahedral Mesh Generator,โ
ACM Transactions on Mathematical Software, vol. 41, no. 2, pp. 1-36, 2015. [69] A. Wittek, N. M. Grosland, G. R. Joldes, V. Magnotta, and K. Miller, โFrom Finite Element Meshes to Clouds of Points: A Review of Methods for Generation of Computational Biomechanics Models for Patient-Specific Applications,โ
Ann Biomed Eng, vol. 44, no. 1, pp. 3-15, 2016. [70] P. Rattanadecho, and P. Keangin, โNumerical study of heat transfer and blood flow in two-layered porous liver tissue during microwave ablation process using single and double slot antenna,โ
International Journal of Heat and Mass Transfer, vol. 58, no. 1-2, pp. 457-470, 2013. [71] P. Yuan, โNumerical analysis of temperature and thermal dose response of biological tissues to thermal non-equilibrium during hyperthermia therapy,โ
Med Eng Phys, vol. 30, no. 2, pp. 135-43, 2008. [72] J. Zhang, N.-D. Bui, W. Cheung, S. K. Roberts, and S. Chauhan, โFast computation of desired thermal dose: Application to focused ultrasound-induced lesion planning,โ
Numerical Heat Transfer, Part A: Applications, vol. 77, no. 6, pp. 666-682, 2020. [73] B. Prasad, J. K. Kim, and S. Kim, โRole of Simulations in the Treatment Planning of Radiofrequency Hyperthermia Therapy in Clinics,โ
J Oncol, vol. 2019, pp. 9685476, 2019. [74] A. Andreozzi, M. Iasiello, and C. Tucci, โAn overview of mathematical models and modulated-heating protocols for thermal ablation,โ 2020. [75] R. Courant, K. Friedrichs, and H. Lewy, โOn the partial difference equations of mathematical physics,โ
IBM journal of Research and Development, vol. 11, no. 2, pp. 215-234, 1967. [76] F. Meister, T. Passerini, V. Mihalef, A. Tuysuzoglu, A. Maier, and T. Mansi, โDeep learning acceleration of Total Lagrangian Explicit Dynamics for soft tissue mechanics,โ
Computer Methods in Applied Mechanics and Engineering, vol. 358, pp. 112628, 2020. [77] Z. A. Taylor, S. Crozier, and S. Ourselin, โA reduced order explicit dynamic finite element algorithm for surgical simulation,โ
IEEE Trans Med Imaging, vol. 30, no. 9, pp. 1713-21, 2011.vol. 30, no. 9, pp. 1713-21, 2011.