Vol.2, No.3, 131-138 (2010) Natural Science
http://dx.doi.org/10.4236/ns.2010.23022
Investigation of nonlinear temperature distribution in
biological tissues by using bioheat transfer equation of
Pennes’ type
Ahmed Lakhssassi1, Emmanuel Kengne1*, Hicham Semmaoui2
1LIMA – Laboratoire d’Ingénierie des Microsystèmes Avancés, Département d’informatique et d’ingénierie, Université du Québec
en Outaouais, Gatineau, Canada;
2Department of Electrical Engineering, Polytechnique, University of Montrea, Montreal, Canada
Received 14 December 2009; revised 8 January 2010; accepted 30 January 2010.
ABSTRACT
In this paper, a two level finite difference scheme
of Crank-Nicholson type is constructed and used
to numerically investigate nonlinear temperature
distribution in biological tissues described by
bioheat transfer equation of Pennes’ type. For
the equation under consideration, the thermal
conductivity is either depth-dependent or tem-
perature-dependent, while blood perfusion is
temperature-dependent. In both cases of depth-
dependent and temperature-dependent thermal
conductivity, it is shown that blood perfusion
decreases the temperature of the living tissue.
Our numerical simulations show that neither the
localization nor the magnitude of peak tempera-
ture is affected by surface temperature; however,
the width of peak temperature increases with
surface temperature.
Keywords: Bioheat Transfer Equation of Pennes
Type; Pennes’ Bioheat Transfer Equation; Modified
Crank-Nicholson Method; Two-Level Finite
Difference Scheme
1. INTRODUCTION
The evaluation of thermal conductivities in living tis-
sues is a very complex process which involves several
phenomenological mechanisms, such as heat transfer
due to perfusion of the arterial-venous blood through
the pores of the tissue (blood convection), heat con-
duction in tissues, metabolic heat generation and ex-
ternal interactions, such as electromagnetic radiation
emitted from cell phones, evaporation, metabolism, etc.
The heat transfer mechanism in biological tissues is
important for therapeutic practices, such as cancer hy-
perthermia, burn injury, brain hypothermia resuscita-
tion, disease diagnostics, cryosurgery, etc.
Many of the bioheat transfer problems have been
modeled using Pennes’ equation [1], which accounts
for the ability of tissue to remove heat by both passive
conduction (diffusion) and perfusion of tissue by blood.
Perfusion is defined as the nonvectorial volumetric
blood flow per tissue volume in a region that contains
sufficient capillaries that an average flow description
is considered reasonable. Pennes’ model was adapted
by many biologists for the analysis of various heat
transfer phenomena in a living body [2-7]. Others,
after evaluations of the Pennes model in specifical
situations, have concluded that many of the hypotheses
(foundational to the model) are not valid. Then these
latter modified and generalized the model to adequate
systems [8-11]. To analyze nonlinear temperature dis-
tribution in living tissues, it is sometime useful to
modify the one-dimensional (1D) Pennes’ bioheat trans-
fer equation by adding nonlinear terms, which gener-
ally account for temperature-dependent variability in
tissue perfusion (see for example [12]).
To treat the system of motion in living bodies, we
have written the transient 1D bioheat transfer type
model in a generalized form as follows [13],

 ,0,,
*
**
Lxtxp
QuTc
x
u
k
xt
u
cmbmb


















(1)
where
x
is the distance from the surface to the body
core (in m), t is the time (in s), kc ,,

are the density
(in gm/m³), specific heat (in cal/gm°C) and thermal
conductivity of tissue (in cal/m×s×°C), respectively and
b
c, is the specific heat of blood (in cal/gm °C), b

is
the density of blood (in gm/m³), m
Q is the metabolic
heat production per volume, m

is the nondirectional
mass flow (in gm×s/m³) associated with perfusion so
that




/
bm b
Tc Tc



is the perfusion coeffi-
cient, and

 txp , is the heat deposited per volume due
