Analytical solution of one-dimensional Pennes’ bioheat equation

Hongyun Wang; Wesley A. Burgei; Hong Zhou

doi:10.1515/phys-2020-0197

Article Open Access

Analytical solution of one-dimensional Pennes’ bioheat equation

Hongyun Wang , Wesley A. Burgei and Hong Zhou

Published/Copyright: December 29, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 18 Issue 1

Abstract

Pennes’ bioheat equation is the most widely used thermal model for studying heat transfer in biological systems exposed to radiofrequency energy. In their article, “Effect of Surface Cooling and Blood Flow on the Microwave Heating of Tissue,” Foster et al. published an analytical solution to the one-dimensional (1-D) problem, obtained using the Fourier transform. However, their article did not offer any details of the derivation. In this work, we revisit the 1-D problem and provide a comprehensive mathematical derivation of an analytical solution. Our result corrects an error in Foster’s solution which might be a typo in their article. Unlike Foster et al., we integrate the partial differential equation directly. The expression of solution has several apparent singularities for certain parameter values where the physical problem is not expected to be singular. We show that all these singularities are removable, and we derive alternative non-singular formulas. Finally, we extend our analysis to write out an analytical solution of the 1-D bioheat equation for the case of multiple electromagnetic heating pulses.

Keywords: Pennes’ bioheat equation; blood flow cooling; surface cooling; analytical solution; pulsed heating

1 Introduction

Heating of biological systems by radiofrequency (RF) energy has many applications, and it has been studied intensively in the literature (see ref. [1] and references therein). The bioheat equation attributed to Pennes is the most commonly used formulation of heat transfer in biological systems such as human skin [2].

In 1978, Foster et al. applied the Fourier transform to Pennes’ bioheat equation and obtained the temperature increase in a semi-infinite domain of tissue irradiated by an incident electromagnetic beam perpendicular to the surface of tissue domain [3]. While an analytical expression of the temperature increase was given as a function of depth and time in ref. [3], no mathematical details regarding derivation were provided. Later, Nelson et al. used the Laplace transform to obtain a formula for temperature at the tissue surface only [4]; they did not give an explicit formula for the temperature profile at locations other than the skin surface. The Laplace transform has also been employed in studying anomalous diffusion equations with the decay exponential kernel [5]. Green’s functions for Pennes’ bioheat equation in various cases were studied in refs. [6,7,8,9]. Exact solutions in the form of infinite series were obtained in the case of one-dimensional (1-D) multiple planar tissue layers or spherical tissue layers, based on the method of separation of variables and Green’s function method [10,11]. However, none of these efforts addressed the fundamental question of how to explicitly obtain the analytical solution of Pennes’ bioheat equation. Therefore, in this work we revisit the 1-D bioheat transfer problem governed by Pennes’ equation. We present a detailed derivation of the analytical solution for RF irradiation of semi-infinite media. The mathematical formulation contains (i) a heating-source term (electromagnetic beam), (ii) a temperature-restoring term (blood flow), and (iii) a boundary condition reflecting surface cooling (radiation and convection). Our approach to the solution relies on transforming this boundary value formulation of non-homogeneous heat equation to an initial value problem of the standard heat equation.

2 Mathematical formulation

Following ref. [3], we consider the problem of microwave heating of skin tissue coupled with the effects of surface cooling and blood flow. Heating the skin by other RF energy such as millimeter waves can be studied in the same formulation with different parameter values. We treat the living tissue as semi-infinite, homogeneous, and isotropic media and thereby only focus on a 1-D case with constant thermal parameters. We first introduce notations for independent variables, functions, and parameters, similar to those used in ref. [3]:

( x ′ , t ′ ) ≡ ( depth, time ) : original independent variables before scaling
( x , t ) ≡ (depth, time) : new independent variables after scaling
T 0 : temperature of the arterial blood entering the tissue, independent of ( x ′ , t ′ )
T e : temperature of the environment outside the skin, independent of ( x ′ , t ′ )
T ′ ( x ′ , t ′ ) : temperature of the tissue as a function of ( x ′ , t ′ )
T ( x ′ , t ′ ) ≡ T ′ ( x ′ , t ′ ) − T 0 : temperature of the tissue relative to T 0 .

Here, the skin surface is at x ′ = 0 (zero depth). The relative temperature T ( x ′ , t ′ ) is governed by an initial boundary value problem (IBVP) of Pennes’ bioheat equation:

(1) μ ∂ T ∂ t ′ = ∂ 2 T ∂ x ′ 2 − λ T + q 0 e − x ′ / L ∂ T ( x ′ , t ′ ) ∂ x ′ x ′ = 0 = α ( T ( 0 , t ′ ) + T 0 − T e ) T ( x ′ , 0 ) = T p ( x ′ ) ,

where μ ≡ ρ C p / k , λ ≡ V s / k , q 0 ≡ P 0 / ( k L ) , and

