Variable thermal conductivity approach for bioheat transfer during thermal ablation

Abstract The aim of the article is to present a mathematical model that predicts tissue temperature during thermal therapy (thermal ablation). The mathematical formulation of variable thermal conductivity of thermal wave bioheat transfer model for a generalized coordinate system due to external heat source during thermal ablation has been given. The governing bioheat transfer equation is simplified by Kirchhoff’s transformation. The numerical solution of the problem has been achieved using a finite difference scheme to discretize in space coordinate; the reduced system of second order ordinary differential equation is solved by the finite element Legendre wavelet Galerkin method (FELWGM). The obtained result from FELWGM is compared with exact analytical solution and shows good agreement. Parametric study is performed to evaluate the effect of variable thermal conductivity on tissue temperature distribution for different physical parameters and illustrated graphically. Particular cases are also deduced from present investigation. The effect of variability of thermal conductivity parameter, variability of time, generalized coordinate system, lagging time and external heat source coefficient on tissue temperature is discussed in detail. Thus, the consideration of the variable thermal conductivity is extremely beneficial for the clinical therapeutic application in the treatment of cancerous cells.


Introduction
For image-guided thermal treatments delivering of heat via focused ultrasound (FUS), radiofrequency (RF), microwave (MW), or laser, temperature induced changes in tissue properties are of relevance in relation to predicting tissue temperature profile, monitoring during treatment, and evaluation of treatment results. These thermal treatments cause reversible or irreversible changes in the properties of tissues. Therefore, for thermal therapies the adequate consideration of tissue properties and their temperature dependence is necessary to achieve accurate results.
Mathematical modelling of thermal therapies has been used extensively to predict and optimize clinical treatments and medical devices. Pennes (1948) model is widely used to model such problems. Kumar, Kumar, and Rai (2015a) studied the heat transfer process in biological tissue during thermal ablation in a generalized coordinate system. Non-linear dual-phase-lag model is utilized by Kumar, Kumar, and Rai (2016) in living tissue during thermal ablation. Kumar and Rai (2017) studied the time fractional DPL model in presence of temperature dependent metabolic and electromagnetic radiation heat source. Kashcooli, Salimpour, and Shirani (2017) investigated the effects of blood vessels on temperature distribution in a skin tissue subjected to various thermal therapy conditions.
Recently, Kumar, Damor, and Shukla, (2018) studied the heat transfer thermal damage in triple layer skin tissue using a fractional bioheat model. For therapeutic treatment of cancerous cells, Kumar, Singh, Sharma, and Rai (2018) used non-linear dualphase-lag model under Dirichlet boundary conditions. Tan, Zou, Dong, Ding, and Zhao (2018) investigated the coupled effects of blood vessels and relaxation time on tissue temperature and thermal lesion region during HIFU hyperthermia.
Very few authors considered the temperaturedependent thermal parameters over wide temperature ranges in tissues. Valvano, Cochran, and Diller (1985) measured thermal conductivity and diffusivity of biomaterials over a wide temperature range. Bhattacharya and Mahajan (2003) presented the temperature-dependent thermal conductivity of sheep collagen and cow liver in vitro. Watanabe, Kobayashi, Hashiszume, and Fujie (2009) investigated the effects of the temperature dependence of thermophysical properties of the liver on the temperature distribution during RFA. Trujillo and Berjano (2013) reviewed the mathematical functions most commonly used for modelling the temperature dependence of electrical and thermal conductivities of biological tissue in RFA. Guntur, Lee, Paeng, Coleman, and Choi (2013) found that the thermal properties of ex vivo porcine liver tissue during RFA vary significantly with temperature. Rossmann and Haemmerich (2014) reviewed temperature-dependent thermal and electrical properties of biological tissues in the supraphysiological temperature regime. Guntur and Choi (2015) examined the effects of temperature-dependent thermal parameters on the temperature distribution in liver tissue exposed to a clinical high-intensity focused ultrasound (HIFU) field, and found that alteration of tissue thermal parameters significantly affects the prediction of temperature distributions and suggest that the dependence of thermal parameters should be considered in estimation of the thermal dose in tissues exposed to a clinical HIFU field.
The present study aims to investigate the effects of variable thermal conductivity during thermal ablation for a thermal wave bioheat transfer model with modified Gaussian external heat source in a generalized coordinate system. The governing bioheat transfer equation is simplified by Kirchhoff's transformation and then finite difference scheme, and Legendre wavelet Galerkin approach is used to solve the problem. The problem is converted into a system of second order differential equation with initial conditions by using the finite difference scheme in space coordinate. Then, this reduced system is solved by the Wavelet Galerkin approach with Legendre wavelet as basis function, which reduces the problem into the Sylvester matrix equation. The solution of the Sylvester matrix equation gives the temperature distribution during thermal ablation with variable thermal conductivity and a significant effect of the variable conductivity has been observed on tissue temperature distribution.

