Sensitivity analysis of temperature in heated soft tissues with respect to time delays

Axisymmetric tissue region heated by an external heat flux is considered. The mathematical model is based on the dual-phase lag equation supplemented by appropriate boundary and initial conditions. This equation, in relation to the Pennes’ equation, has two additional parameters, namely the relaxation time and the thermalization time. The aim of this research is to estimate the temperature changes due to changes of these parameters. To achieve this, sensitivity analysis methods are used. The basic problem and additional ones related to the sensitivity functions are solved using the implicit scheme of the finite difference method. The performed computations show that the temperature changes caused by changes in the relaxation and thermalization times are larger for higher values of the external heat flux and shorter times of its action.


Introduction
In 1948, Pennes [1] presented an equation of heat transfer in human tissue, which included the effects of blood flow on tissue temperature. This equation, known as the 'traditional' or 'classic' or 'Pennes' bioheat equation [2], has the following form where c, ρ, λ are the specific heat, density and thermal conductivity of tissue, respectively, T denotes the temperature, x, t are the geometrical coordinates and time. The source term Q(x, t) is as follows where w b is the blood perfusion rate, c b is the specific heat of blood, T b is the arterial blood temperature, and Q m is the metabolic heat source. Pennes checked the validity of this approximation by comparing temperatures predicted by his equation with experimentally measured temperatures in the human forearm.
It should be noted that the Pennes equation is in fact the Fourier one in which the source components related to perfusion and metabolism are introduced. The Pennes equation (like the Fourier equation) is based on the Fourier law where q(x, t) is the heat flux and ∇T (x, t) is the temperature gradient. Another model of the bioheat transfer is based on the Cattaneo-Vernotte equation [3,4]. This model takes into account the nonhomogeneous inner structure of biological tissues, which causes that the thermal wave propagates through the medium with a finite velocity. The heat flux delay relative to the temperature gradient is introduced where τ q = a/C 2 is the relaxation time, a = λ/(cρ) is the thermal diffusivity of tissue, C is the velocity of thermal wave in the medium. Taking into account the first-order Taylor expansion for q one obtains the Cattaneo-Vernotte equation in the form Experimentally determined values of τ q can be found, among others, in [5,6], but they are questioned by some researchers, e.g., [7]. Mitra et al. [5] reported that relaxation time τ q in processed meat could be as large as 16 s. Kaminski [6] reported that τ q ranges from 10 to 50 s in his experiment with materials with inhomogeneous inner structures. Roetzel et al. [8] in a similar experiment estimated the value of this parameter at 1.77 s. The next model of bioheat transfer is based on the following relationship [9,10] q where τ T is the thermalization time.
Taking into account the first-order Taylor expansions for q and T one obtains the following form of bioheat transfer equation Equation (9) is called dual-phase lag equation (DPLE) in which τ q is the phase lag in establishing the heat flux and associated conduction through the medium, and τ T is the phase lag in establishing the temperature gradient across the medium during which conduction occurs through its small-scale structures [9]. It is also worth mentioning the model based on the theory of porous media proposed by Zhang [11] and are the effective thermal conductivity and effective heat capacity, respectively (subscript b denotes blood subdomain), Q m and Q b are the metabolic heat sources in the tissue and blood sub-domain, respectively, G is the coupling factor, and ε denotes the porosity (the ratio of blood volume to the total volume).
The phase lags for heat flux and temperature gradient and the coupling factor are defined as [11] and while where A [m 2 /m 3 ] is the volumetric heat transfer area between tissue and blood, α is the heat transfer coefficient, and w b is the blood perfusion rate. As can be observed, the phase-lag times are expressed in terms of the properties of blood and tissue, the interphase convective heat transfer coefficient α and the blood perfusion rate w b . These parameters are closely related to the porosity, and larger porosity values determine larger values of these parameters. Zhang [11] computed and argued that reasonable relaxation time τ q and thermalization time τ T should range from 0.464 to 6.825 s.
As can be seen, the key parameters in the non-Pennes bioheat equations are the delay times τ q and τ T . Experimental studies show that delay times vary widely and are different for different types of biological tissues. The purpose of the research is to estimate the temperature changes in heated tissues due to perturbation in time delays. To achieve this, sensitivity analysis methods can be used [12][13][14][15]. In the direct approach of sensitivity analysis, the basic problem described by the dual-phase lag equation and appropriate boundaryinitial conditions is differentiated with respect to the parameter under consideration. Next, these both problems are solved using (in the case considered) the numerical methods. Finally, based on the solutions obtained, the impact of changes in individual parameters on the tissue temperature distribution is assessed. In this way one can analyze the influence of the parameters appearing in the mathematical model on the final form of problem solution; in other words, which of these parameters are more or less important.

