Gradient of the temperature function at the voxel (i,j,k) for heterogeneous bio-thermal model

Determination of the relationship between electromagnetic power absorption and temperature distributions inside highly heterogeneous biological samples based on numerical methods is essential in biomedical engineering (e.g. microwave thermal ablation in clinic). In this paper, the gradient expression is examined and analyzed in detail, as how the gradient operators can be discretized is the only real difficulty to the solution of bio-heat equation for highly inhomogeneous model utilizing implicit scheme.Determination of the relationship between electromagnetic power absorption and temperature distributions inside highly heterogeneous biological samples based on numerical methods is essential in biomedical engineering (e.g. microwave thermal ablation in clinic). In this paper, the gradient expression is examined and analyzed in detail, as how the gradient operators can be discretized is the only real difficulty to the solution of bio-heat equation for highly inhomogeneous model utilizing implicit scheme.


I. INTRODUCTION
Determination of thermal response due to electromagnetic (EM) exposure based on numerical methods is a difficult interdisciplinary research involving mathematics, computer science, physics and other disciplines. The bio-heat transfer equation for homogeneous material model can be easily calculated by using second order finite difference approximation to discretize the spatial derivatives and explicit scheme for time domain discretization. Mr. Gandhi and colleagues solved the bioheat equation for inhomogeneous models with different tissues properties utilizing implicit finitedifference method. [1][2][3][4][5] Whereas we appreciate their research, we would like to address a few issues that may help further clarify or confirm the research.

II. DETAILED ANALYSIS OF THE GRADIENT EXPRESSION IN BIO-HEAT EQUATION
In the Refs. 1-3, instead of showing a discretization of the differential equation by their methods "an implicit finite-difference method that achieves a higher order accuracy of the Crank-Nicolson formulation", however, the following equation (1) was given by solely subscripting all variables of the partial differential equation with voxel indices i, j, k.
The derivation is believed to be correct but is odds with previous published work. First, the expression ∇T i,j,k in (1) is unsuitable since from Taylor's theorem, 6,7 T i,j,k is not a function of space but a numerical approximation of the function "T " at the voxel (i, j, k). It will be elaborated below.
In the three-dimensional (3D) Cartesian coordinate system, the gradient of scalar temperature function T is given by ∇T = ∂T ∂x , ∂T ∂y , ∂T ∂z . For a voxel (i, j, k): The approximation for partial derivatives of temperature function e.g. (2) is obtained from Taylor's theorem: (5) is a numerical approximation of temperature function T at the point (i∆x, j∆y, k∆z). It is a fixed value.

• • •
Following the reference work in mathematics, 7 the gradient of the temperature function T at the voxel (i, j, k) can be expressed as ∇T i,j,k or (∇T ) i,j,k , Similarly, the correct notations for the gradient of T in (10) 4,5 can be expressed as ∇T | i or (∇T ) i The paper 10 they cited presented a method that is for isotropic medium of constant thermal conductivity k, not for highly heterogeneous model with thermal conductivity varying from tissue to tissue as shown in Gandhi and colleagues' articles. [1][2][3][4][5] Readers may expect to see how the gradient operators ∇(k i,j,k ∇T i,j,k ) in (1) and ∇(k i ∇T i ) in (10) can be discretized for their inhomogeneous human head or human leg models, as that is the only difficulty to the solution of bio-heat equation utilizing implicit scheme. Therefore, it would be interesting if Mr. Gandhi and colleagues were to show some discretization by their methods ("The heat conduction equation is solved in rectangular coordinates by a method which is a modification of Douglas and Rachford that achieves higher-order accuracy of the Grank-Nicholson formulation. The unconditionally stable system similar to that of Douglas and Rachford method had the advantage of using larger time steps, while the unknown temperatures are computed explicitly at the advanced time level, similar to that in the Grank-Nicholson formulation having higher-order accuracy", 4 "The heat conduction equation is solved in rectangular coordinates by the standard implicit finite difference technique which is stable for all size time steps" 5 and "an implicit finite-difference method that achieves a higher order accuracy of the Crank-Nicholson formulation" 1,2 ) since these methods are neither common finite-difference schemes nor in the references they cited.
People can get the similar results using different methods in many cases, e.g. if someone does not use "Implicit Finite-Difference method that achieves a higher order accuracy of the Crank-Nicholson Formulation" as the authors mentioned in their articles, however, using explicit methods can also obtain the similar results which only take much longer computational time. With that being said, we don't think it is necessary to discuss the implementation and/or the results reported in their articles. As ∇T i,j,k equals zero, by expanding it in its finite-difference approximation, 8,9 the flawed equation (1) is rewritten as where n is the time step index and δ t is the time step size.
In (11), the new value of a temperature at any voxel depends only on its previous value and has nothing to do with the temperature values at adjacent voxels. Obviously, by solely subscripting all variables of the partial differential equation with voxel indices i, j, k, the bio-heat conduction equation has become a non-conductive false equation.