ρ : mass density of the tissue
C p : specific heat of the tissue
k: heat conductivity of the tissue
V s : product of flow rate, mass density, and specific heat of the blood flow
P 0 : microwave power absorbed into the tissue per unit skin surface area
L: reciprocal of absorption coefficient of the microwave frequency used, describing the depth at which the microwave power deposition is reduced by a factor e.

The effect of surface cooling on temperature evolution is described by coefficient α ≡ h c / k , where h c is the heat transfer coefficient of radiation and convection cooling.

The initial temperature profile of the tissue prior to the start of microwave heating, denoted by T p ( x ′ ) , is the steady state of the skin tissue, corresponding to the balance between the blood flow warming and the surface cooling. Mathematically, T p ( x ′ ) can be obtained through solving a boundary value problem (BVP)

(2) 0 = d 2 T p d x ′ 2 − λ T p d T p ( x ′ ) d x ′ x ′ = 0 = α ( T p ( 0 ) + T 0 − T e ) .

Solving (2) gives us

(3) T p ( x ′ ) = α ( T e − T 0 ) α + λ e − x ′ λ .

To separate out the effect of initial temperature profile T p ( x ′ ) in (1), we introduce

(4) T h ( x ′ , t ′ ) ≡ T ( x ′ , t ′ ) − T p ( x ′ ) .

T h ( x ′ , t ′ ) is the temperature increase attributed to microwave heating. It satisfies

(5) μ ∂ T h ∂ t ′ = ∂ 2 T h ∂ x ′ 2 − λ T h + q 0 e − x ′ / L ∂ T h ( x ′ , t ′ ) ∂ x ′ x ′ = 0 = α T h ( 0 , t ′ ) T h ( x ′ , 0 ) = 0 .

IBVP (5) constitutes the mathematical formulation for our analysis.

3 Transforming the formulation to an initial value problem of the standard heat equation

IBVP (5) contains a non-insulated boundary condition simulating the surface cooling, a source of microwave heating, and a temperature restoring term caused by blood flow. The objective of this section is to transform (5) to an initial value problem (IVP) of the standard heat equation, which allows us to write out an analytical solution. We start by scaling ( x ′ , t ′ , T h ) to simplify the coefficients in (5). We define new variables, temperature function, and parameters as follows:

(6) x ≡ x ′ / L t ≡ t ′ / ( μ L 2 ) T s ( x , t ) ≡ T h ( x ′ , t ′ ) / ( q 0 L 2 ) λ s ≡ λ L 2 α s ≡ α L .

The scaled temperature T s ( x , t ) as a function of scaled independent variables satisfies

(7) ∂ T s ∂ t = ∂ 2 T s ∂ x 2 − λ s T s + e − x ∂ T s ( x , t ) ∂ x x = 0 = α s T s ( 0 , t ) T s ( x , 0 ) = 0 .

Once T s ( x , t ) is solved, T ( x ′ , t ′ ) can be expressed in the form of

(8) T ( x ′ , t ′ ) = ( q 0 L 2 ) T s x ′ L , t ′ μ L 2 + α ( T e − T 0 ) α + λ e − x ′ λ .

Our goal is to transform (7) to an IVP of the standard heat equation. On the road to that destination, our first step is to remove the inhomogeneous term in the partial differential equation (PDE). We apply Duhamel’s principle to move the effect of the heating into the initial value. Duhamel’s principle provides a convenient approach for obtaining a solution to the inhomogeneous heat equation by utilizing the solution to the auxiliary homogeneous problem [12]. We consider the auxiliary IBVP of the homogeneous PDE corresponding to (7):

(9) ∂ G ∂ t = ∂ 2 G ∂ x 2 − λ s G ∂ G ( x , t ) ∂ x x = 0 = α s G ( 0 , t ) G ( x , 0 ) = e − x .

Duhamel’s principle states that solutions of (7) and (9) are related by the equation [12]

(10) T s ( x , t ) = ∫ 0 t G ( x , t − s ) d s .

Thus, the task becomes solving (9). Our next step is to get rid of the temperature restoring term ( − λ s G ) . We write function G ( x , t ) as

(11) G ( x , t ) ≡ e − λ s t u ( x , t ) .

It follows from (9) that u ( x , t ) satisfies

(12) ∂ u ∂ t = ∂ 2 u ∂ x 2 ∂ u ( x , t ) ∂ x x = 0 = α s u ( 0 , t ) u ( x , 0 ) = e − x .

The remaining obstacle in (12) is the surface cooling boundary condition. We write it mathematically as an insulated boundary condition on a transformed function. We map function u ( x , t ) to v ( x , t ) as follows:

(13) v ( x , t ) = 1 1 + α s u ( x , t ) + α s 1 + α s ∫ x + ∞ u ( ξ , t ) d ξ .

