NUMERICAL ANALYSIS OF ARTIFICIAL HYPERTHERMIA TREATMENT

. This paper presents numerical modelling of artificial hyperthermia treatment. Presented model takes into account not only the temperature distributions but also the thermal dose parameter. Obtaining of temperature distributions takes advantage of the generalized dual phase lag equation. For computer calculations the parallelized algorithm was prepared.


Introduction
From the medical point of view hyperthermia is the sudden, rapid rise in a body temperature. Artificial hyperthermia is a treatment, in which the body temperature is purposeful raised, usually to 42 -46 °C. There are three types of artificial hyperthermia: local, regional and whole-body. In this paper, only the local is considered. Local artificial hyperthermia is usually used as cancer treatment often associated with chemoor radiotherapy. When this treatment is used unsupported, in cancer cells it causes lack of oxygen and nutrients, what leads to the apoptosis. Heat shock causes inducing of the heat shock proteins. Rise of temperature results in better blood supply to the organ and therefore drug accumulation. Chemical reactions are faster at higher temperatures. Artificial hyperthermia associated with the radiotherapy perpetuates damage of DNA.
It also should be noted, that the biological tissue is the material with particular nonhomogeneous inner structure and interwoven by blood vessels (Fig. 1). Sensitive influence on the temperature distribution has a volume of the blood vessels and blood velocity. The bioheat transfer process is multiscale, therefore it is necessary to consider delays of heat flux and temperature gradient [1].

Fig. 1. Tissue model
To prevent damage of healthy tissue, as well heating the considered area to desirable temperature, the ability to predict the temperature distribution accurately, in a short calculation time is very important.
A model, which allows one to take into account the tissue porosity and phase lags depending on the parameters of tissue is the generalized dual phase lag model [11]. It should be pointed out that the comparison of various bioheat transfer models was done by authors in [8,9].
It should be remembered that the degree of tissue destruction depends not only on the temperature, but also the exposure time and can be described mathematically by means of the thermal dose parameter [14].

Generalized dual phase lag equation
The tissue, as shown in figure 1, can be treated as a porous medium divided into two regions: the vascular region (blood vessel) and the extravascular region (tissue) [4,11]. To describe temperature field in the heating regions (blood (1) and tissue (2)) the two-equation porous model [12] can be applied where ε denotes the porosity (the ratio of blood volume to the total volume), α is the heat transfer coefficient, v is the blood velocity, A is the volumetric transfer area between tissue and blood, c is the specific heat, ρ is the density, λ is the thermal conductivity, T denotes temperature, t is the time, w is the blood perfusion rate, Q m is the metabolic heat source and Q ex is the capacity of internal heat sources associated with the external heating of tissue [9] while subscripts t and b represent tissue and blood, respectively. Adding both (1) and (2) equations, the following equation can be obtain In this paper is assumed that the coupling factor is equal to G = Aα + wc b and also that before reaching the equilibrium temperature of tissue and blood, the blood temperature changes according to the Minkowycz hypothesis [10]   Based on (4) the temperature of tissue is described as follows Using (5) and (3), after some mathematical operations the equation for blood temperature can be written in the form Now, the effectiveness parameters can be introduced: and Assuming the following form of relaxation time and the thermalization time the equation (6) can be written as follows In equation (11) the unknown is the blood temperature. To determine the equation where only unknown is the tissue temperature the dependence (5) should be transform Based on (11), (12) and after some mathematical operations the equation for tissue temperature is described as follows

Concept of thermal dose
Knowledge of timedependent temperature field during the thermal treatment allows one to determine the thermal dose TD in terms of equivalent minutes at temperature 43°C. In particular, the following equation should be taken into account [14] where t 0 , t F correspond to the initial and final times, respectively, T f is the temperature at the point considered for time t f , Δt is the time step, R = 0 for T ≤ 39 °C., R = 0.25 for 39 °C. < T < 43 °C. and R = 0.5 for T ≥ 43 °C. The TD value required for total necrosis in a case of muscle tissue (this type of soft tissue is considered here) is equal to TD = 240 minutes [14].

Formulation of the problem
Assumed model is shown in figure 2. The domain of healthy tissue Ω 1 is a cube with edge length of 0.05 m and centrally located subdomain of the tumor Ω 2 with edge length of 0.01 m. The considered domain includes the blood vessels arranged in the direction of the X axis.
The thermophysical parameters of tumor and healthy tissue are assumed to be the same, so only one equation describing the temperature field in domain Ω = Ω 1  Ω 2 is considered. The external heating of tissue is a constant function where Q 0 is constant nonzero component and t ex is duration of heating (exposure time).

Methods of solution
Assuming constant value of metabolic heat source and using formula (15) the equation (13) can be written in the form This equation is supplemented by boundary condition where n is the normal outward vector [3]. The initial conditions are as follows , , , where ∆t is the time step [7]. Then, for time t f = f ∆t (f ≥ 2) the following approximate form of equation (16) can be proposed For simplification of notation the subscripts t and e are here omitted. The uniform grid of dimensions n × n × n is introduced and then the finite difference equation for internal node (i, j, k) has the following form [6]  where:  while s = f -1 or s = f -2 and h is the constant grid step. Finally, the temperature at the node (i, j, k) is calculated from ,, It should be pointed out that in the case of explicit scheme application a criterion of stability should be formulated.
The solving system is stable if the coefficients in the difference equations (22) for time t f -1 are non-negative. Hence it results that the following coefficient must be positive