Formulation of the problem
Consider a tissue of finite length (L), which is heated by modified Gaussian external heat source with variable thermal conductivity during thermal ablation in generalized coordinate system. Pennes (1948) performed a series of experiments that measured temperatures on human forearms and derived a thermal energy conservation equation known as bioheat transfer equation as follows: From the principle of conservation of energy to tissue control volume (V), we have where U gain , U storage , and U loss are given below.

Nomenclature
C the number of classify coordinate j tissue thermal diffusivity, m 2 s À1 x b blood perfusion rate, s À1 q tissue density, kgm À3 q b blood mass density, kgm À3 s q heat flux relaxation time, s H dimensionless tissue temperature, in Equation (28) h dimensionless tissue temperature # a dimensionless arterial blood temperature # w dimensionless wall temperature at boundary a 0 scattering parameter, m À1 c specific heat of tissue, jkg À1 c b specific heat of blood, jkg À1 k À1 F 0q dimensionless heat flux relaxation time F 0 dimensionless time k variable thermal conductivity of tissue, Wm À1 K À1 k 0 constant thermal conductivity of tissue, Wm À1 K À1 K 1 dimensionless variable thermal conductivity parameter k 1 variable thermal conductivity parameter, C À1 L tissue slab length, m l Bromwich contour integration line P transmitted power, W P f dimensionless blood perfusion coefficient P m0 dimensionless metabolic heat source coefficient P r0 dimensionless external heat source coefficient q b heat source due to blood perfusion, Wm À3 q m heat source due to metabolic heat generation, Wm À3 q s external heat source, Wm À3 q m0 constant metabolic heat source q r0 reference external heat source, Wm À3 r space coordinate, m r 0 location of tumor, m S antenna power, m À1 s Laplace domain parameter T tissue temperature, C t time, s T 0 initial temperature of body, C T a arterial blood temperature, C T w wall temperature at the boundary, C x dimensionless space coordinate x Ã dimensionless location of tumor where heat is applied The general internal heat generation per unit volume due to metabolic heat (q m ) and the interior heat caused by outer heat source (q s ) is The total rate of stored thermal energy over the control volume is Heat loss to adjacent tissues by convection and conduction is given as where is the rate production through the control volume and the total amount of thermal energy over the tissue control volume. Using Equations (2)-(6) in Equation (1) Applying the divergence theorem to the surface integral and observing that Equation (7) must hold for any arbitrary volume element, we obtain qc @T @t Àq b Àq m Àq s ¼ À @q @r : Cattaneo (1958), Vernotte (1958) proposed the following modified heat flux model that is a linear extension of the Fourier's law in order to take in the wave-like behaviour in heat conduction process, where k(T) is the thermal conductivity with variable temperature. Using Equation (9) in Equation (8) with temperature dependent specific heat yields Equation (10) in generalized coordinate system can be written as where C ¼ 0; 1; 2 corresponds to Cartesian, axis-symmetric, and spherical symmetric coordinate system and The term q b can be expressed as The term q m can be approximated as a linear function of local tissue temperature taken in Gupta, Singh, and Rai (2010), Following (Kumar, Kumar, & Rai, 2015b) modified Gaussian heat source q s is taken as To solve the problem, the following initial and boundary conditions are taken: Initial conditions Symmetric condition @T L; t ð Þ @r ¼ 0: Considering the mapping (Kirchhoff's transformation), following Hetnarski (1986) as An important factor to achieve realistic models is the use of mathematical functions to describe the temperature dependence of thermal properties of tissue. A commonly used approach for modelling the temperature dependence of thermal properties for temperature below 100 C is based on linear equations and employs constant temperature coefficient such as: From the above, we have the following relations: and Equations (11) and (16-18) with the aid of Equations (19-21) in linear form become @# @r 0; t ð Þ ¼ 0:

Solution of the problem
Introducing the non-dimensional variable and similarity criteria With the help of Equation (27), Equations (24-26) reduce to the following form