III. NUMERICAL METHODS TO THE SOLUTIONS OF TIME-DEPENDENT DIFFERENTIAL EQUATIONS
Numerical approximations to the solutions of time-dependent differential equations consist of the explicit, implicit and semi-implicit schemes. A common approach is the explicit finite-difference method due to its robustness and ability to discretize strongly inhomogeneous structures in a straightforward manner. However, the explicit scheme requires a long computation time 11 due to the stability criterion, which limits its practical applications especially for in-vitro studies on small cells. To overcome the disadvantages of the explicit method, the unconditionally stable implicit and semi-implicit schemes have been investigated by researchers. [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] We propose numerical methods that apply the fractional steps method (FSM) 32 to address the bio-heat equation (or other partial differential equations) in 3-dimension as a result of its simple implementation and computational efficiency. The proposed algorithm can be easily implemented within 20 lines of code for heat equation.
The temperature distribution T = T ( r, t) inside the electromagnetically exposed multi-material bio-model can be obtained by using the bio-heat transfer equation as follows: 33 where K is the thermal conductivity [J/(s·m· • C)], A 0 is volumetric heat sources [J/(s·m 3 )] due to metabolic processes, the specific absorption rate (SAR) is the input EM heating source into the bio-heat equation, B is a parameter proportional to the blood perfusion [J/(s·m 3 · • C)], T b is the blood temperature, and C p [J/(kg· • C)] and ρ [kg/m 3 ] are the tissue specific heat and density, respectively.
Using the notations in field theory, bio-heat transfer equation (12) is equivalent to: The semi-implicit Crank-Nicolson formula 12,13 for bio-heat transfer equation (13) with variable coefficient values can be expressed as: where T is temperature, ∆t is the time step size, ∆x, ∆y and ∆z are the spatial meshing sizes along x, y and z axes respectively, i, j and k are the integer-number indices of the computational cells and n is the time step index. It is clear that to solve the above equation (14) in direct way is almost impossible. In the proposed method, the FSM method is used to address the Crank-Nicholson scheme in 3D to overcome the drawback of difficult implementation due to the huge equations in direct Crank-Nicholson calculation. Eq. (13) is solved in each direction as follows (15), (16), (17).
for ∀j, then for all i, k for ∀k, then for all i, j where ∆x i = x i+1 − x i ; ∆y j = y j+1 − y j ; ∆z k = z k +1 − z k . Here, (n + 1 3 ), (n + 2 3 ), and (n + 1) correspond to the time steps of t n+ 1 3 = t n + 1 3 ∆t, t n+ 2 3 = t n + 2 3 ∆t, and t n+1 = t n + ∆t, respectively.  The equations (18) can be simplified to yield an implicit (tridiagonal) updating scheme as follows and easily calculated by using the chase method. Thomas, J. W. shows that Crank-Nicolson scheme is unconditionally stable when applied to u t = au xx . 13 For highly inhomogeneous model, however, the coefficients or spatial meshing sizes along x, y and z axes are different from each other, the above matrix system would not be skew-symmetric and thus the results might get unstable. The temperature updating equation (18) is therefore reformed as Eq. (22) which is unconditionally stable for inhomogeneous models utilizing non-uniform meshes and with multi-tissue.
In the same manner, we can get equations in y and z directions. In this way, three one-dimensional equations are calculated instead of one three-dimensional equation. The solution to the next timelayer is obtained by solving three one-dimension formulas which can be utilized for calculating inhomogeneous model with non-uniform space-increments and various tissues very quickly and accurately. Detailed calculation can be partly found in our previous paper. 31 An in-vitro study on pineal gland at 1800 MHz shown in Figure 1 in a sector of the radial waveguide exposure system 34 was conducted. TABLE I lists the thermal properties of all materials. Figure 2 shows the SAR distributions inside this model. The temperature rise of the target model at 0.8 W/kg average SAR in the isolated pineal gland sample after 5 min exposure was calculated by individually using our proposed solution and the finite volume method 35,36 (FVM) which is fully compatible with the finite-difference method and allows the use of highly heterogeneous objects. Figure 3-4 show the computed temperature distributions from both methods, respectively. The observation plane is at the centre of the in-vitro model. When using the same spatial discretization of smallest space-increment 0.2 mm, the computation time of the explicit method with time step of 0.000375 s is 23 hours while similar results can be achieved in 10 minutes by using our proposed solution with time step of 0.1 s.

IV. CONCLUSION
The bio-heat equation in Refs. 1-5 is examined. The derivation is believed to be correct but is odds with previous published work. We propose numerical methods that apply the fractional steps method to address the heat equation in 3D as a result of its simple implementation and computational efficiency. A numerical example that utilize the FSM method to quickly and accurately solving 3D bio-heat equation for highly heterogeneous pineal gland model is shown.