Formulation of the problem
Axisymmetric tissue region heated by an external heat flux, as shown in Fig. 1, is considered. To describe the temperature field in the tissue dual-phase lag equation (9) for constant thermophysical parameters is applied where Introducing Eq. (2) into Eq. (16) one obtains On the fragment of the tissue surface 1 : r ≤ r D , z = 0 ( Fig. 1), for t ≤ t e , where t e is the exposure time, the Neumann condition is assumed [16] where q 0 is a constant value. For t > t e : q b (r, 0, t) = 0. On the remaining surfaces 2 of the domain, the no-flux condition is accepted. It should be noted that the dual-phase lag model requires the adequate transformation of boundary conditions which appear in the typical macroheat conduction models. In the case considered one has [17] − λ n · ∇T (r, where n is the normal outward vector. The initial conditions are also given (c.f. Eq. (2)) where T p is a constant initial temperature.

Sensitivity analysis
To estimate the changes of tissue temperature due to relaxation and thermalization times, the direct approach of sensitivity analysis is applied [12,15]. This approach is related to the differentiation of governing equations with respect to the parameter considered. Thus, the differentiation of Eq. (18) with respect to the relaxation time τ q gives (the arguments of the T function are omitted here) Using the Schwartz theorem one can obtain where U = ∂ T /∂τ q is the sensitivity function. Boundary condition (20) and initial condition (21) are also differentiated and In a similar manner, an additional problem related to the sensitivity function with respect to the thermalization time τ T can be obtained and where W = ∂ T /∂τ T .

Method of solution
The formulated problem is solved using the implicit scheme of the finite difference method. The temperature where t is the constant time step. The uniform spatial grid of dimensions n × n with a constant grid step h, the same in both directions, is introduced. For the node i, j the temperature is denoted as The following approximate form of Eq. (18) is proposed or where (c.f. Eq. (17))

Boundary condition (20) is written in the form
After mathematical manipulations one has On the remaining surfaces (the no-flux boundary condition) the temperatures are calculated using the equations The system of Eqs. (32), (36)-(39) for transition t f −1 → t f is solved using the iterative method. The iteration process is continued until the condition T ≤ δ for each node has been fulfilled, where k is the number of iteration. It should be noted that the presented algorithm is unconditionally stable [18].

Results of computations
The cylindrical tissue domain (R = 0.02 m, Z = 0.02 m) at the initial temperature T p = 37 • C is considered.
The following values of parameters are assumed: thermal conductivity of tissue λ = 0.5 W/(mK), specific heat c = 4000 J/(kg K), density ρ = 1000 kg/m 3 , blood perfusion rate w b = 0.53 kg/(m 3 s), specific heat of The problem is solved under the assumption that the grid step is equal to h = 0.0002 m (n = 100) and next h = 0.0004 m (n = 50), while the time step is equal to t = 0.0005 s. The system of equations resulting from the FDM algorithm is solved with the tolerance δ = 10 −12 . In Fig. 2 the temperature histories at the selected nodes A(0, 0.0004 m), B(0, 0.002 m) and C(0, 0.004 m) from the domain (for two variants of grid steps) are shown. It can be seen that the temperatures are lower with a denser mesh, especially in the most heated node A.
The time step was also changed, and it turned out that for the first mesh (h = 0.0002 m) the differences in the results for time steps t = 0.0005 s and t = 0.005 s are very small. Therefore, further computations were performed for n = 100 and t = 0.005 s.
In Fig. 3 the temperature field in the fragment of the domain r ≤ R/2, z ≤ Z /2 after 50 s and 100 s is shown. Figures 4 and 5 illustrate the distribution of sensitivity functionsU and W , respectively.
Knowledge of the sensitivity functions allows one, among others, to estimate the changes of temperatures due to the perturbations of the parameters considered. At first, the influence of change of single parameter on where τ q and τ T are the assumed perturbations of parameters. It should be noted that the estimate performed here is a linearization and thus only admissible for relatively small τ q and τ T . In Figs. 6 and 7 the changes of temperature at the points A, B and C under the assumption that τ q = 0.3τ q and τ T = 0.3τ T are shown.
As can be seen, the maximum temperature change is observed at the point A, of course. The maximum of function D 1 T at this point corresponds approximately to the maximum value of boundary heat flux (19) for t = t e /2 = 50 s (Fig. 6). At the cooling stage t > 50 s, the value of the D 1 T function decreases reaching a minimum after 115 s and then tends to zero. The temperature changes at the points B and C are similar, with the maximum and minimum values being lower and appearing later.
In the case of the temperature change caused by a perturbation in the thermalization time (Fig. 7), during the heating stage (t ≤ 50 s), at the point A, the function D 2 T tends to a minimum and then reaches its maximum value after 115 s. Temperature changes at the points B and C are similar, but the minimum and  It should be noted that the maximum temperature change in the whole domain due to the perturbation of relaxation time τ q is less than 1.67 • C, while the temperature change due to the perturbation of thermalization time τ T is less than 1.28 • C.
It is also possible to estimate the changes of temperature due to the simultaneous perturbations of two parameters, namely In Fig. 8 the courses of function DT at the points A, B and C are shown. As can be seen, the maximum temperature change at the point A is about 0.95 • C. The computations were repeated for q 0 = 30,000 W/m 2 , t e = 100 s and q 0 = 30,000 W/m 2 , t e = 50 s, respectively (c.f. boundary condition (19)). In Figs. 9 and 10 the changes of temperature due to the simultaneous perturbations of relaxation time and thermalization time (Eq. (51)) are presented. As shown in Fig. 9, for larger value of boundary heat flux, the changes in temperature are also larger. For example, the maximum temperature change at the point A is about 1.9 • C. Moreover, for a shorter exposure time (Fig. 10) these differences are even larger and the maximum temperature change at point A is about 2.1 • C.

Conclusions
The tissue heating process described by the dual-phase lag equation supplemented with appropriate boundaryinitial conditions has been considered. Sensitivity analysis with respect to delay times has been done. To solve the basic problem and additional ones the implicit scheme of the finite difference method has been used.
First, the influence of a single parameter on the changes in tissue temperature was examined. It turned out that slightly larger changes in temperature are observed in the case of disturbed relaxation time than in the case of disturbed thermalization time. Then, the changes in temperature caused by simultaneous perturbation of both parameters were analyzed. These changes depend on the value of the external heat flux and the exposure time. It turned out that the higher the value of the external heat flux and the shorter the exposure time, the larger the temperature changes, the detail can be found in the previous sub-chapter.
It would seem that the temperature changes are relatively small. However, it should be remembered that in artificial hyperthermia treatments, which consists in increasing the temperature of the diseased tissues to 44-45 • C, the difference of 2 • C is very significant.