Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

52, 4, pp. 947-957, Warsaw 2014

MODELLING OF TISSUE THERMAL INJURY FORMATION PROCESS
WITH APPLICATION OF DIRECT SENSITIVITY METHOD

Marek Jasiński

Silesian University of Technology, Institute of Computational Mechanics and Engineering, Gliwice, Poland

e-mail: marek.jasinski@polsl.pl

In the paper, numerical analysis of thermal processes proceeding in a biological tissue is
presented. The tissue is subjected to the external heat impulse and a 2D problem is taken
into account. In order to determine the influence of variations of thermophysical parameters
of the tissue on the value of tissue injury integral and the area of the lesion, adirect approach
of sensitivity analysis is applied. The process of thermal injury formation is also analyzed.
At the stage of numerical simulation, the boundary element method is used. In the final
part of the paper, an example of numerical simulation is shown.

Keywords: bioheat transfer, Arrhenius scheme, boundary element method, sensitivity
analysis

1. Introduction

The phenomenon of biological tissue damage due to external thermal factors is one of the quite
frequently encountered interaction between a living organism and its surroundings. It may be
either an accidental case, as in thermal burns, or a fully-controlled case, as in the treatment
using thermal effects. The impact of this type of damage on the living organism is dependent on
its extent anddepth. It should be pointed out that an elevated temperature in a biological tissue
doesnotalways lead to irreversible changes.Whenthe influenceof temperature ismoderate (37◦-
45◦C), the main response of the tissue is dilatation of blood vessels, which is not accompanied
by the tissue injury. Of course, an important factor is the exposure time of the thermal pulse,
because the interaction with a suitably long (approximately 6 hours) temperature of 42◦Cmay
lead to epidermis necrosis (Henriques, 1947; Torvi andDale, 1994). Therefore, one can conclude
that if the tissue damage exceeds some threshold value, the tissue injury becomes irreversible,
otherwise the tissuehas theability to returnto its original, healthy state. It shouldbeemphasized
that the change in the area of the tissue thermal injury (its expansion orwithdrawal)may occur
even some time after the cessation of external heat impulse. Thus, determination of the time
after which the wound reaches its final size may be associated with certain difficulties.
In the numerical analysis of the tissue thermal damage process one of the most commonly

used approaches is the so-called Arrhenius injury integral, which allows one to determine the
degree of tissue damage, and furthermore, to define some of the thermophysical parameters
of tissue as dependent of the tissue injury. It is worth to mention the work by Abraham and
Sparrow (2007) inwhich the perfusion coefficient is expressedbyapolynomial function reflecting
the initial increase in blood flow resulting in vasodilatation as well as blood flow decrease as
the vasculature begins to shut down, whereas Glenn et al. (1996) used exponential functions to
model changes in scattering coefficient during laser-tissue interaction.
The classic definition of the tissue injury integral implies that even at a small local increase

in temperature damage of the tissue is irreversible. In the current paper, the algorithmproposed
by Jasiński (2013) is used as amodel of tissue damagewithdrawal at those points of the domain
analyzed, in which the value of the injury integral does not exceed a certain damage threshold.



948 M. Jasiński

The injury formation process analysis (IFPA algorithm, Jasiński, 2012) is applied in order to
determine the time after which the lesion is fully formed and the area of different stages of tissue
damage is based on the Arrhenius integral as well.

One of the problems connected with the application of the mathematical model is the sen-
sitivity of the solution with respect to the parameters appearing in the governing equations.
The sensitivity information may be used, among others, to analyze the influence of the change
of parameters on the final solution of the problem being considered. Such a kind of problems
was byDavies et al. (1997) and Jasiński (2009), whileMajchrzak and Jasiński (2004) presented
sensitivity analysis of the skin burns prediction model proposed by Henriques (1947).

2. Governing equations

The transient heat transfer in the 2D domain of a homogeneous biological tissue of rectangular
shape (Fig. 1.) is described by the Pennes equation in form (Majchrzak et al., 2011; Mochnacki
and Piasecka-Belkhayat, 2013)

x∈Ω : cṪ =λT,ii+QV (2.1)

where λ [Wm−1K−1] is the thermal conductivity, c [Jm−3K−1] is the volumetric specific heat,
QV [Wm

−3] is the internal heat source, while T = T(x, t) and Ṫ denotes temperature and its
time derivative.