CPU parallel algorithm
All the time steps must be performed consecutively so the time loop can't be divided into the parallel calculation. Temperature calculations at all nodes in each time step are executed in three nested loops: in the x direction, in the y direction and in the z direction. These calculations can be easily divided into the parallel because all temperatures at the points are calculated based on the f -1 and f -2 time steps. In figure 3 the example of parallelization of CPU calculations is shown. The number of threads "q" depends on the CPU cores.

GPU parallel algorithm
To achieve acceleration of computing time the CUDA technology of NVidia was implemented in computer program. This platform allows to use the graphics processing unit (GPU) for scientific computing. Unlike the CPU, the graphics processor is made up of hundreds of thousands cores (Fig. 4). Each of these cores may perform calculations independently. Fig. 4. Model of CPU and GPU [5] The main limitation of graphics cards is the speed of copying data from the host to the device and the access time to data in the memory. The idea of using the GPU to accelerate the calculations associated with the finite difference method application is based on the fact that in each successive time steps the calculated data are new data for next time step, so it is not necessary to copy these data from the host to the device. The program algorithm is presented in figure 5.
The next very important factor which haves major impact on application efficiency is the idea of using shared memory. The global memory of the device has long access latencies and finite access bandwidth. Unlike to the previously, shared memory can be accessed at very high speed. The biggest problem in use of shared memory is limited amount of this memory [13]. Data for time steps f, f -1 and f -2 are stored in the memory as one dimensional arrays. In chosen approach each thread block loads a tile of data from those arrays. Proper choice of tile dimensions required some experimentations. It was decided to use a 4 × 4 × 25 block size. Then for calculated all nodes temperature, the two arrays, in each threads blocks, from previous time steps of size 6 × 6 × 27 are required ( fig. 6).
The part of the code of the device kernel function, which copies the data from the global to the shared memory and calculates node temperature is shown on figure 7. After determining the temperature, for every node the thermal dose was also calculated.

Results
In numerical computations the following values of parameters have been assumed: thermal conductivity of blood λ b , thermal conductivity of tissue λ t = 0.5 W/(mK), blood density ρ b = 1060 kg/m 3 , tissue density ρ t = 1000 kg/m 3 , specific heat capacity of blood c b = 3770 J/(kgK), specific heat capacity of tissue c t = 4000 J/(kgK), metabolic heat source (of tissue and blood) Q mb = Q mt = 250 W/m 3 , blood temperature T b = 37°C, initial temperature T p = 37°C, porosity ε = 0.0137 and G = 27097.8 W/(m 3 K). The values of phase lag times τ q and τ T were determined using formulas (9) and (10). The spatial discretization creates 500×500×500 nodes and time step is equal to ∆t = 0.01 s. Following heating condition have been taken into account [9]: 35 s heating with a power density of 1 MW/m 3 .
The program has been running on a computer with the processor Intel Core i7-3960X and the graphic card GeForce GTX 680. The processor has the six cores, each with two threads and the clock speed 3.3 GHz. The graphics processing unit has 1536 CUDA cores with base clock 1006 MHz and 2048 MB of global memory.
In figure 8 the temperature history at the central node of cube is presented. One can see that in the cube the maximum temperature 44°C. occurs and the temperature above 43°C is maintained only by 18 seconds.  Figure 9 illustrates the temperature distribution at the cross section (z = 0) after 10 second. The temperature above 37°C is just inside the tumor region.  In Figure 10 the temperature distribution in the cross section of the cube after 35 second is shown. In healthy tissue the temperature above 37°C appeared. After 60 seconds the temperature began to decrease - figure  11, but in the healthy tissue, there is still the temperature above 37°C. Also, more of the area Ω 1 is occupied by the elevated temperature. After 100 seconds the temperature in the entire domain is less than 39°C (figure 12).   Most important thermal dose distribution is after 100s because then the temperature decreases under 39°C and then the coefficient R in the thermal dose concept is equal to zero. This distribution is presented in figure 15. It should be emphasized that all values of TD above 0 minutes are still inside the tumor region. This means that the healthy tissue is heated, but it does not receive a significant thermal dose.
As previously mentioned, the time step is equal to ∆t = 0.01s, while the analysis time is equal to 100s. Thus, the amount of time steps is equal to 10,000. Assumed number of nodes is equal to 500 × 500 × 500, so the number of temperatures which is necessary to calculate in each time step is equal to 125,000,000.
Computer program, which used only the CPU, all the calculations has been performed during 37 minutes and 43 seconds, at 95% of CPU utilization (Fig. 16).
Computer program which used not only CPU but also GPU, the same calculations has been performed during 14 minutes and 1 second. As can be seen, the acceleration of calculations is very significant.  Table 1 provides a comparison of calculation times for different spatial discretizations. In all variants the same time step has been assumed. When changing the discretization for 100 × 100 × 100 the tile sizes also have to be changed (2 × 2 × 25), of course. The maximum difference between the results of the GPU and CPU calculations was below 1.6•10 -9 .

Results
In this paper the bioheat transfer process in three dimensional domain including the healthy tissue and the tumor region has been considered. Some simplifications were adopted, for example: very regular tumor shape and heating only in domain Ω 2 .
As can be seen in figures 9 -12, the temperatures above 37°C. occur not only in the tumor, but also in the area of healthy tissue. During the treatment it is very important to prevent damage of healthy tissue and to provide adequate thermal dose in the tumor region. For a patient, the long duration of heating at high temperature will induce a feeling of discomfort and pain [8]. Based on the received thermal dose distribution it can be seen that not only the temperature is important, but also the exposure time. Figures 14 and 15 show that the thermal dose values above zero minutes are cumulated inside domain Ω 2 .
Using CUDA platform significantly speeds up the calculations. Calculations on graphics processing units should be implemented to the program, which repeatedly performs the complex analysis.