*E. Kengne dedicates this work to his son Kengneson Fred Jake Sado

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
132
to spatially distributed heating. *
s
uTT is the ele-
vated temperature (in °C) in the
x
-direction, where
T
represents the temperature distribution (in °C) in
the
x
-direction and s T represents the skin’s steady
state temperature (in °C). L is the distance (in m)
between the skin surface and the body core. In this
general form, m

is a function of temperature to in-
clude the specific case of temperature dependent per-
fusion. We assume that the thermal conductivity k
satisfies k

  , where

and

are two posi-
tive constants.
The form of the perfusion coefficient

T  , which
incorporates the heat capacity and density of blood and
tissue, depends on the characteristics of the response.
For a temperature response without an overshoot, the
perfusion is defined as a spatially averaged constant,
i.e.,

0,T

  where 0

is the steady-state per-
fusion at a constant heat flux. With an overshoot re-
sponse, we assume that the perfusion increases with
time because of vasodilation and capillary recruitment
with a higher temperature.
For Pennes’ model

,pxt is constant andm


/,
bb


where b

is the blood perfusion rate. In
this case the second term on the right-hand side of the
state Eq.1 describes the heat transport between the tissue
and the microcirculatory blood perfusion. In biological
modeling, the nondirectional mass flow m

associated
with perfusion is of form m


/
bb

 

F
T * /uand

,pxt is null, where


F
T is a function which can
be chosen as a polynomial function of the temperature.
The model used in this paper is the following
nonlinear bioheat transfer equation
,0,
2
1
*
**
LxQ
c
Q
TT
uc
x
u
k
xt
u
c
m
bb
m
s
bbb




























(2)
which is a special case of Eq.1 In this equation b

and 1

are referred to as the temperature independent
(basal) perfusion component and the temperature de-
pendent (vasodilation and angiogenesis) perfusion
component, respectively. This model is based on a
one-dimensional Pennes bioheat transfer equation and
is modified to account for temperature-dependent
variability in tissue perfusion [12]. Assuming the
metabolic heat production m
Q to be constant and
denoting ,/
*
bbm cQuu

 we obtain the following
simplified form of Eq.2.
,0,
2Lxuu
x
u
k
x
t
u














(3)
where ,/1 c



,/cc bbb





and .0/
1 c



Then the perfusion coefficient is

.uu




 Be-
side the differential equation, initial and boundary
conditions determine the temperature distribution.
Eq.3 can be written, in steady-state form, as
x



u
kx






2,uu

 which, in the case of constant
thermal conductivity k, reads
.
2
2
2
u
k
u
k
dx
ud




 (4)
The first integral of Eq.4 is given by
,
63
223
2
k
u
k
u
kdx
du






 




 (5)
where

is a constant of integration. Eq.5 is an el-
liptic ordinary differential equation. This equation can
be solved using Weierstrass’ elliptic function method
[14]. For either 23/2
  or

3 200

/2 343 ,

its particular solutions are

x
k
x
k
xu






4
cos2
4
cos23
2
2








and
 ,
8
5
,
7
2
14
8
5
,
7
2
59
2
2


























x
k
cn
x
k
cn
xu






respectively, where

 mxcn , denotes the Jacobi ellitic
cosine function with modulus m.
Mathematically, solving bioheat transfer equation
requires knowledge of both the initial and boundary
conditions applicable to the case under study. The ini-
tial conditions describe the state of the biological sys-
tem at time, 0

t while the boundary conditions give
information both on the surface 0xand at depth
.Lx

In this work, we use the following initial and
boundary conditions. The initial condition is that a
biological tissue has a uniform temperature at the
steady-state temperature of the biological tissue as in
 .;0 xuut b 

(6)
Boundary conditions include a constant temperature

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
133
on the tissue surface whereas no heat flow is assumed
at the boundary :Lx 



0,0
,
,,0 


t
x
tLu
utu surf (7)
s
urf
u being the surface temperature of the biological
tissue.
Restricting ourselves to nonnegative