It is straightforward to verify that v ( x , t ) is governed by

(14) ∂ v ∂ t = ∂ 2 v ∂ x 2 ∂ v ( x , t ) ∂ x x = 0 = 0 v ( x , 0 ) = e − x .

The insulated boundary condition for v ( x , t ) allows us to extend it as an even function of x. After the even extension, (14) is converted to an IVP of the standard heat equation:

(15) ∂ v ∂ t = ∂ 2 v ∂ x 2 v ( x , 0 ) = e − | x | .

The solution of (15) can be written in terms of the Gaussian heat kernel as [12]

v ( x , t ) = ∫ − ∞ + ∞ 1 4 π t e − ( ξ − x ) 2 4 t v ( ξ , 0 ) d ξ = 1 4 π t ∫ 0 + ∞ e − ( ξ + x ) 2 4 t + e − ( ξ − x ) 2 4 t e − ξ d ξ = 1 4 π t ∫ 0 + ∞ e − ( ξ 2 + 2 ( 2 t + x ) ξ + x 2 ) 4 t + e − ( ξ 2 + 2 ( 2 t − x ) ξ + x 2 ) 4 t d ξ .

Completing the square in the exponent of each integral, we write v ( x , t ) as

(16) v ( x , t ) = 1 4 π t e t + x ∫ 0 + ∞ e − ( ξ + ( 2 t + x ) ) 2 4 t d ξ + e t − x ∫ 0 + ∞ e − ( ξ + ( 2 t − x ) ) 2 4 t d ξ = e t + x 2 erfc 2 t + x 4 t + e t − x 2 erfc 2 t − x 4 t ,

where erfc ( x ) ≡ 2 π ∫ x + ∞ e − t 2 d t is the complementary error function. In the first integral of (16), we have used change of variables s = ( ξ + ( 2 t + x ) ) 4 t :

∫ 0 + ∞ e − ( ξ + ( 2 t + x ) ) 2 4 t d ξ = 4 t ∫ 2 t + x 4 t + ∞ e − s 2 d s = 4 t π 2 erfc 2 t + x 4 t .

The second integral in (16) is treated in a similar way.

4 Analytical solution

Starting from the analytical solution v ( x , t ) given in (16), we go step by step, to u ( x , t ) , to G ( x , t ) , to T s ( x , t ) , and finally write out T ( x ′ , t ′ ) in original variables ( x ′ , t ′ ) .

4.1 Solution of u ( x , t )

We invert mapping (13) to express u ( x , t ) in terms of v ( x , t ) . Let

(17) w ( x , t ) ≡ ∫ x + ∞ u ( ξ , t ) d ξ .

From equation (13), we get a linear differential equation on w ( x , t ) :

− ∂ w ( x , t ) ∂ x + α s w ( x , t ) = ( 1 + α s ) v ( x , t ) .

We multiply it by the integrating factor e − α s x and then integrate it to solve for w ( x , t ) . This process yields

− ∂ ( e − α s x w ( x , t ) ) ∂ x = ( 1 + α s ) e − α s x v ( x , t ) e − α s x w ( x , t ) = ( 1 + α s ) ∫ x + ∞ e − α s ξ v ( ξ , t ) d ξ w ( x , t ) = ( 1 + α s ) e α s x ∫ x + ∞ e − α s ξ v ( ξ , t ) d ξ .

Differentiating w ( x , t ) with respect to x gives u ( x , t ) in terms of v ( x , t ) :

(18) u ( x , t ) = − ∂ w ( x , t ) ∂ x = ( 1 + α s ) v ( x , t ) − α s e α s x ∫ x + ∞ e − α s ξ v ( ξ , t ) d ξ .

Before we use (18) to derive an expression of u ( x , t ) , we work out a general integration formula for A ≡ ∫ x + ∞ e − c 2 ξ erfc ( c 1 ( ξ + c 0 ) ) d ξ in the case of c 1 > 0 or c 2 > 0 or both. We carry out integration by parts and then complete the square in the exponent to bring the integral into the form:

A ≡ ∫ x + ∞ e − c 2 ξ erfc ( c 1 ( ξ + c 0 ) ) d ξ = − 1 c 2 ∫ x + ∞ erfc ( c 1 ( ξ + c 0 ) ) d ( e − c 2 ξ ) = e − c 2 x c 2 erfc ( c 1 ( x + c 0 ) ) − 2 c 1 c 2 π ∫ x + ∞ e − c 1 2 ( ξ + c 0 ) 2 e − c 2 ξ d ξ = e − c 2 x c 2 erfc ( c 1 ( x + c 0 ) ) − 2 c 1 c 2 π e c 2 2 / ( 4 c 1 2 ) + c 0 c 2 ∫ x + ∞ e − c 1 2 ( ξ + c 0 + c 2 / ( 2 c 1 2 ) ) 2 d ξ .

