## Abstract

Parabolic heat equation based on Fourier's theory (FHE), and hyperbolic heat equation (HHE), has been used to mathematically model the temperature distributions of biological tissue during thermal ablation. However, both equations have certain theoretical limitations. The FHE assumes an infinite thermal energy propagation speed, whereas the HHE might possibly be in breach of the second law of thermodynamics. The relativistic heat equation (RHE) is a hyperbolic-like equation, whose theoretical model is based on the theory of relativity and which was designed to overcome these theoretical impediments. In this study, the three heat equations for modelling of thermal ablation of biological tissues (FHE, HHE and RHE) were solved analytically and the temperature distributions compared. We found that RHE temperature values were always lower than those of the FHE, while the HHE values were higher than the FHE, except for the early stages of heating and at points away from the electrode. Although both HHE and RHE are mathematically hyperbolic, peaks were only found in the HHE temperature profiles. The three solutions converged for infinite time or infinite distance from the electrode. The percentage differences between the FHE and the other equations were larger for higher values of thermal relaxation time in HHE.

## 1. Introduction

Energy-based ablative techniques (EBATs) employ certain types of energy (e.g. radiofrequency (RF), microwave, laser, ultrasounds) to destroy or irreversibly alter a target zone of biological tissue. High-temperature EBATs raise temperature tissue over 50°C to achieve these results. For instance, RF ablation (RFA) is a high-temperature EBAT that uses RF electrical currents to produce coagulative necrosis. It has been broadly used in many procedures such as eliminating cardiac arrhythmias [1] and destroying tumours [2]. Theoretical modelling has been proposed as a fast and inexpensive approach to developing new EBATs and studying the electrical and thermal phenomena involved in the ablative process [3]. To date, most of these models have employed the bioheat equation proposed by Pennes [4] (Fourier's theory heat equation, FHE), in which the heat conduction term is based on Fourier's theory and the heat flux ** q**(

**x**,

*t*) is

*k*is thermal conductivity and

*T*(

**x**,

*t*) temperature. We thus obtain the classical parabolic heat equation

*S*(

**x**,

*t*) is the heat source and

*α*is thermal diffusivity. However, FHE assumes an infinite thermal energy propagation speed, which is physically inadmissible. Although this approach might be suitable to model most EBATS, it has been suggested that under certain conditions, such as very short heating times (miliseconds or shorter), a non-Fourier model should be considered [5–8]. About this non-Fourier approach, Morse & Feshbach [9] proposed a heat equation similar to the telegraph equation, known as Maxwell's equation because it resembles the equation of the propagation of an electromagnetic field, e.g. light, and which is basically a wave equation with a damping term:

*C*is the speed of heat propagation in the medium. Tisza [10] and Landau [11] predicted that the speed of heat in a medium can be different from the speed of sound and called it the

*second sound*[12]. Different arguments suggest that the speed of heat is at least one order of magnitude less than the speed of sound. But within non-homogeneous materials, such as biological tissues or in the presence of severe plastic deformation or a high strain rate, the speeds of heat are much lower than in metals and perfect structures. This is because in metals, both heat and sound are propagated by mechanisms, such as lattice vibrations, which are all decelerated by the structural imperfections of biological tissues and non-metal materials. In addition, some studies have shown that in many cases the values of

*C*and

*α*imply that the first term of equation (1.3) (second-order derivative) does not have to be necessarily negligible with respect to the second term (first-order derivative) [12,13]. With the same purpose, Cattaneo [14] and Vernotte [15] proposed a modification of Fourier's linear flux as follows:

*τ*of the heat conduction for the medium, that is, there is a response time lag between the heat flux vector and the temperature gradient. Tzou [16] provided a macroscopic formulation of the problem considering that due to this delay the temperature gradient at time

*t*results in a heat flux at a later time

*t*+

*τ*

**(**

*q***x**,

*t*+

*τ*) for its first-order Taylor expansion. Then we obtain the corresponding modification of the heat flux law by direct integration:

*S*(

**x**,

*t*)=0), this equation is identical to the Maxwell equation taking