,,uxt we
follow Zhao et al. [15] and construct a two level finite
difference scheme for (3) which has the same order of
accuracy. As in [15] our scheme requires only one ini-
tial condition and is also unconditionally stable and
convergent. The rest of the paper is organized as fol-
lows: in Section 2 we construct a finite difference
scheme for the modified Pennes Eq.3 and prove its
solvability. Numerical experiments are reported in
Section 3, while a brief summary of the results is done
in Section 4.
2. NUMERICAL SCHEME
The main difficulties in numerical solving the bioheat
transfer Eq.3 are the nonlinearity due to the perfusion
term and the different material properties of the tissue.
Using a modified Crank-Nicholson method (see the
description below), we are able to integrate the heat
equation efficiently. Within this approach, the approxi-
mate elevated temperature n
u at time n t is con-
structed by a combination of previous elevated tempera-
ture 1n u at time 1.
n
t
For the numerical computation of the nonlinear heat
transfer Eq. 3 with initial and boundary conditions (6)
and (7) we use h to represent the space mesh and

to
represent the time step so that
x
will be incremented by
h and
t
by .

We choose an integer
N
and h
such that .NhL We construct a finite difference
scheme for (3) by using the second-order central differ-
ence scheme in space and a scheme of Crank-Nicholson
[16] type in time. In what follows, n
j
u denotes the ap-
proximate elevated temperature at time n
t at depth :
j
x

,.
nn
jj
uuxt If
x

denotes the central difference
operator, the discretization scheme is then

1
11 1
22 2
nn
jj nnnnn
jjj jj
uu uuuuu




 

  


1
1/2 1/2
,,,
nn
xjnxjxjnxj
kxt ukxt u 





1
00
0,
, 0.
nn
nNN
jN
uu
uu h




The difference scheme then takes the form


, ,
,
22
22
1
22
2
22
1
2
1
100
1
1
2
1
1
1
2
1
1
1
1
1
1
2
1
1
2
1
1
2
1
1
1
11
2
1











































n
N
n
N
n
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
n
j
uuuu
u
h
k
u
h
kk
u
uu
h
k
u
h
k
h
kk
uu
h
k






(8)
where


.,  ntjhxkk n
j
n
j The truncation error of
this difference schemes in of order


22
 hO for each
interior grid point


1,, ntxn
j and .0 Nj

 In
matrix form, the difference scheme (8) reads
..., 2, 1,n ,uu1 n
R
n
Lnn QQ (9)
where the vector


