An efficient parallel numerical modeling of bioheat transfer in realistic tissue structure

International Journal of Heat and Mass Transfer 95 (2016) 843–852

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

An efficient parallel numerical modeling of bioheat transfer in realistic tissue structure Zhi-Zhu He a,⇑, Jing Liu a,b a Key Laboratory of Cryogenics, and Beijing Key Laboratory of Cryo-Biomedical Engineering, Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, Beijing 100190, China b Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China

a r t i c l e

i n f o

Article history: Received 3 July 2015 Received in revised form 15 November 2015 Accepted 14 December 2015

a b s t r a c t A fast and accurate numerical method of solving bioheat transfer problems is critical for many biomedical applications, such as efficient implementation of thermal therapy planning in a clinical setup, and evaluation of the thermal effects induced by specific energy absorption rate (SAR) from external electromagnetic field. This paper developed a parallel finite-difference method based on alternating direction explicit (ADE) scheme to solve three dimensional transient bioheat equation for the realistic tissue structure. An optimized processing pipeline for accurate presentation of the irregular boundary condition was established to considerably reduce the staircase errors induced by cubic voxel. The voxels order aligned along diagonal direction was designed to allow fast parallel strategy for the lower and upper matrix solution involved in ADE scheme for irregular tissue. The test cases including a real head model have indicated that the developed method can support large time step and cubic voxel approximation of irregular interface and boundary, and also produce significant parallel efficiency. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The heat transfer in biological tissue is one of the most fundamental physical principles, which significantly affects the biochemical reaction and physiological process. The prediction of the temperature evolution in living tissue is of great importance for modeling of many clinical practices including hypothermia therapy for cancer treatments [1], brain cooling [2,3], or evaluation of the thermal injury caused by external device such as the electromagnetic field of magnetic resonance imaging [4] and cellular phones [5]. Variants of the Pennes bioheat transfer model have been successfully applied to estimate the tissue temperature evolution, which have been validated by experiments [6]. In a simulation framework, however, one of the central questions is to develop a fast and accurate numerical method for solving bioheat transfer equation. Recently, the simulation approaches for bioheat transfer have been reviewed in [7]. Among the existing simulation methods, the finite difference method (FDM), or denoted by finitedifference time-domain (FDTD), is widely applied for biomedical applications, especially in the context of increasing concern over ⇑ Corresponding author at: Technical Institute of Physics and Chemistry, Chinese Academy of Sciences, P.O. Box 2711, Beijing 100190, China. Tel./fax: +86 10 82543766. E-mail address: [email protected] (Z.-Z. He). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.12.028 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

electromagnetic radiation exposure [4]. One of the reasons is that FDM is easily implemented in terms of coding, where no complex mesh generation process involved in finite element method (FEM) is demanded. In addition, cubic voxel used in FDM can be directly obtained from the medical imaging such as the magnetic resonance imaging (MRI) and the computed tomography (CT), which are also represented by cubic voxel. However, FDM suffers considerably from staircase errors induced by curvilinear boundaries [8], and time-step limitation in explicit scheme [7], which have significant impact on the computational accuracy and efficiency in particular. Samaras et al. [9] proposed a conformal scheme to reduce the staircase effect, which modifies the boundary area to improve the calculation accuracy of heat flux at irregular boundary. The alternating-direction implicit (ADI) approach allow a larger time step to rapidly compute temperature evolution in biological tissue [10]. Our previous work [7] has developed a boundary-source method to efficiently correct the heat transfer of boundary, which was further combined with ADI method to improve the computational efficiency. The immersed boundary method [6] has also been applied to solve complex heat transfer of cryo-freezing tissue. In order to rapidly evaluate temperature increase induced by electromagnetic field, the spatial filter method [11] based on fast Fourier transform has been developed to significantly reduce computational time. However, it cannot capture the accurate thermal effects of irregular boundary condition, which is often important for many applications, such as scalp cooling for thermal regulation

844

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

of brain injury [2]. In addition to developing efficient numerical scheme, parallel method could also speed up computational process. Zeng et al. [12] has developed an efficient parallel ADI method for application in cryoablation. However, the parallel strategy has not been widely used for bioheat transfer simulation because of its complexity, respectively, considering complex structure. This paper is aimed to develop a simple and efficient parallel alternating direction explicit (ADE) scheme to solve three dimensional transient bioheat equation for the real tissue structure. ADE method has been used for solving heat conduction equation only for the simple cubic geometry [13,14], while it has not been applied for bioheat transfer solution with irregular tissue structure. More importantly, we here developed an efficient parallel strategy for ADE with complex tissue structure based on special voxel alignment and rearrangement method. In addition, the optimized processing pipeline for the irregular boundary condition was established, which could be directly combined with both medical image reconstruction and the developed ADE scheme. The paper is organized as follows: firstly we describe a modified bioheat transfer equation and its ADE discretization considering irregular boundary. Afterwards the parallel method for ADE is presented. The detailed numerical simulation tests are subsequently implemented to validate the numerical scheme developed in current work. 2. Methodology 2.1. Modified bioheat equation considering irregular boundary condition In our previous work [7], we have constructed a boundarysource Qw = Tw–AwT to correct the heat transfer from irregular boundary. Thus a modified Pennes’ type bioheat equation considering irregular boundary condition is given by He et al. [7]