Expressing the second part in terms of erfc() and noting c 1 2 = | c 1 | , we obtain

(19) A = e − c 2 x c 2 erfc ( c 1 ( x + c 0 ) ) − sign ( c 1 ) e c 2 2 / ( 4 c 1 2 ) + c 0 c 2 c 2 erfc | c 1 | x + c 0 + c 2 2 c 1 2 .

The condition of c 1 > 0 or c 2 > 0 is necessary for the convergence of the improper integral A . When both c 1 and c 2 are negative, the integrand does not decrease to zero as ξ → + ∞ .

Formula (19) appears to have a singularity as c 2 → 0 . However, for c 1 > 0 , the singularity is removable. When c 1 > 0 , we rewrite (19) as

(20) A = 1 c 2 e − z x z = 0 z = c 2 erfc ( c 1 ( x + c 0 ) ) − 1 c 2 e z 2 / ( 4 c 1 2 ) + c 0 z erfc c 1 x + c 0 + z 2 c 1 2 z = 0 z = c 2 .

In the limit of c 2 → 0 , both terms in (20) converge to derivatives of the corresponding functions with respect to z at z = 0 . Thus, we have

lim c 2 → 0 A = − ( x + c 0 ) ⋅ erfc ( c 1 ( x + c 0 ) ) + 1 c 1 π e − c 1 2 ( x + c 0 ) 2 .

With the help of formula (19), we evaluate the integral in (18) by substituting expression (16) of v ( x , t ) into the integral:

∫ x + ∞ e − α s ξ v ( ξ , t ) d ξ = ∫ x + ∞ e − α s ξ e t + ξ 2 erfc 2 t + ξ 4 t + e t − ξ 2 erfc 2 t − ξ 4 t d ξ = e t 2 ∫ x + ∞ e − ( α s − 1 ) ξ erfc 2 t + ξ 4 t + e − ( α s + 1 ) ξ erfc 2 t − ξ 4 t d ξ .

We apply formula (19) on each part inside the brackets. The first part has c 1 = 1 / 4 t , c 0 = 2 t , and c 2 = ( α s − 1 ) ; the second has c 1 = − 1 / 4 t , c 0 = − 2 t , and c 2 = ( α s + 1 ) . The integral becomes

(21) ∫ x + ∞ e − α s ξ v ( ξ , t ) d ξ = e t 2 e − ( α s − 1 ) x α s − 1 erfc 2 t + x 4 t + e − ( α s + 1 ) x α s + 1 erfc 2 t − x 4 t + 2 e − ( 1 − α s 2 ) t 1 − α s 2 erfc 2 α s t + x 4 t .

Putting expression (21) into equation (18), we find the expression of u ( x , t ) :

(22) u ( x , t ) = 1 + α s 2 ( 1 − α s ) e t + x erfc 2 t + x 4 t + 1 2 e t − x erfc 2 t − x 4 t − α s 1 − α s e α s x + α s 2 t erfc 2 α s t + x 4 t .

4.2 Several terms in the calculation of T s ( x , t )

Equations (10) and (11) give T s ( x , t ) in terms of u ( x , t ) :

(23) T s ( x , t ) = ∫ 0 t e − λ s ( t − s ) u ( x , t − s ) d s = ∫ 0 t e − λ s τ u ( x , τ ) d τ .

We now use u ( x , t ) obtained in (22) to calculate T s ( x , t ) .

(24) T s ( x , t ) = ∫ 0 t e − λ s τ u ( x , τ ) d τ = ( 1 + α s ) 2 ( 1 − α s ) e x ∫ 0 t e ( 1 − λ s ) τ erfc 2 τ + x 4 τ d τ + 1 2 e − x ∫ 0 t e ( 1 − λ s ) τ erfc 2 τ − x 4 τ d τ − α s 1 − α s e α s x ∫ 0 t e ( α s 2 − λ s ) τ erfc 2 α s τ + x 4 τ d τ ≡ ( 1 + α s ) 2 ( 1 − α s ) I 1 + 1 2 I 2 − α s 1 − α s I 3 .

We perform integration by parts on each term in (24). For I 1 , we have

(25) I 1 = e x ∫ 0 t erfc 2 τ + x 4 τ d e ( 1 − λ s ) τ 1 − λ s = e x e ( 1 − λ s ) t 1 − λ s erfc 2 t + x 4 t + e x ∫ 0 t e ( 1 − λ s ) τ 1 − λ s 2 π e − ( 2 τ + x ) 2 4 τ d 2 τ + x 4 τ ≡ I 1 A + I 1 B 1 − λ s .

In I 1 B , we combine all exponents and then complete the square, which gives