*C*

^{2}=

*α*/

*τ*.

However, the formulation of HHE has a problem. There are differences between FHE and HHE from a thermodynamic point of view. When they are analyzed within the framework of the traditional irreversible thermodynamics, which is based on the assumption of local equilibrium, FHE yields an expression for the entropy production rate which is always non-negative. By contrast, equation (1.4) yields an expression for the entropy production rate which, under certain circumstances [18], can assume negative values, which contradicts the second law of thermodynamics. It should be emphasized that a generally accepted theory of irreversible thermodynamics, in which local equilibrium is released, is not yet available. Of course, this does not necessarily mean that HHE violates a non-equilibrium second law, if such a law can be established, but it cannot be ruled out. Moreover a non-equilibrium equation that yields an expression of the entropy production rate which is non-negative is likely to satisfy a more general non-equilibrium second law.

To overcome the problems involved in FHE and HHE, Ali & Zhang [12] proposed a theoretical model based on a weak interpretation of the theory of relativity. In this setting, relativistic does not mean ‘concerned with objects moving at speeds comparable to the speed of light’, but ‘concerning objects that move or propagate at speed comparable with a limiting speed characteristic of the field or medium involved, denoted by *C*’. For example, the speed of the photon in electromagnetism, speed of light in a vacuum, speed of the graviton in cosmology, speed of sound in wave propagation in fluids and the limiting speed of heat in heat conduction problems in gas and solid media. Consequently, this theory ensures a finite thermal energy propagation speed.

Ali and Zhang's relativistic model adds the time coordinate *t* to the three Euclidean coordinates (*x*,*y*,*z*). In relativistic physics, it is usual to work in a four-dimensional pseudo-Euclidean space–time (*ζ*,*x*,*y*,*z*), where *ζ* is the so-called ‘spacelike-time’, has a length dimension and, moreover, *ζ*=*iCt*, where i is the imaginary unit constant and *C* is a constant with a velocity dimension, ensuring in this way that *ζ* remains orthogonal to the other three spatial dimensions, with respect to Minkowski metrics

In this setting, Fourier's flux law takes the form
*C* (assumed to be constant) is the speed of heat propagation in the medium.

This RHE in a four-dimensional space–time, turns out to be identical to the Maxwell equation, is derived without any microscopic or material-specific considerations, and overcomes thermodynamic contradictions, as has been shown by Ali & Zhang [12] . Moreover, except for the terms involving the heat source *S*(**x**,*t*), RHE (1.10) is mathematically identical to HHE (1.7), but the HHE heat flux is based on (1.4) and the RHE heat flux on (1.7). That is why RHE has good thermodynamic behaviour.

Since this difference in the flux computations involves different mathematical expressions in the boundary conditions of RHE and HHE, the corresponding solutions could also be different. It is necessary to emphasize that Ali & Zhang [12] only proposed the mathematical formulation of the relativistic problem and they did not solve it in any specific case. Our objective was therefore to solve the relativistic problem for the concrete case of thermal ablation of biological tissues (including an internal heat source) in order to compare the temperature distributions in the tissue computed theoretically from the three heat equations (Fourier, hyperbolic and relativistic). This case was considered since experimental data have suggested that the thermal relaxation time in biological tissues could be significant [13]. The analytical solutions for the HHE and FHE were previously obtained [19]. In this study, we first obtained the analytical solution to RHE for RFA, and then compared the corresponding solutions in terms of computer temperature distributions.

## 2. Material and methods

### (a) Description of the model for relativistic heat equation

In order to achieve an analytical solution, we considered a model based on a spherical electrode of radius *r*_{0} completely embedded and in close contact with the biological tissue, which has an infinite dimension. As a result, the model has a radial symmetry, and a one-dimensional approach is possible. We can then write *T*(*r*,*t*) instead of *T*(**x**,*t*). For thermal ablations based on RFA, the source term in the tissue is the Joule heat produced per unit volume of tissue, *S*(*r*), which can be expressed as follows [20]:
*P* is the total applied power (W). In our case, we use spherical coordinates and all the involved functions are independent of the spherical angles. The resulting RHE is