@T qc ¼ r  ½jrT þ xb C b ðT b  TÞ þ Q m þ Q e þ ðT w  AW TÞ @t

ð1Þ

where T is the temperature and t denotes the time, q is the tissue density, c is the specific heat capacity, k is the thermal conductivity, xb is blood perfusion from capillary blood, subscript b denotes blood, Qm is metabolic heat generation, and Qe is exterior heat generation such as SAR induced by electromagnetic field exposure. The detailed parameters (Tw and Aw) description involved in boundary heat source will be discussed at next section. The current modified bioheat transfer model requires the adiabatic condition on the tissue boundary, which can be automatically met by enforcing zero thermal conductivity to voxels outside tissue domains. This method allows cubic voxel approximation of the complex tissue structure without explicitly representing the irregular boundary. 2.2. Optimized processing pipeline for irregular boundary condition For

general

thermal

boundary

condition

expression

a@T=@n þ bT ¼ c, the parameters Tw and Aw involved in boundary source term can be estimated by He et al. [7]

AW ¼

Sjw b ; Vðbl þ aÞ

Tw ¼

Sjw c Vðbl þ aÞ

ð2Þ

where parameters a, b and c are determined by the detailed boundary conditions, such as a = 1 and b = c = 0 denoting adiabatic condition. The geometric parameters S, V and l denote, respectively, the intersection surface area, volume and distance for the boundary voxel, which are completely determined by the boundary geometry information. Here, we give an optimized processing pipeline for the

above three parameters determination, which has four main steps summarized in Fig. 1. 2.2.1. Step 1: representing the boundary with smooth index The geometric structure of tissue can be obtained from segmentation and reconstruction of medical images such as MRI. Image segmentation provides volumetric quantification of the tissue shape. The binary mask (denoted by index / = 1 and / = 1) for each tissue can be determined based on labeled voxels according to tissue type (see Fig. 1A). We can smooth the binary value defined at the voxel vertex to represent the tissue boundary with zero level value [15] (see Fig. 1B). It is noteworthy that the level set based segmentation method can be directly used for boundary representation with zero level set [16]. 2.2.2. Step 2: voxel classification The voxels can be classified as three types: tissue voxel defined by / 6 0, background voxel for / > 0, and boundary voxel for other condition (illustrated in Fig. 1C). 2.2.3. Step 3: geometric parameters calculation For each boundary voxel, the intersection surface can be estimated by isosurface of / = 0. Firstly, the intersection point (such as x1 and x2 in Fig. 1D or 3D case in Fig. A1 of Appendix A) determined by / = 0 at the edge of boundary voxel is calculated, and P the point on the intersection plane xc ¼ M1 M i¼1 xi is estimated, where M is the number of intersection points for all the edges of the voxel. n ¼ r/=jr/j denotes the normal direction of the plane. Then we can obtain the plane equation by n  (X–Xc) = 0. Furthermore we could estimate the above geometric parameters based on the analytical relations connecting with plane equation. We should enforce l = 0 for l < 0. The detailed calculation of S and V are given in Appendix A. 2.2.4. Step 4: boundary voxel rearrangement Not all the boundary voxels directly participate in heat transfer calculation. The boundary voxel with small volume (here defined as the smaller than 20% voxel volume) should be discarded to ensure numerical scheme stability. In order to correcting the heat transfer from boundary condition, however, we rearrange the discarded boundary voxel to the adjacent tissue or boundary voxel with large volume (larger than 20% voxel volume), which are both called as active voxel (illustrated in Fig. 1E). For each discarded boundary voxel, the occupied S and V are redistributed to the adjacent active voxel according to the normal direction of the intersection plane. The parameter S would affect the heat flux calculation, and V has impact on the heat generation (such as metabolic heat generation), which would be scaled by the increment volume ratio. A sample of the rearranged increment calculation is given in Fig. 1E. For voxel with multiple tissue surfaces, we could divide the voxel into many sub-voxels (such as 27 smaller sub-voxels) to calculate the surface and volume of each tissue segment. The geometric parameters S, V and l need to be calculated only once at the boundary voxels. Thus the computational time can be ignored. 2.3. Discrete scheme of diffusion term The second-order partial derivative of T with respect to x adopts a central-differencing scheme:

@ @x

    n jiþ1=2;j;k T niþ1;j;k  T ni;j;k  ji1=2;j;k T ni;j;k  T ni1;j;k @T  j ¼ 2 @x i;j;k Dh



ð3Þ

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

845

Fig. 1. The schematic diagram of optimized processing pipeline for boundary condition: (A) the binary mask for each tissue determined based on labeled voxels according to tissue type; (B) the tissue structure represented by the smooth index; (C) the voxels classified as three type: tissue voxel, background voxel and boundary voxel; (D) the intersection surface estimated by isosurface of / = 0 for each boundary voxel; (E) the active boundary rearrangement.

where index (i, j, k) denote the center of voxel, and Dh is the voxel width. The thermal conductivity jiþ1=2; j;k is defined at the voxel boundary (i + 1/2, j, k). For the voxel without tissue interface intersection, jiþ1=2;j;k can be considered as ji;j;k . However, the weighted harmonic interpolation method is used at voxel boundary with tissue interface intersection

jiþ1=2;j;k ¼

ji;j;k jiþ1;j;k

hx jiþ1;j;k þ ð1  hx Þji;j;k

ð4Þ

The weighted parameter is given by hx ¼ DxI =Dh, where DxI denotes the interface distance from voxel center (see Fig. 2A). It can be estimated through / interpolation

hx ¼

j/i;j;k j j/i;j;k j þ j/iþ1;j;k j

/i;j;k /iþ1;j;k < 0

ð5Þ

It is noteworthy that / here has been smoothed by Step 1 in Section 2.2. For the small difference between thermal conductivities, approximation jiþ1=2;j;k ’ ji;j;k cannot induce the large calculation

Fig. 2. The schematic diagram of the weighted parameter determination (A) and (B) the voxels aligned along diagonal order.

846

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

error [8]. However, it is more accurate to adopt Eq. (4) for the larger thermal conductivities difference. For voxel with multiple tissue interfaces, we firstly calculated the occupied volume ratio for each tissue. The interpolation parameter hx in Eq. (4) is replaced by volume ratio of the tissue. In addition, the heat generation for the interface voxel adopts linear interpolation based on the volume ratio. The second-order partial derivative of T with respect to y and z also adopts the above discretization scheme. 2.4. ADE scheme for bioheat transfer equation

ð6Þ

where a ¼ ðxb qb cb þ Aw Þ=qc and q ¼ ðxb qb cb T b þ Q m þ Q e þ T W Þ=qc, T denotes the temperature set (Ti,j,k for all indexes of i, j and k) from all the voxels, and I is the unit vector with value one. A is symmetric seven-diagonal matrix with responding to three dimensions diffu^ þ D þ U, ^ where D is sion term in Eq. (1). A is decomposed as A ¼ L ^ and U ^ are the strictly lower and upper trithe diagonal part of A, and L ^ þ D=2 and angular part of A, respectively. If we define L ¼ L ^ þ D=2, we could express ADE scheme [13] of Eq. (6) as U¼U

8  

1  

a ^ > 1  2a Dt T n I þ DtT n U þ DtqI > > T ¼ 1 þ 2 Dt I  DtL <  

 

 ¼ 1 þ a Dt I  DtU 1 1  a Dt T n I þ DtT n L þ DtqI T 2 2 >   > > : T nþ1 ¼ 1 T ^þT  2

All the simulations are performed at Intel Xeon E5-2630 station with twelve-core processors and shared memory (2.3 GHz and 64G RAM). For the cubic region, FDTD is an accurate numerical scheme for bioheat transfer simulation, which is used as the comparison method in case 1. In order to validate the computational accuracy of ADE for irregular boundary, FEM method was also applied to obtain the base solution in case 2 and 3. 3.1. The validation of ADE method

Based on the discrete scheme of diffusion term Eq. (3), Eq. (1) can be written as vector–matrix type:

@T ¼ AT þ qI  aT @t

3. Results

ð7Þ