(26) I 1 B = 2 π ∫ 0 t e − λ s τ − x 2 4 τ d 2 τ + x 4 τ = e λ s x 2 π ∫ 0 t e − 2 λ s τ + x 2 τ 2 1 + λ s 2 λ s d 2 λ s τ + x 2 τ + 1 − λ s 2 λ s d 2 λ s τ − x 2 τ = 1 + λ s 2 λ s e λ s x 2 π ∫ 0 t e − x + 2 λ s τ 2 τ 2 d x + 2 λ s τ 2 τ − 1 − λ s 2 λ s e − λ s x 2 π ∫ 0 t e − x − 2 λ s τ 2 τ 2 d x − 2 λ s τ 2 τ = − 1 + λ s 2 λ s e λ s x erfc x + 2 λ s t 2 t + 1 − λ s 2 λ s e − λ s x erfc x − 2 λ s t 2 t .

In a similar manner, I 2 is calculated. The noticeable difference is that in the integration by parts of I 2 , at t = 0 + , we have erfc ( − ∞ ) = 2 , instead of erfc ( + ∞ ) = 0 .

(27) I 2 = e − x ∫ 0 t erfc 2 τ − x 4 τ d e ( 1 − λ s ) τ 1 − λ s = e − x e ( 1 − λ s ) t 1 − λ s erfc 2 t − x 4 t − e − x 2 1 − λ s + e − x ∫ 0 t e ( 1 − λ s ) τ 1 − λ s 2 π e − ( 2 τ − x ) 2 4 τ d 2 τ − x 4 τ ≡ I 2 A + I 2 B 1 − λ s .

In I 2 B , we collect all exponents, complete the square and rearrange the terms, which yields

(28) I 2 B = 2 π ∫ 0 t e − λ s τ − x 2 4 τ d 2 τ − x 4 τ = e λ s x 2 π ∫ 0 t e − 2 λ s τ + x 2 τ 2 × 1 + λ s 2 λ s d 2 λ s τ − x 2 τ + 1 − λ s 2 λ s d 2 λ s τ + x 2 τ = 1 − λ s 2 λ s e λ s x 2 π ∫ 0 t e − x + 2 λ s τ 2 τ 2 d x + 2 λ s τ 2 τ − 1 + λ s 2 λ s e − λ s x 2 π ∫ 0 t e − x − 2 λ s τ 2 τ 2 d x − 2 λ s τ 2 τ = − 1 − λ s 2 λ s e λ s x erfc x + 2 λ s t 2 t + 1 + λ s 2 λ s e − λ s x erfc x − 2 λ s t 2 t .

For I 3 , the calculation is again similar to that of I 1 .

(29) I 3 = e α s x ∫ 0 t erfc 2 α s τ + x 4 τ d e ( α s 2 − λ s ) τ α s 2 − λ s = e α s x e ( α s 2 − λ s ) t α s 2 − λ s erfc 2 α s t + x 4 t + e α s x ∫ 0 t e ( α s 2 − λ s ) τ α s 2 − λ s 2 π e − ( 2 α s τ + x ) 2 4 τ d 2 α s τ + x 4 τ ≡ I 3 A + I 3 B α s 2 − λ s .

In I 3 B , we gather all exponents and complete the square, which produces

(30) I 3 B = 2 π ∫ 0 t e − λ s τ − x 2 4 τ d 2 α s τ + x 4 τ = e λ s x 2 π ∫ 0 t e − 2 λ s τ + x 2 τ 2 × α s + λ s 2 λ s d 2 λ s τ + x 2 τ + α s − λ s 2 λ s d 2 λ s τ − x 2 τ = α s + λ s 2 λ s e λ s x 2 π ∫ 0 t e − x + 2 λ s τ 2 τ 2 d x + 2 λ s τ 2 τ − α s − λ s 2 λ s e − λ s x 2 π ∫ 0 t e − x − 2 λ s τ 2 τ 2 d x − 2 λ s τ 2 τ = − α s + λ s 2 λ s e λ s x erfc x + 2 λ s t 2 t + α s − λ s 2 λ s e − λ s x erfc x − 2 λ s t 2 t .

4.3 Solution of T s ( x , t )

Incorporating results (25)–(30) into equation (24), we obtain T s ( x , t ) in the form:

(31) T s ( x , t ) = ( 1 + α s ) 2 ( 1 − α s ) e x + ( 1 − λ s ) t 1 − λ s erfc 2 t + x 2 t − ( 1 + α s ) 2 ( 1 − α s ) ( 1 + λ s ) e λ s x 2 ( 1 − λ s ) λ s erfc x + 2 λ s t 2 t + ( 1 + α s ) 2 ( 1 − α s ) ( 1 − λ s ) e − λ s x 2 ( 1 − λ s ) λ s erfc x − 2 λ s t 2 t + 1 2 ⋅ e − x + ( 1 − λ s ) t 1 − λ s erfc 2 t − x 2 t − e − x 1 − λ s − 1 2 ⋅ ( 1 − λ s ) e λ s x 2 ( 1 − λ s ) λ s erfc x + 2 λ s t 2 t + 1 2 ⋅ ( 1 + λ s ) e − λ s x 2 ( 1 − λ s ) λ s erfc x − 2 λ s t 2 t − α s 1 − α s e α s x + ( α s 2 − λ s ) t α s 2 − λ s erfc 2 α s t + x 2 t + α s 1 − α s ( α s + λ s ) e λ s x 2 ( α s 2 − λ s ) λ s erfc x + 2 λ s t 2 t − α s 1 − α s ( α s − λ s ) e − λ s x 2 ( α s 2 − λ s ) λ s erfc x − 2 λ s t 2 t .