The initial and boundary conditions are
*T*_{0} is the initial temperature. To write the boundary condition in *r*=*r*_{0}, we shall adopt a simplification used in Erez & Shitzer [20], assuming the thermal conductivity of the electrode to be much greater than that of the tissue (i.e. assuming that the boundary condition at the interface between electrode and tissue is mainly governed by the thermal inertia of the electrode). At every time *t*, the modulus of the heat flux along the surface electrode by unit time is
*ν*_{0} and *c*_{0} are the density and specific heat of the active electrode, respectively. Then the boundary condition in *r*=*r*_{0} must be

### (b) Comparison of models

Once the RHE analytical solution was obtained, we compared the three models (FHE, HHE and RHE) by plotting temperature distributions and evolution with Mathematica v. 6.0 (Wolfram Research Inc. Champaign, IL, USA). Even though some results were plotted using dimensionless variables, we had to consider specific values for some parameters of the model: electrode radius *r*_{0}=1.5 mm, constant power *P*=1 W, initial temperature *T*_{0}=37°C, tissue density *ν*=1200 kg m^{−3}, tissue-specific heat *c*=3200 J kg^{−1} K^{−1}, tissue thermal conductivity *k*=0.7 W m^{−1} K^{−1}, electrode density *ν*_{0}=21 500 kg m^{−3} and electrode-specific heat *c*_{0}=132 J kg^{−1} K^{−1}. These values have been used in previous modelling studies [19]. The dimensionless thermal relaxation time of the biological tissue for HHE was *λ*=1.2963, which is equivalent to *τ*=16 s. This value was initially suggested by Mitra *et al.* [13] in processed meat. In RHE appears a parameter *C* which depends on the medium properties. To compare with the HHE, we chose *C* in such a way that *C*^{2}=*α*/*τ* so that the coefficients of the governing equation would be the same in both models, even though this procedure did not ensure the exact replication of the physical reality. It is noteworthy that with the values of the parameters considered in this study, *C* takes a theoretical value as small as *C*=9.07×10^{−5} m s^{−1}, in good agreement with the qualitative indications given in Ali & Zhang [12]. To date, there are no experimental data about the value of *C*.

## 3. Results

### (a) Relativistic heat model

To solve analytically the RHE for a spherical electrode, we chose dimensionless variables
*B*:=*A*,*α*/*r*_{0}.

The dimensionless problem was to find *V* (*ρ*,*ξ*) such that

Taking Laplace transforms *ξ*, we obtained

Using the new function *M*_{1}(*s*) and *M*_{2}(*s*) are functions to be determined in order to verify the boundary conditions.

To satisfy (3.3), we needed to choose

The computation of the inverse Laplace transforms *H*(*x*) is the Heaviside function, *x*≤0 and *x*>0,
*s*_{1} is defined in appendix A.

### (b) Hyperbolic and Fourier heat models