where Dt is the time step, subscripts n + 1 and n denote time, T1 and T2 are the intermediate temperature. ADE scheme has both second order accurate in space and time, which is unconditionally stable. For each time-step update, only two simple triangular matrixes solving are used for ADE compared with implicit FDM. The runtime of ADE agrees with the explicit FDM in the same order for each time-step update. However, the larger time step can be used for calculation acceleration compared with explicit FDTD. 2.5. Parallel ADE (PADE) scheme From Eq. (7), it is found that the lower and upper triangular matrix need to be solved explicitly for every time-step update, which hold the most computationally intensive portion in ADE scheme. Considering a two-dimensional voxel O shown in Fig. 2(B), the calculation of T1,o need the value of T1,A, T1,B, Tn,C and Tn,D, where Tn,C and Tn,D are the known quantity, T1,A and T1,B must be updated previously. Inspired by the work of [17], the voxels are aligned along diagonal order (see Fig. 2(B)) so that the parallel method is conveniently applied. Let L = i + j + k, we find from Fig. 2(B) that for all the voxels with the same level L updating T1(L) only depends on T1(L) and Tn(L + 1). For T2(L) calculation, this process is just opposite. Hence it is a good choice to compute T1(L) along ascending order of level L, and updating T2(L) along its descending order. For each level, the parallel strategy is directly used to accelerate T1 (or T2) calculation, which means that all the voxels with level L are uniformly distributed to different threads. For irregular tissue, however, only active voxels are used for computation. In order to balance the computation task for different threads, the active voxels with the same level are reassembled and uniformly distributed to threads. The above parallel strategy is easily implemented based on OpenMP language for multiple processors.

3.1.1. Case 1 The first case is to calculate the transient and steady-state temperate filed in a cubic box with size 50 mm  50 mm  50 mm and voxels number 2013. The boundary conditions are assumed as: adiabatic boundary for x = 0 (and 50 mm) and y = 0 (and 50 mm), Tjz¼0 ¼ 37  C and j@T=@x ¼ hðT 0  TÞjz¼50mm , T 0 ¼ 25  C and h = 10 W m2 °C1. The thermal properties of tissue are assumed as: qc ¼ 4:2 M J m3  C1 and k = 0.5 W m1 °C1. For explicit

FDTD, the allowed maximum time-step is dt ¼ Dh qc=6j, which is considered as base time-step. The initial temperature is given by 37 °C. The tissue temperature would decease to the stable state due to the cooling effects from the boundary z = 50 mm. For the transient bioheat transfer process, the computational results from ADE and FDTD are both given in Fig. 3A. Fig. 3B shows the iteration convergence of ADE for steady-state temperature field simulation. Considering the same physical evolution time, the impact of the time step on the computational time (CPU time) is presented in Fig. 3C. In order to evaluate the parallel efficiency on the regular geometry, the total iteration number of 60,000 with time step dt is considered, the computational time and speedup ratio are both shown in Fig. 3D. 2

3.1.2. Case 2 The second case is to evaluate thermal ablation of a spherical tumor with diameter 20 mm, which locates at the cubic tissue with the same size considered in case 1. Its center position is (25, 25, 35) mm. The boundary at z = 50 mm is modified as T = 15 °C. The blood perfusion of the both two tissues are assumed as xb = 0.0005 s1. The nano-particles (such as gold particles) are often injected to the tumor region to enhance the heat transfer to improve the treatment output. The thermal conductivity of tumor ks is chosen larger than that (kt = 0.5 W m1 °C1) of the surrounding tissue. We first compare the accuracy of the harmonic interpolation method (HIM) based on Eq. (4) and the isotropic interpolation method (IIM, ki+1/2,j,k = ki,j,k [8]) for the steady-state temperature field. FEM method is also applied to obtain the base solution, where the mesh size of 0.1 mm is used to the tumor region and 1 mm for the surrounding tissue. Fig. 4A shows the temperature distribution at the section y = 25 mm for ks = 40kt and Dh = 1 mm through HIM-ADE. The quantitative comparison results are presented in Fig. 4B. For ks = 40kt, the maximum temperature bias is about 0.53 °C for IIM-ADE and 0.1 °C for HIM-ADE compared with the base solution. For ks = 2kt, the both HIM-ADE and IIM-ADE keep well agreement with FEM. In order to modeling the thermal ablation process, the external heat generation Qe = 500,000 W/m3 is imposed on the tumor regions (ks = 2kt). The steady-state temperature field without external heat is considered as initial condition. The influence of the time step on the temperature evolution during thermal ablation process is presented in Fig. 5. Fig. 5A shows the temperature distribution of the section y = 25 mm at heating time t = 600 s for time step 10dt (dt ¼ Dh qc=6js ). The corresponding absolute differences of the temperature for time step dt and 10dt are given 2

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

847

Fig. 3. The comparison results of FDTD and ADE: (A) temperature distribution on line (25, 25, 0:50) mm for different cooling time; (B) the iteration convergence of steadystate temperature computation at the position (25, 25, 47.5) mm; (C) computational time for physical evolution during the transient process; (D) the speedup performance of PADE.

Fig. 4. The influence of the tissue interface on the temperature distribution: (A) the temperature distribution at the section y = 25 mm for ks = 40kt and Dh = 1 mm by HIMADE; (B) the temperature distribution along the line [0:50, 25, 25] mm compared with HIM-ADE, IIM-ADE and FEM for ks = 40kt and ks = 2kt.