Combining like terms in (31), we express the final expression for T s ( x , t ) as

(32) T s ( x , t ) = 1 λ s − 1 e − x + 1 + α s 2 λ s − α s ⋅ e λ s x erfc x + 2 λ s t 2 t − 1 + α s 2 ( λ s + α s ) ⋅ e − λ s x erfc x − 2 λ s t 2 t − 1 2 e − x + ( 1 − λ s ) t erfc 2 t − x 2 t − ( 1 + α s ) 2 ( 1 − α s ) e x + ( 1 − λ s ) t erfc 2 t + x 2 t + α s ( λ s − 1 ) ( 1 − α s ) ( λ s − α s 2 ) e α s x + ( α s 2 − λ s ) t erfc 2 α s t + x 2 t .

4.4 Analytical solution in original variables ( x ′ , t ′ ) and comparison with Foster’s result

Based on the expression of T s ( x , t ) we just derived in (32) and the scaling relations given in (6) and (8), we arrive at the following formula for T ( x ′ , t ′ ) in terms of the original variables ( x ′ , t ′ ) :

(33) T ( x ′ , t ′ ) = q 0 λ − 1 / L 2 e − x ′ / L + 1 / L + α 2 ( λ − α ) ⋅ e λ x ′ erfc x ′ 2 μ t ′ + λ t ′ μ − 1 / L + α 2 ( λ + α ) ⋅ e − λ x ′ erfc x ′ 2 μ t ′ − λ t ′ μ − 1 2 e − x ′ L + 1 L 2 − λ t ′ μ erfc 1 L t ′ μ − x ′ 2 μ t ′ − ( 1 / L + α ) 2 ( 1 / L − α ) e x ′ L + 1 L 2 − λ t ′ μ erfc 1 L t ′ μ + x ′ 2 μ t ′ + α ( λ − 1 / L 2 ) ( 1 / L − α ) ( λ − α 2 ) e α x ′ + ( α 2 − λ ) t ′ μ erfc α t ′ μ + x ′ 2 μ t ′ + α ( T e − T 0 ) α + λ e − x ′ λ .

The only difference between Foster’s result and ours is the factor ( λ − 1 / L 2 ) in the second to last term in (33). In ref. [3], the corresponding factor was stated as ( λ + 1 / L 2 ) .

5 Removable singularities in (32)

Expression (32) appears to have singularities when α s = 1 or λ s = 1 or λ s = α s 2 . These singularities are removable. We discuss these cases one by one.

5.1 Case 1: λ s = 1 and α s ≠ 1

We calculate lim λ s → 1 T s ( x , t ) . Let ε ≡ λ s − 1 . We write λ s = 1 + ε and write equation (32) in terms of ε

(34) T s ( x , t ) = 1 ε ( 2 + ε ) e − x + 1 + α s 2 ( 1 − α s + ε ) ⋅ e x + ε x erfc x + 2 t + 2 ε t 2 t − 1 + α s 2 ( 1 + α s + ε ) ⋅ e − x − ε x erfc x − 2 t − 2 ε t 2 t − 1 2 e − x − ε ( 2 + ε ) t erfc 2 t − x 2 t − ( 1 + α s ) 2 ( 1 − α s ) e x − ε ( 2 + ε ) t erfc 2 t + x 2 t + α s ε ( 2 + ε ) ( 1 − α s ) ( 1 − α s 2 + ε ( 2 + ε ) ) × e α s x + ( α s 2 − 1 ) t − ε ( 2 + ε ) t erfc 2 α s t + x 2 t .

In the limit as ε → 0 , we have

(35) T s ( x , t ) λ s = 1 ≡ lim λ s → 1 T s ( x , t ) = 1 + α s 4 ( 1 − α s ) x + 2 t − 1 1 − α s e x erfc x + 2 t 2 t + 1 4 x − 2 t + 1 1 + α s e − x erfc x − 2 t 2 t + t e − x − 1 ( 1 − α s ) π t e − x 2 4 t + t + α s ( 1 − α s ) ( 1 − α s 2 ) e α s x + ( α s 2 − 1 ) t erfc 2 α s t + x 2 t .

5.2 Case 2: λ s = 1 and α s = 1