The component QV comprises the information of the internal heat sources and is described
as

QV =Qperf +Qmet = cBGB(TB −T)+Qmet (2.2)

where GB [(m
3
blood/s)/(m

3
tissue)], cB [Jm

−3K−1] and TB correspond to the perfusion coeffi-
cient, the volumetric specific heat of blood and the artery temperature, respectively, while
Qmet [Wm

−3] is the internal metabolic heat source (Majchrzak et al.,2003; Torvi and Dale,
1994).

Fig. 1. The considered domain

Equation (2.1) is supplemented by appropriate boundary conditions:
— on the external surface of tissue (Γ0)

x∈Γ0 : q(x, t)=

{

q0 t¬ texp

α(T −Tamb) t> texp
(2.3)

— on the remaining parts of the boundary (Γc)

x∈Γc : q(x, t)= 0 (2.4)



Modelling of tissue thermal injury formation process... 949

where q0 [Wm
−2] is the knownboundaryheatflux, α [Wm−2K−1] is the convective heat transfer

coefficient and Tamb is the temperature of surroundings, while texp is the exposure time. The
initial distribution of temperature is also known.

It should be pointed out that the heat flux along the boundary Γ0 is assumed to be irregular
as in the majority of non-controlled cases of high temperature-biological tissue interactions. In
this paper, it is assumed as a polynomial function of the 7th degree with the distribution of the
heat flux visible in Fig. 1.

The damage of a biological tissue resulting from temperature elevation is usually modeled
by the so-called Arrhenius injury integral defined as (Abraham and Sparrow, 2007; Oden et al.,
2007)

θ(x)=

tF
∫

0

Aexp
(

−
∆E

RT

)

dt (2.5)

where A is thepre-exponential factor [s−1], ∆E is the activation energy [Jmole−1], R is the uni-
versal gas constant [Jmole−1K−1], and [0, tF ] is the considered time interval, while the criterion
for tissue necrosis is

θ(x)­ 1 (2.6)

According to thenecrotic changes in tissue, thebloodperfusioncoefficient isdefinedas (Abraham
and Sparrow, 2007)

GB =GB(θ)=























GB0 θ=0

(1+25θ−260θ2)GB0 0<θ¬ 0.1

(1−θ)GB0 0.1<θ¬ 1

0 θ> 1

(2.7)

where GB0 is the initial perfusion coefficient. The values of coefficients for the interval from
0 to 0.1 correspond to the increase in the perfusion coefficient caused by vasodilatation up to
the value θ=0.05 (maximum of the function) and the beginning of narrowing of blood vessels
(between 0.05 and 0.1). The interval 0.1 to 1 reflects the blood flowdecrease as the vasculature
going to shut down.

On the basic on the injury integral in form (2.5), the damage fraction FD is calculated as
(Glenn et al., 1996; Oden et al., 2007)

FD(x)= 1−exp(−θ) (2.8)

3. Modeling of the thermal damage of tissue

The assumption of theArrhenius scheme is that the damage of the tissue is irreversible. In order
to consider that the tissue could get back to its native state after the thermal impulse is ceased,
the following algorithm is proposed (Jasiński, 2013), see Fig. 2.

Let us assume that for the time interval [0, tF ] being under consideration and divided in-
to F subintervals [tf−1, tf] (where f = 1,2, . . . ,F), the values T0(x), . . . ,Tf(x) as well as
θ0(x), . . . ,θf−1(x) at the point x∈Ω are known.At the same time, the recovery threshold θrec
is accepted.

If the injury integral at the point x for time tf achieves the value equal or greater than θrec
then the injury of the tissue becomes irreversible. Otherwise, the function denoted as θapp(x,T)



950 M. Jasiński

Fig. 2. The proposed algorithm of tissue injury calculation

is introduced in order tomodel thewithdrawal of the tissue injury. In currentpaper, it is assumed
as a linear one between (T0,θ0) and (Tf−1,θf−1)

θapp(x,T)=m1+m2T (3.1)

For the time tf+1, if Tf+1 <Tf and θf <θrec, again the function θapp defined for the time t
f

is used.

The damage fraction FD (Eq. (2.8)) is used in the thermal injury formation process analysis
(IFPA) algorithm. In this algorithm, five intervals of values for FD (denoted as tha) have been
distinguished:

• tha 1: [0, 0.01),

• tha 2: [0.01, 0.05),

• tha 3: [0.05, 0.63),

• tha 4: [0.63, 0.99),

• tha 5:­ 0.99.

The values in the intervals are interpreted as:

• 0.01: up to this value, the tissue is in its normal state, so the value could be named as the
border of thermally untouched tissue,

• 0.05: the border of vasodilatation – arises from the polynomial function for GB (c.f. Eq.
(2.7)); at this value of FD the perfusion coefficient has its maximum value,

• 0.63: corresponds to the criterion of tissue necrosis (c.f. Eq. (2.6)),

• 0.99: could be treated as the criterion of complete tissue destruction.

InFig. 3, the concept of IFPAalgorithm is presented.At first, the comparison of tha intervals
achieved on an element of the domain considered for two successive time steps is made. Next,
only the elements which have changed the intervals are selected. Finally, the bar chart with the
number of elements achieving the individual interval is drawn.



Modelling of tissue thermal injury formation process... 951

Fig. 3. The concept of the IFPA algorithm

It should be pointed out that if we add up the number of elements which have swapped from
one interval to another in successive time intervals, we can obtain the knowledge of dynamics
of the thermal injury formation process. Additionally, changes of the sum of elements with tha
interval greater than 1 give information about the proliferation of the lesion.

4. Sensitivity analysis

To determine the influence of thermophysical parameters on the value of the injury integral, a
direct approach of the sensitivity analysis has been applied (Dems, 1986; Kleiber, 1997).

According to the rules of the direct method, Eq. (2.1) is differentiated with respect to the
parameter ps, where s = λ,c,GB0 or Qmet. So, the number of additional sensitivity tasks
corresponds to the number of parameters to which the sensitivity analysis is done (Majchrzak
and Jasiński, 2002; Majchrzak and Kałuża, 2006). If we denote

Us =
∂T

∂ps
(4.1)

as a sensitivity function, then one can write, c.f. Eq. (2.1), (Mochnacki andMajchrzak, 2003)

x∈Ω : cU̇s =λUs,ii+Q
s
V (4.2)

where

U̇
s =
∂Us

∂t
U
s
,ii =
∂T,ii

∂ps
(4.3)

while

QsV =
[cBGB0f(θ)

λ

∂λ

∂ps
− cBf(θ)

∂GB0

∂ps
− cBGB0

∂f(θ)

∂ps

]

(T −TB)

− cBGB0f(θ)U
s+
(c

λ

∂λ

∂ps
−
∂c

∂ps

)

Ṫ −
Qmet

λ

∂λ

∂ps
+
∂Qmet

∂ps

(4.4)

The variation of θ is calculated as (c.f. Eq. (2.5))

∂θ

∂ps
=

tF
∫

0

A
∆EUs

RT2
exp
(

−
∆E

RT

)

dt (4.5)



952 M. Jasiński

The change of the injury integral due to changes of the parameters ps can be estimated using
the following formula

∆θ(x)=

√

√

√

√

n
∑

s=1

(∂θ(x)

∂ps
∆ps

)2

(4.6)

Finally, using Eq. (4.6), the value of Arrhenius integral is calculated as

θ(x,ps±∆ps)= θ(x,ps)±∆θ(x) (4.7)

5. Boundary element method

The primary and also the additional problems resulting from the sensitivity analysis have been
solved using the 1st scheme of the BEM for 2D transient heat diffusion problems (Brebia and
Dominquez, 1992). So, the following equation will be considered

x∈Ω : cḞ =λF,ii+S (5.1)

where F =F(x, t) denotes the temperature or functions resulting from the sensitivity analysis
(c.f. equations (2.1) and (4.2)), while S = S(x, t) is the source function (c.f. Eqs. (2.2) and
(4.4)).
For the timegridwitha constant time step ∆t, theboundary integral equation corresponding

to the transition tf−1 → tf is of the form

B(ξ)F(x, tf)+
1

c

tf
∫

tf−1

∫

Γ

F
∗(ξ,x, tf, t)J(x, t) dΓ dt=

1

c

tf
∫

tf−1

∫

Γ

J
∗(ξ,x, tf, t)F(x, t) dΓ dt

+
1

c

tf
∫

tf−1

∫∫

Ω

F∗(ξ,x, tf, tf−1)F(x, tf−1) dΩdt+
1

c

tf
∫

tf−1

∫∫

Ω

S(x, t)F∗(ξ,x, tf, t) dΩdt

(5.2)

In equation (5.2), F∗ is the fundamental solution

F∗(ξ,x, tf, t)=
1

4πa(tf − t)
exp
[

−
r2

4a(tf − t)

]

(5.3)

where r is the distance from the point under consideration x to the observation point ξ, while

J
∗(ξ,x, tf, t)=−λF∗(ξ,x, tf, t),ini (5.4)

and B(ξ) is the coefficient from the interval (0,1).
In this paper, constant boundary elements have been used. Details concerning numerical

realization of theBEMcanbe found, amongothers, inMajchrzak (2013) andPiasecka-Belkhayat
(2011).

6. Results

A domain of rectangular shape (c.f. Fig. 1) of dimensions 0.05m×0.015m is considered. The
interior of the domain has been divided into 6000 internal constant cells, while the external
boundary into 320 constant elements.



Modelling of tissue thermal injury formation process... 953

In computations, the following values of thermophysical tissue parameters have been as-
sumed: λ = 0.3Wm−1K−1, c = 3.647MJm−3K−1, GB0 = 0.00125(m

3
blood/s)/(m

3
tissue),

Qmet = 245Wm
−3, while for the blood cB = 3.9962MJm

−3K−1 and TB = 37
◦C. The pa-

rameters of the Arrhenius injury integral are: A = 3.1 · 1098 s−1, ∆E = 6.27 · 105Jmole−1,
R=8.314Jmole−1K−1 and θrec =0.05.

In the boundary condition (c.f. Eq. (2.3)), themaximumvalue of the heat flux q0 is assumed
as 20kWm−2, the exposure time is 30s, while α=10Wm−2K−1 and Tamb = 20

◦C (Jasiński,
2012).

The coordinates of the points are (c.f. Fig. 1): C1(0.01575,0.001875), D1(0.03475,0.000625),
C2(0.01575,0.003125) and D2(0.03475,0.000875), the time step equals to ∆t=1s.

As has been previouslymentioned, the sensitivity analysis is performedwith regard to ther-
mal conductivity, volumetric specificheat, initial perfusion coefficient andmetabolic heat source.
It is assumed that for all parameters ∆ps =0.1ps.

Figure 4 shows courses of the temperature as well as courses of the injury integral θ. At two
of these points, C1 and D1, the value of the injury integral is above the recovery threshold θrec.
At the point C1, the value of the injury integral is much greater than 1, so the tissue is fully
damaged, while at the point D1 the value of injury integral is 0.168 which means partly da-
maged tissue. TheArrhenius integral value at points C2 and D2 have not reached the recovery
threshold, so the functions θapp (see the algorithm in Section 3) are defined for the stage of
temperature lowering.

Fig. 4. Courses of temperature and the injury integral θ

In Fig. 5, courses of the Arrhenius integral found on the basis of sensitivity analysis are
shown (c.f. Eq. (4.7)). It is clearly visible that at every point the courses for θ−∆θ are very
close to zero, so one can say that the tissue is not destructed in this case.

One can note that at the point D2 for the case θ+∆θ, some values calculated on the basis
of formula (4.7) are above the recovery threshold and then decrease. Although these results are
correct from the sensitivity analysis point of view, the “real” course of the injury integral should
correspond to the dashed line inthe right-hand side figure. This is due to sensitivity calculation
without taking into account that according to the algorithm of Arrhenius injury calculation
presented in Section 3, the value of integral cannot be reduced when it is above the recovery
threshold.

Figure 6 presents courses of the perfusion coefficient GB (c.f. Eq. (2.7)) obtained on the
basis of sensitivity analysis of the Arrhenius integral. It is clearly visible that according to the
necrotic changes in the tissue domain, the perfusion coefficient is changing, decreasing to zero
for the point C1 in which the injury integral is greater than 1. At the point D2 marked by the
dashed line, it corresponds, as previously, to the non-reduced value of the injury integral. The
effect of vasodilatation is visible for all four points.



954 M. Jasiński

Fig. 5. Courses of the injury integral

Fig. 6. Courses of the perfusion coefficient GB

The differences in values of the injury integral (and in consequence, the damage fraction
– c.f. Eq. (2.8)) visible in Fig. 5 are substantial from the point of view of the thermal injury
formation process analysis (c.f. Section 3). The results obtained by this approach are shown in
Fig. 7. The results are presented for two cases: the Arrhenius integral calculated for the basic
thermophysical parameters (LHS) and for sensitivity analysis case θ+∆θ (RHS). One could
see that for t > texp some elements are classified into tha 1 (so, they go back into the healthy
state) only in the first one.

Fig. 7. The number of elements achieving individual intervals (LHS: basic values of thermophysical
parameters, RHS: sensitivity analysis, θ+∆θ)



Modelling of tissue thermal injury formation process... 955

In Fig. 8, courses of the thermal injury proliferation process are shown. The curves are
obtained by adding the number of elements which are in tha interval greater than 1. The
recovery stage, thatmeans the process of reducing the lesion area is the best visible for the case
of the injury integral calculated directly. The figure is scaled in the number of elements, but it
should be pointed out that these data could be very easily recalculated in the cross-section area
of the lesion using the field of single internal element (for the geometrical grid assumed in the
paper: 1.25 ·10−7m2).

Fig. 8. Proliferation of the thermal injury

7. Conclusion

Theproposedmethodof tissue injury calculation is closer to the real conditions of the interaction
between the tissue and a high-temperature impulse. Its main advantage is the possibility of
estimation of the tissue recovery area, which plays an important role in the case of modelling of
the fully controlled case of the interaction.

In the sensitivity analysis, the direct approach has been used. The impact of a 10% change
in the values of tissue thermophisical parameters to the Arrhenius integral has been examined.

It has turned out that in some points of the domain considered, the sensitivity parameter
∆θ may vary by much greater than 1, what e.g. at the point C1 determines the difference
between the irreversibly damaged tissue and the totally destructed tissue (from the point of
view of the IFPA algorithm, the difference between tha 4 and tha 5). The value of ∆θ obtained
at the point D2 (approximately 0.05) causes that the point could be classified into different
tha interval depending on the value calculated on the base of Eq. (4.7). It is clearly visible
in Fig. 5. Furthermore, the whole algorithm of tissue injury calculation should be taken into
account, otherwise the results of sensitivity analysis could give misleading information about
tissue recovery at the point, in which the injury integral value is above the recovery threshold.

It is also important due to the course of the perfusion coefficient, expressed in the present
paper as a polynomial function of the Arrhenius integral (c.f. Fig. 6), its value increases or
decreases in accordance with the increasing or decreasing injury integral.

The tissue returns to the native state at two points considered (calculations for the basic
values of thermophysical parameters, c.f. Fig. 5): at the point C1 after 383 seconds, while at
point D2 after 197 seconds.

The process of proliferation of the thermal injury shows that the most distinct recovery
stage occurred for the basic value of the injury integral (c.f. Fig. 8). The maximal area of the
lesion is 489 elements and decreases down to 437 elements after 157 seconds. The results of the
sensitivity analysis show that for the case θ−∆θ the lesion area is very small (37 elements) and



956 M. Jasiński

achieved after 13 seconds, while for the θ+∆θ, the maximal area is 605 elements between 76
and 81 seconds, then reduces to its final area equal to 591 elements.
The sensitivity analysis in combination with tha intervals and the new algorithm of tissue

injury calculation seems to be a quite convenient means of analysis of the thermal injury forma-
tion process and could give more precise data about the depth and cross-section area of injury.
It could be very important, especially in cases of a controlled coagulation process like, e.g., in
some thermotherapies.

References

1. Abraham J.P., Sparrow E.M., 2007, A thermal-ablation bioheat model including liquid-to va-
por phase change, pressure- and necrosis-dependent perfusion, andmoisture-dependent properties,
International Journal of Heat and Mass Transfer, 50, 2537-2544

2. Brebia C.A., Dominquez J., 1992, Boundary elements, an introductory course,Computational
Mechanics Publications, McGraw-Hill Book Company, London

3. Davies C.R., Saidel G.M., Harasaki H., 1997, Sensitivity analysis of one-dimensional heat
transfer in tissue with temperature-dependent perfusion, Journal of Biomechanical Engineering,
Transactions of The ASME, 119, 77-80

4. Dems K., 1986, Sensitivity analysis in thermal problems. Part I: Variation of material parameters
within fixed domain, Journal of Thermal Stresses, 9, 303-324

5. Glenn T.N., Rastegar S., Jacques S.L., 1996, Finite element analysis of temperature con-
trolled coagulation in laser irradiated tissue, IEEE Transactions on Biomedical Engineering, 43,
79-87

6. Henriques F.C., 1947, Studies of thermal injuries, V. The predictability and the significance of
thermally induced rate process leading to irreversible epidermal injury,Archives of Pathology, 43,
489-502

7. JasińskiM., 2009, Sensitivity analysis of transient bioheat transferwith perfusion rate dependent
on tissue injury,Computer Assisted Mechanics and Engineering Science, 16, 267-277

8. Jasiński M., 2012, Numerical modeling of burn wound formation process, 10th World Congress
on Computational Mechanics (WCCM 2012), Sao Paulo, Brasil, 1-10

9. Jasiński M., 2013, Investigation of tissue thermal damage process with application of direct sen-
sitivity method,MCB: Molecular and Cellular Biomechanics, 10, 3, 183-199

10. KleiberM., 1997,Parameter Sensitivity in NonlinearMechanics, J.Willey& Sons Ltd, Chicester

11. Majchrzak E., 2013, Application of different variants of the BEM in numerical modeling of
bioheat transfer processes,MCB: Molecular and Cellular Biomechanics, 10, 3, 201-232

12. Majchrzak E., Jasiński M., 2002, Numerical analysis of bioheat transfer processes in tissue
domain subjected to a strong external heat source, [In:] Boundary Elements Techniques, Z. Yao,
M.H. Aliabadi (Edit.), Tsinghua University Press, Springer

13. Majchrzak E., Jasiński M., 2004, Sensitivity analysis of burn integrals, Computer Assisted
Mechanics and Engineering Science, 11, 125-136

14. Majchrzak E., Kałuża G., 2006, Sensitivity analysis of biological tissue freezing process with
respect to the radius of spherical cryoprobe, Journal of Theoretical and Applied Mechanics, 44, 2,
381-392

15. Majchrzak E., Mochnacki B., Jasiński M., 2003, Numerical modelling of bioheat transfer in
multi-layer skin tissue domain subjected to a flash fire,Computational Fluid and Solid Mechanics,
1/2, 1766-1770

16. Majchrzak E., Mochnacki B., Dziewoński M., Jasiński M., 2011, Numerical modelling of
hyperthermia and hypothermia processes,Advanced Materials Research, 268/270, 257-262



Modelling of tissue thermal injury formation process... 957

17. Mochnacki B., Majchrzak E., 2003, Sensitivity of the skin tissue on the activity of external
heat sources,CMES: Computer Modeling in Engineering and Sciences, 4, 3/4, 431-438

18. MochnackiB.,Piasecka-BelkhayatA., 2013,Numericalmodeling of skin tissue heating using
the interval finite difference method,MCB: Molecular and Cellular Biomechanics, 10, 3, 233-244

19. Oden J.T., Diller K.R., Bajaj C., Browne J.C., Hazle J., Babuska I., Bass J., Biduat
L., Demkowicz L., Elliott A., Feng Y., Fuentes D., Prudhomme S., Rylander M. N.,

Stafford R. J., ZhangY., 2007,Dynamic data-drivenfinite elementmodels for laser treatment
of cancer,Numerical Methods for Partial Differential Equations, 23, 904-922

20. Piasecka-BelkhayatA., 2011, Intervalboundaryelementmethod for transientdiffusionproblem
in two-layered domain, Journal of Theoretical and Applied Mechanics, 49, 1, 265-276

21. Torvi D.A., Dale J.D., 1994, A finite elementmodel of skin subjected to a flash fire, Journal of
Mechanical Engineering, 116, 250-255

Manuscript received April 2, 2014; accepted for print April 30, 2014