in Fig. 5C. Fig. 5B and D, respectively, quantitatively show the impact of time step on the temperature distribution and evolution. The maximum temperature error happens at the tumor center regions, which are 0.15 °C, 0.56 °C and 2.6 °C for time steps of 10dt, 20dt and 30dt at heating time t = 600 s in particular.

3.1.3. Case 3 The third case is to consider the thermal effect of a large blood vessel with diameter 10 mm, which is embedded in the cubic box with same size and outer boundary condition in the first case. The center line position of the vessel passes through the point (35, 25)

848

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

Fig. 5. The impact of the time step on the temperature evolution during thermal ablation process: (A) the temperature distribution of the section y = 25 mm at t = 600 s for 10dt; (B) the temperature distribution of the line [0:50, 25, 25] mm at t = 600 s for different time-steps; (C) the temperature difference of the section y = 25 mm at t = 600 s for dt and 10dt; (D) the temperature evolution of tumor center for different time step.

mm and along y direction. For many biomedical applications, the thermal effect of blood vessel can be assumed as convective boundary condition j@T=@n ¼ hb ðT b  TÞ, with constant blood temperature Tb = 37 °C. The convective heat transfer coefficient hb can be estimated from blood velocity [18], which is here selected as hb = 100 W m2 °C1 and 1000 W m2 °C1 to evaluate the method for irregular boundary treatment. In order to validate the accuracy of ADE, the base solution is also obtained through FEM with conformal tetrahedral meshes, where the mesh size of 0.2 mm is used to the vicinity of vessel and 1 mm set for other domains. Fig. 6A shows the computational temperature contour at the section y = 25 mm by ADE method with Dh = 0.5 mm for hb = 100 W m2 °C1. The weak disorder of the temperature field near the blood vessel is observed. The quantitative results are given in Fig. 6B, where the temperature along line [25, 25, 0:50] mm and crossing the vessel is obtained from FEM and ADE. We also consider a larger hb = 1000 W m2 °C1 [18], which may be more suitable for charactering the thermal effects of large vessel. From Fig. 6C and D, we can find that the computational results from ADE and FEM have good agreement. 3.2. Application in the real tissue In order to furthermore evaluate the performance of PADE for biomedical application, a realistic head model for scalp cooling

evaluation is considered. Fig. 7 shows the tissue structure of the head model, which is reconstructed based on MRI T1 and T2 (Achieva 3T NX, Philips, Netherlands) of a male with 26 years-old. The whole computational domain has size 210 mm  155 mm  235 mm and about 7.65  106 voxels with size 1 mm3, while the number of active voxels is about 3.81  106. The thermal parameters of different tissues can be found in [11]. For the skin surface, the convective boundary condition is considered as j@T=@n ¼ hðT 0  TÞ, T0 = 25 °C and h = 10 W m2 °C1. In order to characterize the thermal effect of the scalp cooling method, the convective boundary condition is also considered for scalp cooling, where h = 3000 W m2 °C1 and cooling temperature Ts = 5 °C (or 10 °C, 15 °C) are assumed. All the initial temperatures are given as 37 °C. The steady-state temperature field is obtained by controlling the updated temperature difference within L2-normal error 107. As comparison, the base state of head temperature without scalp cooling is firstly computed. The results are given in Fig. 8A–C, which show the temperature distribution at sagittal plane y = 72 mm (Fig. 8A), transverse plane z = 140 mm (Fig. 8B) and coronal plane x = 100 mm (Fig. 8C), respectively. For scalp cooling with Ts = 5 °C, the temperature contours are given in Fig. 8D–F with the same position in base-state case (shown in Fig. 8A–C). In order to track the influence of the cooling temperature on the inner head temperature distribution, Fig. 9 shows the temperature at line L1 of

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

849

Fig. 6. The validation of the irregular boundary treatment for ADE considering two cases hb = 100 W m2 °C1 for (A) and (B), and hb = 1000W m2 °C1 for (C) and (D): (A) and (C) the temperature contour of the section y = 25 mm from ADE with Dh = 0.5 mm, (B) and (D) the comparison of temperature computation on the line [25, 25, 0:50] mm with ADE and FEM.

y = 80 mm of coronal plane x = 100 mm, and line of L2 of y = 80 mm of transverse plane z = 140 mm for different cooling temperature.