The dimensionless problem of the FHE is to find *V*_{F}(*ρ*,*ξ*) such that
*m*=*ν*_{0}*c*_{0}/*νc* is the dimensionless electrode thermal inertia. The solution of this problem was previously presented in López Molina *et al.* [19], and we can express it as follows:

Likewise, the dimensionless problem for the HHE is to find *V*_{H}(*ρ*,*ξ*) such that
*H*(*ξ*) is the Heaviside function, *δ*(*ξ*) is the Dirac's delta function, *m* is the dimensionless electrode thermal inertia as above. The analytical solution of HHE was also obtained in López Molina *et al.* [19], where was already demonstrated that

### (c) The Fourier's theory heat equation solution as a limit to the relativistic heat equation solution

It was therefore important to study the behaviour of the solution versus the values of *C*. In this respect, we will now see that *C* and the FHE solution is the limit to the RHE solution. This is equivalent so saying that

To prove this fact, we use the asymptotic expansion

given in Watson [21]. Then we have
*λ* there are no poles of the function

In consequence

### (d) Comparison of the solutions

Figure 1 shows the dimensionless temperature profiles during RF heating of biological tissue along the radial axis for two different dimensionless times (*ξ*=1 and *ξ*=7) obtained from the FHE, HHE and RHE equations. Figure 2 shows the dimensionless temperature evolution during RF heating of biological tissue at two dimensionless locations: *ρ*=1.5 and *ρ*=2.5. The most important result was the evolution of the values obtained for each solution: RHE temperature values were always lower than those of the FHE, while HHE values were higher than the FHE except in the early stages of heating and at points away from the electrode (see for instance temperature at *ξ*=1 for *ρ*>2.2 in figure 1, or temperature at *ρ*=2.5 for *ξ*<1.5 in figure 2). The three solutions converged for infinite time or infinite distance from the heat source. The RHE solution converged towards FHE quicker than the HHE.

Figure 2 shows the peaks in the temperature profile obtained from the HHE, which were due to the hyperbolic nature of (3.13), which includes an essential part of the irregularities of the source term in its solutions (i.e. presence of the Heaviside and Dirac distributions). Such peaks are also detected in numerical treatment of these problems, but they are difficult to represent in numerical plots where they appear smoothed, and this fact is reported in [22]. By contrast, although (3.1) is also hyperbolic, the source term is much more regular (without Heaviside and Dirac distributions), and hence the solution to the RHE equation does not have these peaks (figure 2).

Figure 3 shows the temperature profiles and temperature evolution in dimensional values as an indication of the absolute differences between the three equations. These profiles were obtained for a thermal relaxation time of 16 s in the case of HHE and a value for the speed of heat propagation of *τ* and *C* on the differences in the temperature profiles, in figure 4 we plotted the percentage differences between the solutions from HHE and FHE, and between RHE and FHE. These differences were higher at the start of heating and for locations close to the electrode surface. The differences were seen to increase for higher values of thermal relaxation time in HHE or equivalently for lower values of the speed of heat propagation in RHE (note the decrease in the maximum percentage difference from 6 to 3% in the case of HHE and from −7 to −4.5% in the case of RHE when *τ* was reduced from 16 to 8 s).

## 4. Discussion

The objective of this study was to compare the temperature profiles in biological tissues during RF ablation computed from the Fourier-based, hyperbolic and relativistic heat conduction models. As far as we know, this is the first time these three models have been compared. The relativistic approach was originally described by Ali & Zhang [12] as an alternative heat conduction model, since the FHE assumes an infinite heat propagation speed, and also because there are grounds for believing that the hyperbolic heat conduction model may violate a more general non-equilibrium second law of thermodynamics. Our study did not try to elucidate which model is true, but to qualitatively describe how each approach works for the mathematical modelling of heating of biological tissue, and especially to identify the differences in the computed temperature profiles. In fact, Ali & Zhang [12] did not compare the thermal performance of the relativistic model with the others models. In this respect, our results showed that RHE temperature values were always lower than those of the FHE, where in general the computed temperature in the HHE is higher than the FHE, except at the early stages of heating and at points away from the electrode. Note that if there is no heat sources and the temperature at the boundary would be constant, the relativistic solution would coincide with the hyperbolic one taking *C*^{2}=*α*/*τ*.

Although most experimentally validated theoretical models consider FHE [23,24], some studies have suggested that HHE results match experimental results better than FHE when modelling short pulse laser irradiation [25]. However, the temperature peaks predicted by HHE for some specific points and times have not been usually reported in experimental studies, probably because it is difficult that the those peaks are taken in the time and spatial sampling process. Moreover, it could be also difficult to observe those temperature peaks in numerical studies due to the numerical oscillations around of sharp discontinuities. Finally, it is important to note that the RHE temperature profiles did not show these irregularities.