In order to evaluate lim α s → 1 ( lim λ s → 1 T s ( x , t ) ) , we write equation (35) in terms of η ≡ α s − 1

(36) lim λ s → 1 T s ( x , t ) = 2 + η − 4 η x + 2 t + 1 η e x erfc x + 2 t 2 t + 1 4 x − 2 t + 1 2 + η e − x erfc x − 2 t 2 t + t e − x + 1 η π t e − x 2 4 t + t + 1 + η η 2 ( 2 + η ) e x + η x + η ( 2 + η ) t erfc 2 t + 2 η t + x 2 t .

Taking the limit as η → 0 yields

(37) T s ( x , t ) | λ s = 1 , α s = 1 ≡ lim α s → 1 ( lim λ s → 1 T s ( x , t ) = 1 4 ( x + 2 t ) 2 + 2 t − 1 2 e x erfc x + 2 t 2 t + 1 4 x − 2 t + 1 2 e − x erfc x − 2 t 2 t + t e − x − t 2 π ( x + 2 t + 1 ) e − x 2 4 t + t .

5.3 Case 3: λ s = α s 2 with α s ≠ 1 and α s > 0

To determine lim λ s → α s 2 T s ( x , t ) , we introduce ε ≡ λ s − α s and substitute λ s = α s + ε into equation (32) to obtain

(38) T s ( x , t ) = 1 α s 2 − 1 + O ( ε ) × e − x + 1 + α s 2 ε ⋅ e α s x + ε x erfc x + 2 α s t + 2 ε t 2 t − 1 + α s 4 α s e − α s x erfc x − 2 α s t 2 t − e ( 1 − α s 2 ) t × 1 2 e − x erfc 2 t − x 2 t + ( 1 + α s ) 2 ( 1 − α s ) e x erfc 2 t + x 2 t + α s ( α s 2 − 1 + ε ( 2 α s + ε ) ) ( 1 − α s ) ε ( 2 α s + ε ) × e α s x − ε ( 2 α s + ε ) t erfc 2 α s t + x 2 t + O ( ε ) .

As ε approaches zero, we find that

(39) T s ( x , t ) λ s = α s 2 ≡ lim λ s → α s 2 T s ( x , t ) = 1 α s − 1 e − x α s + 1 − t π e − x 2 4 t + α s 2 t − 1 4 α s e − α s x erfc x − 2 α s t 2 t − e ( 1 − α s 2 ) t × e − x 2 ( 1 + α s ) erfc 2 t − x 2 t + e x 2 ( 1 − α s ) erfc 2 t + x 2 t + 1 2 x + 2 α s t − 3 α s 2 + 1 2 α s ( α s 2 − 1 ) e α s x erfc x + 2 α s t 2 t .

5.4 Case 4: λ s = 0 and α s = 0

In computing lim λ s → 0 ( lim α s → 0 T s ( x , t ) ) , we calculate lim α s → 0 T s ( x , t ) and write it in terms of ε ≡ λ s

(40) lim α s → 0 T s ( x , t ) = 1 ε 2 − 1 e − x + 1 2 ε e ε x erfc x + 2 ε t 2 t − 1 2 ε e − ε x erfc x − 2 ε t 2 t − 1 2 e t − x e r f c 2 t − x 2 t − 1 2 e x + t erfc 2 t + x 2 t + O ( ε 2 ) .

Letting ε → 0 , we see that

(41) T s ( x , t ) | λ s = 0 , α s = 0 ≡ lim λ s → 0 ( lim α s → 0 T s ( x , t ) = 1 2 e t − x erfc 2 t − x 2 t + 1 2 e x + t erfc 2 t + x 2 t − e − x − x ⋅ erfc x 2 t + 2 t π e − x 2 4 t .

6 Solution for a sequence of heating pulses

In formulation (5), the microwave power is on and remains constant after t > 0 . Mathematically, after scaling transformation (6), the heating term q ( x , t ) in IBVP (7) is expressed in the form:

q ( x , t ) = H ( t ) e − x ,

where H ( t ) is the Heaviside step function ( H ( t ) = 1 for t > 0 ). In the case of constant power heating, the analytical solution of T s ( x , t ) for t > 0 is derived in Section 4 and is given explicitly in (32). To facilitate the discussion below, we extend function T s ( x , t ) to t ∈ ( − ∞ , + ∞ ) as follows:

(42) T s ( x , t ) = 0 , for t ≤ 0 given by ( 32 ) , for t > 0 .

We consider the case where the microwave power is a sequence of rectangular pulses. Each heating pulse is specified by a starting time t j ( s ) and an ending time t j ( e ) . For a sequence of M pulses, the heating term q ( x , t ) in IBVP (7) takes the form:

(43) q ( x , t ) = ∑ j = 1 M ( H ( t j ( s ) ) − H ( t j ( e ) ) ) e − x .

Let T s ( x , t ; { t j ( s ) , t j ( s ) , 1 ≤ j ≤ M } ) denote the temperature increase caused by the sequence of microwave pulses. The analytical solution of T s ( x , t ; { t j ( s ) , t j ( s ) } ) is given by a superposition of T s ( x , t )

(44) T s ( x , t ; { t j ( s ) , t j ( s ) , 1 ≤ j ≤ M } ) = ∑ j = 1 M ( T s ( x , t − t j ( s ) ) − T s ( x , t − t j ( e ) ) ) .

7 Conclusion

The bioheat equation attributed to Pennes [2] is applied most often to study RF irradiation heating of biological systems such as human skin. In ref. [3], Foster et al. analyzed the 1-D Pennes equation describing the heating of a semi-infinite domain of skin tissue irradiated by an incident microwave beam. They used the Fourier transform to obtain an analytical expression of the resulting temperature increase as a function of time and depth from the skin surface. However, no detail regarding the mathematical derivation was given in their paper. Here, we revisited the 1-D bioheat transfer problem to present a detailed mathematical derivation of the analytical solution based on integrating the PDE directly. Our analytical expression corrected an error in Foster’s result which might be a typo [3]. The analytical expression has apparent singularities for several combinations of parameter values where the physical problem is not expected to have any singularity. We found that all these apparent singularities are removable and derived an alternative non-singular formula for each case. Finally, we extended the analytical solution to the case where the skin tissue is heated by a sequence of RF pulses each specified by a starting time and an end time.

Acknowledgments

The authors acknowledge the Joint Intermediate Force Capabilities Office of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the US Government.

Conflict of interest: The authors declare that there is no conflict of interest.

References

[1] Foster KR, Ziskin MC, Balzano Q. Thermal response of human skin to microwave energy: a critical review. Health Phys. 2016;111(6):528–41.10.1097/HP.0000000000000571Search in Google Scholar PubMed

[2] Pennes H. Analysis of tissue and arterial blood temperatures in the resting human forearm. J Appl Physiol. 1948;1:93–122.10.1152/jappl.1948.1.2.93Search in Google Scholar PubMed

[3] Foster KR, Kritikos HN, Schwan HP. Effect of surface cooling and blood flow on the microwave heating of tissue. IEEE Trans Biomed Eng. 1978;25(3):313–6.10.1109/TBME.1978.326313Search in Google Scholar PubMed

[4] Nelson DA, Walters TJ, Ryan KL, Emerton KB, Hurt WD, Ziriax JM, et al. Inter-species extrapolation of skin heating resulting from millimeter wave irradiation: modeling and experimental results. Health Phys. 2003;84:608–15.10.1097/00004032-200305000-00006Search in Google Scholar PubMed

[5] Yang X-J, Feng Y-Y, Cattani C, Inc M. Fundamental solutions of anomalous diffusion equations with the decay exponential kernel. Math Meth Appl Sci. 2019;42(11):4054–60.10.1002/mma.5634Search in Google Scholar

[6] Gao B, Langer S, Correy PM. Application of the time-dependent Green’s function and Fourier transforms to the solution of the bioheat equation. Int J Hypertherm. 1995;11:267–85.10.3109/02656739509022462Search in Google Scholar PubMed

[7] Yeung C, Atalar E. A Green’s function approach to local RF heating in interventional MRI. Med Phys. 2001;28:826–32.10.1118/1.1367860Search in Google Scholar PubMed

[8] Deng Z, Liu J. Analytical study on bioheat transfer problems with spatial or transient heating on skin surface or inside biological bodies. J Biomech Eng. 2002;124:638–49.10.1115/1.1516810Search in Google Scholar PubMed

[9] Wang H, Burgei W, Zhou H. A concise model and analysis for heat-induced withdrawal reflex caused by millimeter wave radiation. Am J Oper Res. 2020;10:31–81.10.4236/ajor.2020.102004Search in Google Scholar

[10] Durkee JW, Antich PP, Lee CE. Exact-solutions to the multiregion time-dependent bioheat equation 1. Solution development. Phys Med Biol. 1990;35:847867.10.1088/0031-9155/35/7/004Search in Google Scholar PubMed

[11] Durkee JW, Antich PP. Exact-solution to the multiregion time-dependent bioheat equation with transient heat-sources and boundary-conditions. Phys Med Biol. 1991;36:345368.Search in Google Scholar

[12] Fritz J. Partial differential equations. 4th edn. New York: Springer-Verlag; 1982.Search in Google Scholar

Received: 2020-08-11

Revised: 2020-09-16

Accepted: 2020-10-15

Published Online: 2020-12-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2020-0197

Keywords for this article

Pennes’ bioheat equation; blood flow cooling; surface cooling; analytical solution; pulsed heating

Creative Commons

BY 4.0