Finite difference scheme for spatialdiscretization
The space coordinate domain [0,1] divided into n þ 1 sub intervals of equal step length (h) i.e. 0 ¼ x 0 <x 1 <x 2 < . . . <x n <x nþ1 ¼ 1, where x iþ1 ¼ x i þ h. The first and second order derivatives are approximated by using central difference formula i.e.
1 i n: Using this spatial discretization scheme, Equation (28) becomes 1 i n; F 0 >0, with initial condition boundary condition The Equation (34) is converted into a system of second order ordinary differential equation in vector matrix form as follows: subjected to initial condition where n ¼ 2 kÀ1 M, I is the identity matrix of order n, A is a square matrix of order 2 kÀ1 M Â 2 kÀ1 M; B, Hð0Þ, and dHð0Þ dF 0 are the column matrix of 2 kÀ1 M Â 1 which are given as follows
Using Equations (39), (42), and (46) in Equation (37), we have f ÀF 0q P m0 dÞPÀðP m0 dÀP 2 f ÞP 2 I n is the identity matrix of order n and d 1 is n Â 1 matrix such that d T 1 wðF 0 Þ ¼ 1. Equation (47) is Sylvester matrix equation, after solving this equation, we obtain the element of matrix C T . Substituting the value of C T in Equation (46) gives tissue temperature distribution HðF 0 Þ. The tissue temperature in terms of T can be obtained from Equation (23).
Particular cases: i. If we take K 1 ¼ 0 in Equation (47), the problem is reduced to the non-Fourier bioheat transfer model during thermal ablation and the corresponding results are the same as obtained by Kumar, Kumar, and Rai (2015a). ii. If we take s q ¼ 0 in Equation (47), then we obtain results for the Pennes' bioheat transfer model with variable thermal conductivity during thermal ablation. iii. If we take s q ¼ 0 and K 1 ¼ 0 in Equation (47), then we obtain results for the Pennes' bioheat transfer model in the absence of variable thermal conductivity during thermal ablation. Yadav, Kumar, and Rai (2014) showed that the Legendre wavelet Galerkin method produces a higher degree of accuracy and minimizes error. That is why we used this approach for numerical procedure. This approach combines both the multi-resolution and multi-scale computational property of Legendre wavelet with element-wise analysis. Multi-resolution analysis of Legendre wavelet localizes small scale variation of solution and fast switching of functional bases. A brief algorithm of the numerical method is given in Table 1. The problem is solved analytically by Laplace transform in the absence of modified Gaussian external heat source, which is given as follows:
q ¼ 1000 kgm À3 ; C ¼ 4:18 Â 10 3 Jkg À1 K À1 ; ; L ¼ 0:05 m; q r0 ¼ 7:58 Â 10 3 Wm À3 ; q m0 ¼ 1:091 Â 10 3 Wm À3 ; W b ¼ 0:5 kgm À1 s À1 ; k 0 ¼ 0:5 Wm À1 K; a 0 ¼ 1 Â 10 2 m À1 ; r 0 ¼ 0:025 m; s q ¼ 20 s; t ¼ 6 min: MATLAB(R2016a) software has been used for all computational work. Thermal ablation refers to the destruction of tissue by extreme hyperthermia (elevated tissue temperatures) or hypothermia (depressed tissue temperatures). The temperature change is concentrated to a focal zone in and around the tumour. The high temperature occurs at the target region for four to six minutes during the process. Figure 1 shows the effect of variable thermal conductivity parameter K 1 on tissue temperature distribution. It is noted that as the value of K 1 increases, the temperature increases and the temperature attains maxima at the heating position x Ã ¼ 0:5. Thus, to achieve more accuracy in clinical treatments variable thermal conductivity plays an important role. Figure 2 shows the effect of time on tissue temperature along the space coordinate when K 1 ¼ 0:1. It is evident that temperature increases as the time increases from four minutes to six minutes during thermal ablation. Thus, it is ned that time duration affects the treatment process. Figure 3 shows the effect of a generalized coordinate system on the treatment of tumour during thermal ablation for K 1 ¼ 0 and K 1 ¼ 0:1. It is noted that for all three coordinate systems, tissue temperature remains unaffected at the heating position x Ã ¼ 0:5 for both the cases, but for x>x Ã and x<x Ã , it affects the tissue temperature. Figure 4 shows the comparison of temperature distribution for Pennes (s q ¼ 0) and thermal wave model ðs q ¼ 20sÞ. It is observed that temperature profile for the thermal wave model is less than the temperature profile predicted by the Pennes' model for both the cases of K 1 ¼ 0 and K 1 ¼ 0:1. But, independently, both models predict higher temperature with the variable thermal conductivity parameter. Thus, lagging time is affecting the temperature distribution during the treatment. Figure 5 shows the effects of external heat source on temperature distribution for K 1 ¼ 0 and K 1 ¼ 0:1. It is clear that as the value of K 1 and q r0 increases, temperature also increases. Thus, the consideration of external heat source is important during the treatment. Figure 7 depicts the absolute error of the solutions obtained by finite element Legendre wavelet  Galerkin method (FELWGM) and analytical method. Temperature (T) for time (t) 0 to 500 s is evaluated by both the methods and results are given in Table 2. Therefore, we can say that the consideration of variable thermal conductivity in bioheat transfer model is very beneficial for the clinical therapeutic application in the treatment of cancerous cells.

Conclusions
The problem of modelling the variable thermal conductivity of thermal wave bioheat transfer model in a generalized coordinate system due to external heat source during thermal ablation has been investigated. Kirchhoff's transformation is used to formulate the governing bioheat transfer equation. The numerical solution of the problem is achieved by using a finite difference scheme and Legendre wavelet Galerkin approach. FELWGM converts the problem into a system of algebraic equation, so it reduces the computational complexity and avoids large storage requirements, as in the case of conventional discretize schemes. The results obtained from the FELWGM are compared with the analytical results and give good accuracy in the absence of external heat source. The significant effect of variable thermal conductivity parameters, variability of time, lagging times, and external heat source is observed on tissue temperature distribution. Therefore, the temperature dependent changes need to be considered to enhance the accuracy of predicting and monitoring pre-clinical and clinical treatments and consideration of this model is useful in image-guided thermal treatments (e.g. hyperthermia and thermal ablation) and other biomedical treatments. It will give new direction in clinical therapeutic application.