4. Discussions The body temperature distribution is mainly determined by blood perfusion rate, metabolism and heat transfer with environment. In addition, the effect of the heat conduction is to homogenize the temperature distribution. For the simulation scheme of Pennes’ model, the calculation of heat conduction term is the most time-consuming. The conventional explicit FDTD has to adopt the smaller time step. The implicit scheme can allow the larger time step, however, the complicated and expensive matrix solving is implemented for each time-step update. ADE method presented here eliminates the stability restriction and avoids complicated matrix solving, which only needs to explicitly solve the simple lower and upper triangular matrix for every time-step update. It is partially conducted from Fig. 3A–C. For the transient bioheat transfer process in a simple geometry, ADE method allows the large time step 16dt, while it keeps the same accuracy with FDTD method with time step dt (Fig. 3A). Removing time-step limitation makes simulation more efficient in steady-state temperature computation shown in Fig. 3B, where only iteration number of 760 is obtained by ADE (using a larger time step 800dt) compared with FDTD, which needs iteration number about 30,000. Although ADE has a similar explicit solving process with FDTD, the computational

time of ADE is about twice longer than that of FDTD because of both T1 and T2 calculation in ADE for each time step update. It also indicates that computational speed of ADE (for time step 16dt) is about 8 times faster than that of explicit FDTD (for dt) (see Fig. 3C). For the human tissue, the thermal-conductivity discrepancy is so small that its induced thermal effects are omitted for many medical applications [8,11]. However, this discrepancy would be enlarged in some applications such as nanoparticles-enhanced heat transfer during thermal ablation, which must be carefully treated. For large thermal conductivity jump across the tissue interface, Fig. 4 clearly indicates that Eq. (4) could capture the temperature distribution more accurate compared with the conventional method, where the maximum temperature bias is about 0.53 °C for IIM-ADE and 0.1 °C for HIM-ADE compared with the base solution. In order to further examine the influence of large time step on the transient temperature evolution during the thermal ablation process, a detailed time-step dependence relation is presented in Fig. 5. The results indicate that the maximal temperature error of 0.56 °C induced by time-step 20dt is still accepted compared with the maximal temperature 50.2 °C at heating time t = 600 s. Although time-step 30dt cannot accurately track the history of temperature evolution, it could be used to rapidly predict the steady-state temperature field. Another drawback of FDTD [8] is hard to treat accurately irregular boundary of living tissue structure. An optimized processing pipeline for irregular boundary condition based on the boundarysource method [7] is presented in Section 2.2. The method can

850

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

Fig. 7. The real head model for scalp cooling: (A) the whole head model, (B) the main structure of the head including white matter, gray matter, cerebrospinal fluid, skull, fat and skin; (C) the three-dimensional brain structure.

Fig. 8. The temperature distribution at the base state (A–C) and scalp cooling with Ts = 5 °C (D–F): (A) and (D) sagittal plane y = 72 mm; (B) and (E) transverse plane z = 140 mm; (C) and (F) coronal plane x = 100.

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

851

Fig. 9. The head temperature change induced by scalp cooling, (A) the temperature distribution at line L1 of y = 80 mm of coronal plane x = 100, (B) L2 of y = 80 mm of transverse plane z = 140 mm for different cooling temperature.

directly combine with the medical image segmentation. In addition, the whole process is easily implemented, respectively, the analytical expression of the involved geometrical parameters are given in Appendix A. The accuracy of boundary condition treatment is validated in Fig. 6. We can find that the maximal error of the temperature near vessel is about 0.11 °C for ADE with Dh = 1 mm, within about 0.3% of the maximal temperature, which would decrease to 0.025 °C for the smaller Dh = 0.25 mm (see Fig. 6B). Our test indicated that the temperature gradient near boundary has important effect on the accuracy of temperature estimation. A smaller temperature gradient can make the computation more accurate, which is demonstrated in Fig. 6C and D, where hb = 1000 W m2 °C1 is considered. The above results and discussion have indicated that ADE is more efficient than conventional FDTD method for both transient and steady-state temperature field calculation. In order to furthermore improve the computation efficiency, a special parallel strategy has been developed in Section 2.5. It is noteworthy that ADE and PADE have the same coding complexity. For a simple regular geometry, the speedup ratio of PADE can reach to 10.6 for 12 threads with parallel efficiency 88.3%, which means that the computational time can decrease from 109 min to 10.3 min (see Fig. 3D). For irregular geometry, however, the parallel efficiency would be seriously weakened because of the task-loading imbalance. An improved method is to