T
n
N
nnnnuuuu ,...,,, 210
u represents the
numerical solutions at the time level )),,(( n
j
n
j
ntxut 
n
L
Q (left-hand side matrix depending on n) and n
R
Q
(right-hand side matrix depending on n) are both
square tridiagonal matrices of dimension N+1 of the
form





















, and other for ,0
,1 if ,)2/(
,
2
22
1
,1 if ,)2/(
},1,1{ if ,1
2
2
1
1
1
1
2
qp
pqkh
h
kk
u
pqkh
Nqp
Q
n
p
n
p
n
p
n
p
n
p
Ln

































. and other for ,0
,1 if ,)2/(
, if ,
2
22
1
,1 if ,)2/(
,1 if ,1
,1 if ,0
12
2
11
1
1
1
1
1
2
qp
pqkh
pq
h
kk
u
pqkh
Nqp
qp
Q
n
p
n
p
n
p
n
p
n
p
Rn




In other words,
0
22
1
2
0
01
1
11
2

















n
j
n
j
Luk
h
Qn


A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
134

,
10
22
0
1
22
1

















n
j
n
j
n
jk
hh
kk


0
22
1
2
0
00
1
1
1
1
2



















n
j
n
j
Ruk
h
Qn


.
10
22
0
1
1
22
11
1

















n
j
n
j
n
jk
hh
kk


Theorem (on the solvability of system (9)). System
(9) is unconditionally solvable for each time step .n
Proof. It is sufficient to show that matrix n
L
Q is
inversible. From the expression for n
L
Q we have
,0
22
1
;01Q ;01
1
1
1
,1
1
!11
L
1
2
111
n












n
p
N
pqq
pq
L
pp
L
N
q
qN
L
NN
N
q
q
L
uQQ
QQQ
nn
nn n
L

So that matrix n
L
Q is diagonally dominant. By Ger-
shgorin’s theorem [17] we conclude that the matrix n
L
Q
is invertible. This proves the theorem.
3. NUMERICAL EXPERIMENTS
For the numerical experiments we use mL 01208.0
and 4200000/1

[18]. Because our main aim in this
work is the impact of the nonlinear term on the tem-
perature distribution, we will work with different choices
of nonlinear coefficient .

As initial elevated tempera-
ture

,xuB we use one of the functions


 1
0
uxuB
 
u
xx/exp and  

 ,/exp1/
0
uB xxuxu where
0
u and u x are two constants. The thermal conductiv-
ity of tissue k is taken to be either a function of ,
x
namely,
 9.01.11210000 2 xxxk or a linear func-
tion of the elevated temperature, namely, ,
10ukkk


with 0.4574k0 and 0.001403k1 (see [19]). To
analyze the effects of the blood perfusion, we used
the following three couples



, of parameters



 ,01.0,009.0 ,010.0005,0.0 and  .02.0,005.0
3.1. Effects of Blood Perfusion, Thermal
Conductivity, and Initial Elevated
Temperature
To analyze the effects of the blood perfusion, we used
the following three couples



, of parameters



 ,01.0,009.0 ,010.0005,0.0 and  .02.0,005.0 The ca-
lculation results of the influences of the blood perfu-
sion, the thermal conductivity, and the initial elevated
temperature on the temperature distribution are shown
in Figures l-4. Figures 1 and 3 and Figures 2 and 4
are obtained with the initial conditions


01
B
ux u


1exp /
u
x
x



, respectively, with1/200,
u
x
01 6,u

and 02 u

12.928. For all these figures, we used


0.
surf B
uu Plot (a) gives the elevated temperature
profile along the
x
direction at time 5.2,ts while
plot (b) shows the elevated temperature profile along the
t
direction at depth .01148.0 mx

Plots (c) and (d)
are obtained with

 
, 0.005,0.02   and show the
elevated temperature along the
x
-direction for different
times
t
and along the
t
-direction for different depths
x
, respectively. Plots (e) and (f) show the elevated tem-
perature for temperature-dependent (plot (1)) and
depth-dependent (plot (2)) thermal conductivity .k
Here, we used the parameter
 
, 0.005,0.02  
Figures 1 and 2 give the elevated temperature when
the thermal conductivity of tissue k is a function of dis-
tance x. In Figure 1(a) and (b), as well as in Figure 2(a)
and (b), the results show that the higher the blood perfu-
sion is, the lower is the elevated temperature. In other
words, to decrease the temperature in the living tissue, it
is sufficient to increase the blood perfusion.
As in the case of
x
-dependent thermal conductivity,
plots (a) and (b) of Figures 3 and 4 show that, in the
case of temperature-dependent thermal conductivity, the
higher the blood perfusion is, the lower is the elevated
temperature. We may then conclude from Figures 1-4
that the effect of the blood perfusion is to decrease the
temperature in the living tissue.
Figures 1(c) and (d) and Figure 3(c) and (d) show
that the peak temperature is a decreasing function of
both depth and time. When the initial temperature of the
living tissue is a decreasing function of depth x (Figures
1 and 3), a phenomenon of heating (in the sense that the
temperature is almost above the initial temperature) is
observed both in depth and time, and the peak tempera-
ture at any depth is above the initial temperature (see
plots (c) and (d) of Figures 1 and 3). To the contrary, a

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
135
Figure 1. Effect of blood perfusion on the temperature in the case of depth-dependent ther-
mal conductivity of tissue ,)]200exp(1[6)()0,( Cxxuxu B with decreasing initial
temperature. Plots of the first row show the temperature elevation for different blood perfu-
sion while plots of the second give the temperature elevation for the same blood perfusion at
different time t (c) and at different depth
x
(d).

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
136
Figure 2. Effect of blood perfusion on the temperature in the case of depth-dependent thermal
conductivity of tissue with increasing initial temperature


928.12)0,( 

x
u
x
uB
.)]200exp(1[ 1Cx  Plots of the first row show the temperature elevation for different blood
perfusion while plots of the second and third row give the temperature elevation for the same
blood perfusion at different time
t
((c) and (e)) and at different depth
x
((d) and (f)). Plots (e) and (f)
of the third row indicate the elevated temperature obtained when using temperature-dependent and
depth-dependent thermal conductivity, respectively.
Figure 3. Effect of blood perfusion on the temperature in the case of temperature dependent
thermal conductivity of tissue with decreasing initial temperature
  xuxu B )0,(
.
)
]200exp
(
1[928.12
C
x


Plots of the first row show the temperature elevation for different
blood perfusion while plots of the second and third row give the temperature elevation for the
same blood perfusion at different time
t
((c) and (e)) and at different depth
x
((d) and (f)). Plots
(e) and (f) of the third row indicate the elevated temperature obtained when using tempera-
ture-dependent and depth-dependent thermal conductivity, respectively.

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
137
Figure 4. Effect of blood perfusion on the temperature in the case of temperature-dependent ther-
mal conductivity of tissue with increasing initial temperature .)]200exp(1/[928.12)( CxxuB


Plots of the first row show the temperature elevation for different blood perfusion while plots of
the second row give the temperature elevation for the same blood perfusion at different time
t
(c)
and at different depth
x
(d).
Figure 5. Effect of surface temperature on the temperature distribution, location, magnitude,
and width of peak temperature.
cooling phenomenon (in the sense that the temperature is
almost below the initial temperature) is observed when
the initial temperature increases as a function of depth
x
(Figures 2 and 4).
Plots (e) and (f) of Figure 2 and Figure 3 show that
the elevated temperature, numerically obtained in the
case of temperature-dependent thermal conductivity (pl-
ots (1)) is larger than the one obtained in the case of

A. Lakhssassi et al. / Natural Science 2 (2010) 131-138
138
depth-dependent thermal conductivity (plots (2)).
3.2. Effect of Surface Heating
The effects of surface heating (when the surface of the
living tissue is maintained at different constant tempera-
ture) on the temperature distribution and on the location
of the peak temperature have been studies and observed
when the initial temperature is either a decreasing func-
tion of depth
x
(Figure 5(a) and (b)) or an increasing
function of depth
x
(Figure 5(c) and (d)). The plots
of this figure are obtained with 1/ 4200000

 and
(, )(0.005,0.02).

 Figure 5(a) gives the depth at
which the temperature peak occurs at a given time while
Figure 5(b) gives the time at which the temperature peak
at a given depth occurs. Figure 5(a) shows that the tem-
perature peak does not depend on the surface tempera-
ture. These two plots also show that the magnitude of
peak temperature does not depend on the surface tem-
perature: the magnitude of peak temperature is the same
for all three surface temperatures. As we can see from
Figures 5(c) and (d), the depth (Figure 5(c) and the
time (Figure 5(d) from which the variation of the tem-
perature suddenly changes do not depend on the surface
temperature. Many other information about the tempera-
ture distribution in the living tissue can be obtained from
Figure 5. For example, it is seen from plots (a) and (b)
that the width of the peak temperature increases with the
surface temperature: the higher the surface temperature is,
the higher is the width of the peak temperature.
4. CONCLUSIONS
In this work, we present a numerical investigation of a
1D bioheat transfer equation with either depth-dependent
or temperature-dependent thermal conductivity and with
temperature-dependent blood perfusion. An implicit un-
conditional numerical scheme of the Crank-Nicholson
type is constructed and used to solve the nonlinear bio-
heat transfer equation with given initial and boundary
conditions. We found that blood perfusion decreases the
temperature in living tissue. It is also shown that the lo-
calization and the magnitude of peak temperature do not
depend on the surface temperature, while its width in-
creases with surface temperature. The computation pre-
sented in this paper can be used to predict the tempera-
ture distribution in living tissue during many bioheat
transfer processes. The method used in this paper has
potential to provide the temperature distribution in tissue
in the absent of heat flux.
REFERENCES
[1] Pennes, H.H. (1948) (1998) Analysis of tissue and arterial
blood temperatures in the resting human forearm. Journal
of Applied Physics, 1, 93-122; 85, 5-34.
[2] Meyer, C., Philip, P. and Troltzsch, F. (2004) Optimal
control of semilinear PDE with nonlocal radiation
interface conditions. IMA Preprint Series, 2002.
[3] Yamamoto, M. and Zou, J. (2001) Simultaneous rec-
onstruction of the initial temperature and heat radiative
coefficient. Investment Problems, 17, 1181-1202.
[4] Ganzler, T., Volkwein, S. and Weiser, M. (2006) SQP
methods for parameter identification problem arising in
hyperthermia. Optimization Methods and Software, 21,
869-887.
[5] Arkin, H., et al. (1986) Thermal pulse decay method for
simultaneous measurement of local thermal conductivity
and blood perfusion: A theoritical snalysis. Journal of
Biomechanical Engineering, 108, 208-214.
[6] Deuflhard, P. and Seebass, M. (1998) Adaptive multilevel
FEM as decisive tools in clinical cancer therapy hyper-
thermia. Konrad-Zuse-Zentrum für Information stechn,
Berlin Takustr, Berlin, 7.
[7] Hill, J. and Pincombe, A. (1992) Some similarity temp-
erature profiles for the microwave heating of a half-space.
Journal of the Australian Mathematical Society, Series B,
33, 290-320.
[8] Chato, J.C. (1980) Heat transfer to blood vessels. Journal
of Biomechanical Engineering, 102, 110-118.
[9] Weinbaum, S. and Jiji, L.M. (1985) A two simplified
bioheat equation for the effect of blood flow on average
tissue temperature. Journal of Biomechanical Engineer-
ing, 107, 131-139.
[10] Chen, M.M. and Holmes, K.R. (1980) Microvascular
contributions in tissue heat transfer. Annals of the New
York Academy of Sciences, 335, 137-150.
[11] Charney, C.K. (1992) Mathematical models of bioheat
transfer. Advanced Heat Transfer, 22, 19-155.
[12] Davies, C.R., Saidel, G.M. and Harasaki, H. (1997)
Sensitivity analysis of 1-D heat transfer in tissue with
temperature-dependent perfusion. Journal of Biomechan-
ical Engineering, 119, 77.
[13] Lang, J., Erdmann, B. and Seebass, M. (1999) Impact of
nonlinear heat transfer on temperature control in regional
hyperthermia. IEEE Transactions on Biomedical Engin-
eering, 46, 1129-1138.
[14] Weierstrass, K. (1915) Mathematische Werke V, New
York, Johnson, 4-16; Whittaker, E.T. and Watson, G.N.
(1927) A Course of Modern Analysis. Cambridge
University Press, Cambridge, 454.
[15] Zhao, J.J., Zhang, J., Kang, N. and Yang, F. (2005) A two
level finite difference scheme for one dimensional
Pennes’ bioheat equation. Applied Mathematics and
Computation, 171, 320-331.
[16] Press, W.H., Flannery, B.P., Teukolsky, S.A. and
Vallerling, W.T. (1992) Numerical Recipes in Fortran.
Cambridge University Press, Cambridge.
[17] Heath, M.T. (2002) Scientific Computing, an Introductory
Survey, second edtion, McGraw-Hill, New York.
[18] Liu, J., Chen, X. and Xu, L.X. (1999) New thermal wave
aspects on burn evaluation of skin subjected to ins-
tantaneous heating. IEEE Transactions on Biomedical
Engineering, 46, 420-428.
[19] Tunç, M., Çamdali, Ü., Parmaksizoglu, C. and Cikrik çi,
S. (2006) The bioheat transfer equation and its applica-
tions in hyperthermia treatments. Engineering Computa-
tions, 23, 451-463.