This study has a limitation related to the simplicity of the model geometry. We considered a spherical electrode totally surrounded by infinite and homogeneous tissue, with thermal and electrical characteristics constant with temperature. All this allowed a one-dimensional scenario, and hence an analytical solution for the three models studied. However, we should point out that our goal was not to obtain accurate and realistic biological tissue heating models from the three equations (which can be achieved with numerical methods), but to conduct a preliminary assessment of the thermal performance of each thermal model. In spite of the power and accuracy of current numerical methods to realistically model the EBATs, we consider it essential to first use a comparative approach with an analytical solution, as reported here.

Our modelling study did not include the blood perfusion term, which is one of the most important characteristic of the bioheat equation. In spite of this, this term acts removing heat such as has been observed in computer simulations [26] and hence implies that the temperatures profiles are lower than in the case without perfusion. Consequently, we think that our comparison among models would not change in presence of the blood perfusion term.

As regards the most suitable heat model for RFA, in the light of our results future experimental and theoretical studies should be conducted that could consider not only FHE and HHE, but also RHE. The present results also confirm that the non-Fourier heat equations are suitable for modelling temperature for short-time heating close to the heat source, since that the solutions of all the models converged for infinite time or infinite distance from the electrode.

There is still a controversy around the value of thermal relaxation time in biological tissues. Although Mitra *et al.* [13] measured values around 16 s in processed meat, Ali & Zhang's [12] reported a smaller value. Likewise, Kaminski [7] reported a relaxation time between 20 s and 30 s for meat products and Vedavarz *et al.* [27] found an interval between 1 s and 100 s for some biological tissues at room temperature. In this respect, the relativistic heat theory considers a similar parameter *C*, known as speed of heat propagation, which is somehow inversely proportional to *τ*. Ali & Zhang [12] found experimentally that the speed of heat (second sound) is about one order of magnitude less than the speed of sound. If we consider a value of the sound propagation in biological tissues around 1500 m s^{−1}, this would imply a value of *C* around 150 m s^{−1}, which for a value of *α* around 1.5×10^{−7} m^{2} s^{−1}, would estimate *τ* as *α*/*C*^{2}, i.e. a value around of 6.7×10^{−12} s. In practical terms, this means that the relativistic and hyperbolic heat models would only be useful for very short heating times in biological tissues. However, Ali & Zhang [12] also point out that these relations are for metals and perfect structures. On the other hand, non-homogeneous materials, such as biological tissues, have very long relaxation times, and consequently would have very low heating speeds [13,7]. Unfortunately, no experimental data have been published to date on the speed of heat propagation in biological tissues.

## 5. Conclusion

Our findings suggest that temperature values obtained from a model with RHE are always lower than those of the FHE, while HHE values are higher than the FHE, except for the initial heating stage and at points away from the electrode. Both HHE and RHE are mathematically hyperbolic, but temperature profile peaks are only observed with HHE. The three solutions converged for infinite time or infinite distance from the electrode. The percentage differences between the FHE and the other solutions are larger for higher values of thermal relaxation time in HHE or equivalently for lower values of the speed of heat propagation in RHE.

## Funding statement

This work received financial support from the Spanish Government (Ministerio de Ciencia e Innovación, Ref. TEC2011-27133-C02-01).

## Appendix A. Computation of the inverse Laplace transform in the relativistic heat equation

**(a) Computation of **

By the translation theorem for Laplace transforms and formula 36, chapter 5, section 5.6 from Erdélyi *et al.* [28], we have

**(b) Computation of **

In an analogous way, we obtain

Remark that
*x*≤0 and *x*>0.

**(c) Computation of **

Computation of

As above we have

On the other hand,

The last inverse can be computed with help of a Bromwich's contour. There are two branch points *s*=0 and *s*=−1/*λ* and two possible poles verifying
*s*_{1}<−1/*λ*,*s*_{2}<−1/*λ*, 1+*Bs*_{1}>0 and 1+*Bs*_{2}<0. Then only *s*_{1} is a pole. Using that

Then

Now we have

We obtain

Then by (A 3)–(A 5), we have

- Received July 18, 2014.
- Accepted October 7, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.