reassemble active voxels, which are then uniformly distributed to different threads for the computation task balance. For a real head model, the speedup ratio of PADE can reach to 9.1 for 12 threads with parallel efficiency 75.8% (see Fig. 10). For steady-state temperature computation, the iteration number is another important parameter for computational time. ADE method removes the time-step limitation, while the larger time step cannot reduce the iteration number for irregular tissue due to the inhomogeneous boundary-source terms. We find an optimal time step about 50 s for the real head model, which can obtain a speedup ratio about 30 compared with that for time step 1.2 s (an allowable maximal time step in FDTD). Adopting the optimal time step and parallel scheme, the steady-state temperature field is obtained in a short computational time about 55 s (see Fig. 10). Scalp cooling method has been widely applied to protect the brain and reduce the extent of neurologic induced by brain injuries such as stroke [3], or reduce hair loss during the administration of cytotoxic drugs [2]. Predicting the temperature variation of the inner head is very important for designing scalp cooling strategy. However, the most of previous works focused on a simple spherical model [2], which is not easily to obtain the real temperature field state. Here we have estimated the temperature distribution based on PADE method for a real head model. Compared with the base state (Fig. 8A–C), the temperature of superficial region evidently decreases for a lower cooling temperature Ts = 5 °C (Fig. 8D–F). However, the temperature near brain core is hard to be affected due to the high blood perfusion from the white and gray matter, which can be also conducted from Fig. 9. Benefiting from significant performance of PADE method, we can implement a study on parameterization including cooling temperature and thermal physiology. A detailed discussion will be presented in our near future work. 5. Conclusion

Fig. 10. The parallel efficiency of ADE in the real head model.

In summary, we have developed a fast and accurate parallel numerical scheme for solving bioheat transfer problems. The computational accuracy and efficiency of the method are demonstrated in the test cases including a real head model. The results indicate that PADE method can accurately estimate the temperature filed of the tissue even for irregular structure. In addition, the large time step and efficient parallel performance significantly decrease the computational time, which is important for many biomedical applications, such as developing an efficient treatment planning tool for thermal ablation [20]. PADE method can be directly combined with electromagnetic computation with FDTD method to

852

Z.-Z. He, J. Liu / International Journal of Heat and Mass Transfer 95 (2016) 843–852

where we assume

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8 > jmj ¼ m21 þ m22 þ m23 > < V 1 ¼ m21 =ðmaxð6m2 m3 Þ; 1020 Þ > > : S1 ¼ m1 =ðmaxð2m2 m3 Þ; 1020 Þ

ðA2Þ

The above equation is conducted for unit voxel. So it should be 3

2

scaled by Dh for V and Dh for S. For the other case ni < 0, xi should be firstly transformed by 1  xi to obtain the corresponding to positive component [19]. For voxel with multiple tissue boundary, we could divide the voxel into many sub-voxel so that each sub-voxel has only one segment, and then calculate its surface and volume for each sub-voxel and accumulate them. References

Fig. A1. The schematic diagram of a unit voxel with intersection plane.

evaluate the thermal effect induced by SAR absorption, where both simulation are performed on the same cubic voxel system. Acknowledgments This work is supported by the National Natural Science Foundation of China under Grant Nos. 51476181 and 51376102, Beijing Municipal Science and Technology Project (Z141100000514005), Youth Innovation Promotion Association CAS (2015020), and Key Laboratory of Cryogenics, TIPC, CAS (No. CRYOQN201306). Appendix A. In this appendix, we give the detailed computation of S and V for intersection plane within boundary voxel. The analytical relation connecting linear interface and volume fractions in cubic voxel has been given in [19]. Here, we list the result of volume computation, and add the area equation of intersection plane. Considering a unit voxel shown in Fig. A1 with intersection plane B1B2B3B4, S denotes the area of plane SB1B2B3B4 and V equals to VA1A2A3A4B1B2B3B4. According to discussion in Section 2.2, the intersection plane can be rewritten by m1x1 + m2x2 + m3x3 = a, where three components mi is defined by mi = ni/(n1 + n2 + n3), and have the relations m1 6 m2 6 6 m3. ni denotes component of normal direction defined by Step 3 of Section 2.2. We considered firstly the case that all the three components ni of normal direction are positive. We let m12 = m1 + m2 and m = min(m12, m3). V has odd symmetry with respect to the point (V, a) = (0.5, 0.5) [19]. For a 6 0.5, we can give the analytical expression of S and V as following: 8 3 2 V ¼ 6m1am2 m3 S ¼ jmj 2m1am2 m3 > > >   > > aðam Þ > > S ¼ jmj m2am3  S1 > V ¼ 2m2 m13 þ V 1 > > > > > > V ¼ a2 ð3m12 aÞþm21 ðm1 3aÞþm22 ðm2 3aÞ > > 6m1 m2 m3 > < að2m12 aÞm21 m22 S ¼ jmj 2m1 m2 m3 > > > > a2 ð32aÞþm21 ðm1 3aÞþm22 ðm2 3aÞþm23 ðm3 3aÞ > > V ¼ > 6n1 n2 n3 > > > > 2ðaa2 Þm21 m22 m23 > > S ¼ jmj > 2m1 m2 m3 > > > : 12 V ¼ 2am S ¼ jmj m13 2m3

0 6 a < m1 m1 6 a < m2

m2 6 a < m

m ¼ m3 6 a 6 12 m ¼ m12 6 a 6 12 ðA1Þ

[1] R.B. Roemer, Engineering aspects of hyperthermia therapy, Annu. Rev. Biomed. Eng. 1 (1999) 347–376. [2] F.E.M. Janssen, G.M.J. Van Leeuwen, A.A. Van Steenhoven, Modelling of temperature and perfusion during scalp cooling, Phys. Med. Biol. 50 (2005) 4065–4073. [3] K.R. Diller, L. Zhu, Hypothermia therapy for brain injury, Annu. Rev. Biomed. Eng. 11 (2009) 135–162. [4] C.M. Collins, W. Liu, J. Wang, R. Gruetter, J.T. Vaughan, K. Ugurbil, M.B. Smith, Temperature and SAR calculations for a human head within volume and surface coils at 64 and 300 MHz, J. Magn. Reson. Imag. 19 (2004) 650–656. [5] G.M.J. Van Leeuwen, J.J.W. Lagendijk, B.J.A.M. Van Leersum, A.P.M. Zwamborn, S.N. Hornsleth, A.N.T.J. Kotte, Calculation of change in brain temperatures due to exposure to a mobile phone, Phys. Med. Biol. 44 (1999) 2367–2379. [6] M.Y. Ge, K.J. Chua, C. Shu, W.M. Yang, Analytical and numerical study of tissue cryofreezing via the immersed boundary method, Int. J. Heat Mass Transfer 83 (2015) 1–10. [7] Z.Z. He, X. Xue, J. Liu, An effective finite difference method for simulation of bioheat transfer in irregular tissues, J. Heat Transfer Trans. ASME 135 (2013) 071003–071010. [8] T. Samaras, A. Christ, N. Kuster, Effects of geometry discretization aspects on the numerical solution of the bioheat transfer equation with the FDTD technique, Phys. Med. Biol. 51 (2006) 221–229. [9] E. Neufeld, N. Chavannes, T. Amaras, N. Kuster, Novel conformal technique to reduce staircasing artifacts at material boundaries for FDTD modeling of the bioheat equation, Phys. Med. Biol. 52 (2007) 4371–4381. [10] V. Singh, A. Roy, R. Castro, K. McClure, R. Dai, R. Agrawal, R.J. Greenberg, J.D. Weiland, M.S. Humayun, G. Lazzi, On the thermal elevation of a 60-electrode epiretinal prosthesis for the blind, IEEE Trans. Biomed. Circuits Syst. 2 (2008) 289–300. [11] G. Carluccio, D. Erricolo, S. Oh, C.M. Collins, An approach to rapid calculation of temperature change in tissue using spatial filters to approximate effects of thermal conduction, IEEE Trans. Biomed. Eng. 60 (2013) 1735–1741. [12] P. Zeng, Z.S. Deng, J. Liu, Parallel algorithms for freezing problems during cryosurgery, Int. J. Inf. Eng. Electron. Bus. 2 (2011) 11–19. [13] S. Leung, S. Osher, An alternating direction explicit (ADE) scheme for timedependent evolution equations, UCLA Comput. Appl. Math. Rep. 11 (2005) 1– 21. [14] A. Fell, G. Willeke, Fast simulation code for heating, phase changes and dopant diffusion in silicon laser processing using the alternating direction explicit (ADE) method, Appl. Phys. A Mater. Sci. Process. 98 (2010) 435–440. [15] V. Lempitsky, Surface extraction from binary volumes with higher-order smoothness, IEEE Conf. Comput. Vision Pattern Recogn. 2010 (2010) 1197– 1204. [16] C. Li, R. Huang, Z. Ding, J.C. Gatenby, D.N. Metaxas, J.C. Gore, A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI, IEEE Trans. Image Process. 20 (2011) 2007–2016. [17] F. Gibou, C. Min, On the performance of a simple parallel implementation of the ILU-PCG for the Poisson equation on irregular domains, J. Comput. Phys. 231 (2012) 4531–4536. [18] L. Consiglieri, I. dos Santos, D. Haemmerich, Theoretical analysis of the heat convection coefficient in large vessels and the significance for thermal ablative therapies, Phys. Med. Biol. 48 (2003) 4125–4134. [19] R. Scardovelli, S. Zaleski, Analytical relations connecting linear interfaces and volume fractions in rectangular grids, J. Comput. Phys. 164 (2000) 228–237. [20] Z.Z. He, J. Liu, An efficient thermal evolution model for cryoablation with arbitrary multi-cryoprobe configuration, Cryobiology 71 (2015) 318–328.