# Modelling of thermal damage process in soft tissue subjected to laser irradiation

### Marek Jasiński

Journal of Applied Mathematics and Computational Mechanics |
Download Full Text |

MODELLING OF THERMAL DAMAGE PROCESS IN SOFT TISSUE SUBJECTED TO LASER IRRADIATION

Marek Jasiński

Institute of Computational Mechanics and
Engineering

Silesian University of Technology

Gliwice, Poland

marek.jasinski@polsl.pl

Received: 12 April 2018; Accepted: 18 June
2018

**Abstract.** The
numerical analysis of thermal damage process proceeding in biological
tissue during laser irradiation is presented. Heat transfer in the tissue is
assumed to be transient and two-dimensional. The internal heat source resulting
from the laser irradiation based on the solution of the diffusion equation is
taken into account. The tissue is regarded as a homogeneous domain with
perfusion coefficient and effective scattering coefficient treated as dependent
on tissue injury. At the stage of numerical realization, the boundary element
method and the finite difference method have been used. In the final part of
the
paper the results of computations are shown.

*MSC 2010:** 65M06, 65M38, 80A20*

*Keywords: **bioheat transfer,
optical diffusion equation, Arrhenius scheme, boundary
element method, finite difference method*

1. Introduction

It is known that biological tissues are characterized by a strong scattering and weak absorption in the so-called therapeutic window (wavelengths 650¸1300 nm) [1-4]. Because of this, to describe the light propagation in biological tissues the different mathematical models can be taken into account. One of them is the transport theory which concerns the transport of light through scattering and absorbing media. In this case the radiative transport equation should be considered [4-6]. To solve this equation, different approaches are used: the modifications of discrete ordinates method, the statistical Monte Carlo method or the optical diffusion equation [2, 5, 7].

Interactions between tissue and the laser beam often lead to the temperature elevation that can cause irreversible damage of the tissue. This in turn could cause the alteration of thermophysical and optical properties of tissue, e.g the perfusion coefficient is often treated as the main indicator of tissue injury. In addition, an increase of scattering due to tissue denaturation has a visual effect of tissue “whitening”. Consequently, parameters applied in mathematical models of heat transfer in biological tissue domain can be regarded as tissue damage-dependent [8-11].

The most popular model of the tissue damage process is the Arrhenius injury integral [12-16]. This approach basically refers only to the irreversible tissue damage, however, there are models that allow one to take into account the withdrawal of tissue injury in the case of temporary, small local increasing of temperature [17].

The purpose of this paper is to analyze the phenomena occurring in the laser-treated biological tissue. The analysis is based on the Pennes bioheat transfer equation which is still the most frequent used model to determine the temperature distribution in biological tissues [9, 10, 13, 18], while the light distribution in biological tissue is estimated on the base of the optical diffusion equation [4, 5, 16]. The degree of tissue damage is calculated by using of Arrhenius scheme with TTIW algorithm and two parameters of the tissue are assumed to be dependent on the injury integral (perfusion coefficient and the effective scattering coefficient).

2. Governing equations

The transient heat transfer in the 2D homogeneous biological tissue domain of rectangular shape W (Fig. 1) is described by the Pennes bioheat transfer equation with adequate boundary-initial conditions [11, 18]

where λ [Wm^{-1}K^{-1}] is the thermal conductivity, *c*
[Jm^{-3}K^{-1}] is the volumetric
specific heat, *T = T*(**x**, *t*) is the
temperature while denotes its time derivative.
The components *Q _{perf}*,

_{ }

*Q*and

_{las}*Q*[Wm

_{met}^{-3}] are the internal source functions containing information connected with the perfusion, metabolism and laser irradiation. In boundary-initial conditions, a [Wm

^{-2}K

^{-1}] is the convective heat transfer coefficient,

*T*is the temperature of surroundings while

_{amb}*T*denotes the initial distribution of temperature.

_{p}The metabolic heat source *Q _{met}*
[Wm

^{–3}] is assumed as a constant value while the perfusion heat source is described by the formula

where *c _{B}* [Jm

^{–3}K

^{–1}] is the volumetric specific heat of blood,

*T*corresponds to the arterial temperature while

_{B}*w*(Y) [s

^{–1}] is the perfusion coefficient defined as [8, 9, 17]

where *w*_{0} is the initial
perfusion coefficient while Y denotes the Arrhenius injury integral, in form [9, 10, 12]

where *R *[J mole^{–1}K^{–1}]
is the universal gas constant, *E* [J mole^{–1}] is the
activation energy, *P* [s^{-1}] is the pre-exponential factor. A value of
integral Y(**x**) = 1 corresponds to a 63% probability of cell
death at a specific point **x**, while Y(**x**) = 4.6 corresponds to a 99% probability of cell death at this point. Both values are
used as the necrosis criteria.

Fig. 1. The domain considered

According to the values of formula , the values of coefficients for the interval from 0 to 0.1 in equation respond to the initial increase of the perfusion coefficient caused by vasodilatation, while values of coefficients for the interval 0.1 to 1 correspond to the decrease of the perfusion during the damage process of the vascular network [9].

The main assumption of the Arrhenius formula is that the damage of tissue is irreversible, so even in the case of very little rise and lowering of temperature the tissue remain damaged. On the other hand, at the initial tissue heating, when the temperature is moderate (i.e. between 37°C and 45¸55°C), the blood vessels in the tissue become dilated without being thermally damaged. Because of this, the TTIW algorithm (thermal tissue injury withdrawal algorithm) has been applied which allows to take into account the possibility of withdrawal of tissue injury when the thermal impulse is ceased [17].

The TTIW algorithm requires the assumption of
a certain value Ψ* _{rec}* which
is defined as recovery threshold. The withdrawal of tissue damage is possible
only for those points of the domain considered in which the value of injury
coefficient Ψ is below Ψ

*. If the tissue damage at the point*

_{rec}**x**Î W achieves the value greater than Ψ

*then the injury becomes irreversible, so, it will be calculated on the base of the Arrhenius formula . The full description of the TTIW algorithm could be found in [17].*

_{rec}The thermal damage of the tissue also affects the optical properties of the tissue. During the process of coagulation, the changes in tissues lead to higher scattering while the absorption remains the same, thus in the current work, the effective scattering coefficient of tissue is described as [10]

where m¢* _{s nat}* and m¢

*[cm*

_{s den}^{-1}] denote the effective scattering coefficient for native and destructed (denaturated) tissue, respectively.

The source function *Q _{las}*
(c.f. equation ) associated with the laser heating is defined as follows [2]

where m* _{a}* [m

^{–1}] is the absorption coefficient, f

*and f*

_{c}*are the collimated and diffuse parts of the total light fluence rate, respectively while*

_{d}*p*(

*t*) is the function equal to 1 when the laser is

*on*and equal to 0 when the laser is

*off.*

The collimated fluence rate part f* _{c}* is determined on the basis of the Beer-
-Lambert law, namely [2, 11, 16]

where f_{0} [Wm^{–2}]
is the surface irradiance of laser, *d* is the diameter of laser beam and m¢* _{t}* [m

^{–1}] is the attenuation coefficient given as

To determine the diffuse fluence rate f* _{d}* the optical diffusion equation should be solved [1-3]

where **n** is
the outward unit normal vector.

3. Method of solution

At the stage of numerical realization, the 1^{st}
scheme of the boundary element method (BEM) for 2D transient heat diffusion has
been used while the optical
diffusion equation has been solved by the finite difference method (FDM).

For the transient heat diffusion problem, for
the time grid with constant step D*t*, the boundary
integral equation corresponding to transition *t ^{ f}*

^{–1}®

*t*is of the form

^{ f}where *Q _{V}* denotes the sum of
the internal heat function associated with the perfusion, metabolism and laser
irradiation (c.f. equation ),

*T*

^{*}is the fundamental solution

where *r* is the distance from the point
under consideration **x** to the observation point **x**, *a* =
l/*c*,
while

and *B*(**x**) is the coefficient
from the interval (0, 1).

In this paper, the constant boundary element has been used. Details concerning numerical realization of the BEM can be found, among others, in [11, 19-21].

In order to determine the source function *Q _{las}*
at the internal nodes (c.f. equation ), at each time step additionally the optical diffusion equation must be solved. As it was mentioned previously, it is done by using
the finite difference method. The global and local numeration of the nodes is
shown in Figure 2.

The difference equation for the central node of stencil can be written in the form

where *R*_{0e} are
resistances between the central node and remaining nodes of the stencil -
Figure 2.

More details concerning numerical realization of the FDM, can be found in [16, 22, 23].

Fig. 2. Five-point stencil and discretization

4. Results of computations

The aim of the research was to analyse the destructive changes in the tissue domain of the size of 4´4 cm during laser irradiation (c.f. Fig. 1). The interior of the domain has been divided into 1600 internal constant cells while the external boundary into 160 constant elements.

The most commonly used types of laser
impulses in medical treatments are multiple impulses of duration ranging from a
few to several dozen seconds with appropriate pauses between subsequent
impulses. The purpose of this is not to overheat, and thus not to cause thermal
damage to healthy tissue adjacent to
the treated area. Thermography techniques are also used for similar reasons.
They
allow one to maintain the tissue surface temperature in a certain range and to
control the laser action during the medical procedure [1, 5, 10, 13].
Because of this, in current paper, the two cases of laser irradiation were
taken under consideration. In example 1, the tissue is subjected to multiple
laser impulse (100 seconds *on* and 100 seconds *off*) while in
example 2 the laser action depends on the tissue surface temperature**,**
meaning the laser is *on* when the temperature drop below *T _{ctrl}*
– D

*T*and

_{ctrl}*off*when the temperature reaches above

*T*+ D

_{ctrl}*T*.

_{ctrl}In computations, the following values of
parameters have been assumed: l _{ }= _{ }0.609 _{ }Wm^{-1}K^{-1}, *c* _{ }= _{ }4.18 _{ }MJm^{-3}K^{-1}, *w*_{0} _{ }= _{ }0.00125 _{ }s^{–1},
m_{a}_{ }= _{ }0.4 _{ }cm^{–1},
m¢_{s nat}_{ }= _{ }10 _{ }cm^{–1}, m¢_{s den}_{ }= _{ }40 _{ }cm^{–1}, *Q _{met}*

_{ }=

_{ }245

_{ }Wm

^{–3},

*c*

_{B}_{ }=

_{ }3.9962

_{ }MJm

^{-3}K

^{-1},

*T*

_{B}_{ }=

_{ }37°C,

*P*

_{ }=

_{ }3.1×10

^{98}

_{ }s

^{–1},

*E*

_{ }=

_{ }6.27×10

^{5}

_{ }Jmole

^{–1},

*R*

_{ }=

_{ }8.314

_{ }Jmole

^{–1}K

^{–1}, Y

*= 0.05, ϕ*

_{rec}_{0}= 30 Wcm

^{–2},

*d*= 2 mm, a = 10 Wm

^{–2}K

^{–1},

*T*= 20°C,

_{amb}*T*= 37°C,

_{p}*T*± D

_{ctrl}*T*= 70±3°C, D

_{ctrl}*t*= 1 s [10, 11, 17].

The optical parameters of tissue (m* _{a}*, m¢

*and m¢*

_{s nat}*), which were assumed in simulations, are typical for near-IR irradiation on soft tissue like e.g. Nd:YAG laser of 1064 nm which is used for prostate coagulation.*

_{s den}The results containing information about
courses in time are presented in points *N*_{0}(0,0), *N*_{1}(0.0045,0),
*N*_{2}(0.0065,0), *N*_{3}(0.0085,0).

In Figure 3, distribution of the diffuse
fluence rate obtained on the basis of the optical diffusion equation is shown. It is visible that initially the
area of scattering is larger, although the values of f* _{d}* are lower than they are in the end of the simulation. The results
are presented for example 1 only, but it should be pointed out that in example
2 they were almost the same.

Fig. 3. Distribution
of diffuse fluence rate ϕ* _{d}* [kWm

^{–2}] for 0 and 500 s (example 1)

Figures 4 and 5 are associated with tissue
temperature. In Figure 4, the temperature history for the point *N*_{0}
and both examples analyzed are presented. For multiple laser impulse (example
1), it is visible that for each of the subsequent laser’s cycle (*on*/*off*),
the reached temperatures are higher and higher. Lower temperatures obtained for
the example 2 are also visible in Figure 5.

Fig. 4. Tissue
temperature history at the node *N*_{0} (examples 1 and 2)

Fig. 5. Temperature distribution [°C], time 100, 300, 500 s (up: example 1, down: example 2)

Fig. 6.
Arrhenius injury integral history at nodes *N*_{1}, *N*_{2},
*N*_{3} (examples 1 and 2)

In Figure 6, the time courses of the Arrhenius integral. The
decreases of the
Arrhenius integral value resulting from the application of the TTIW algorithm
were noticed: for example 1 at point *N*_{2} while for example 2
at point *N*_{1}. In example 2, they are not big, but they are
still visible.

Figure 7 shows the distribution of the Arrhenius injury integral for times 100, 300 and 500 s for both examples. The differences are visible mainly at the early stage of the tissue damage process.

Fig. 7. Arrhenius injury integral distribution, time 100, 300, 500 s (up: example 1, down: example 2)

Obviously, all changes in the value of the Arrhenius integral have reflection in changes of tissue damage-dependent parameters assumed in the current model, i.e. perfusion coefficient and effective scattering coefficient. The distribution of these parameters is presented in the next two figures.

In Figure 8, the grey zone refers to the zone in which the value of the perfusion coefficient is lowered while the check zone is the so-called hyperemic ring - the area in which the value of perfusion is raised. As it can be seen, this zone is almost adjacent to the coagulation zone. In Figure 9, the distribution of the effective scattering coefficient is presented.

In Figure 10, the information about expansion of tissue damage is presented. The number of elements on the scale of the figure means that tissue injury is treated as a sum of element on which the injury integral is above 0.01, so this value could be treated as the border of thermally untouched tissue.

Fig. 8. Perfusion coefficient
distribution *w*(Y) ´ 1000 [s^{–1}], time 100, 300, 500 s
(example 2)

Fig. 9. Effective scattering
coefficient m¢* _{s}*(Y) [cm

^{–1}], time 100, 300, 500 s (example 2)

Fig. 10. Expansion of the thermal injury

6. Conclusion

The heat transfer process in the tissue subjected to laser irradiation is analyzed. The analysis is based on the bioheat transfer equation in the Pennes formulation while for calculation of light distribution in tissue, the Beer-Lambert law (collimated part of fluence rate) and optical diffusion equation (diffuse part of fluence rate) are used. The perfusion coefficient and effective scattering coefficient are assumed to be dependent on the tissue damage which was estimated on the base of Arrhenius injury integral with the TTIW algorithm. Because the calculation in this way is closer to real conditions of the tissue destruction process during laser-tissue interaction, the estimation of the total size of tissue damage is more precise.

Comparing the process of increasing of tissue damage, it is visible that for the example 2 the expansion of the lesion occurs slightly more evenly. The final areas of tissue damage are almost the same in both cases (c.f. Figure 10), however in the example 2 it was achieved at a lower local temperature value (c.f. Figures 4 and 5).

It should be pointed out that for the description of light propagation in tissue, different models could also be taken into account e.g. the Monte Carlo approach [7]. It is also possible to use another equation of bioheat transfer such as Cattaneo- -Vernotte equation or dual phase lag equation [4, 5, 14, 15, 20-22, 24]. Furthermore, due to individual characteristics of biological objects, the various values of thermophysical and optical parameters should be considered. For this purpose, for example, the sensitivity analysis can be used.

Acknowledgements

*The research is funded from the projects
Silesian University of Technology, Faculty of Mechanical Engineering.*

References

[1]
Welch, A.J., & van Gemert, M.J.C. (Eds.) (2011). *Optical-thermal Response of Laser Irradiated
Tissue*. Springer.

[2]
Niemz, M.H. (2007). *Laser-tissue Interaction*.
Berlin, Heidelberg, New York: Springer-Verlag.

[3]
Tuchin, V.V. (2007). *Tissue Optics: Light
Scattering Methods and Instruments for Medical
Diagnosis *(Vol. PM166). Bellingham: SPIE Press.

[4]
Dombrovsky, L.A. (2012). The use of
transport approximation and diffusion-based models
in radiative transfer calculations. *Computational Thermal Sciences*, *4*(4),
297-315.

[5]
Dombrovsky, L.A.,
Randrianalisoa, J.H., Lipinski, W., & Timchenko, V. (2013).
Simplified
approaches to radiative transfer simulations in laser induced hyperthermia of
superficial tumors. *Computational Thermal Sciences*, *5*(6),
521-530.

[6]
Jacques, S.L., & Pogue, B.W. (2008). Tutorial on diffuse light transport. *Journal of Biomedical Optics*, *13*(4), 1-19.

[7]
Banerjee, S., & Sharma, S.K.
(2010). Use of Monte Carlo simulations for propagation of light
in biomedical tissues. *Applied Optics*, *49*, 4152-4159.

[8]
Jasiński, M. (2011). Numerical modeling of
tissue heating with application of sensitivity
methods. *Mechanika 2011: Proceedings of the 16th International Conference*,
137-142.

[9]
Abraham, J.P., & Sparrow, E.M.
(2007). A thermal-ablation bioheat model including liquid-to-vapor phase
change, pressure- and necrosis-dependent perfusion, and moisture-dependent
properties. *International Journal of Heat and Mass Transfer*, *50*(13-14),
2537-2544.

[10] Glenn, T.N., Rastregar, S., & Jacques, S.L. (1996). Finite element analysis of temperature controlled coagulation in
laser irradiated tissue. *IEEE Transactions on**
Biomedical **Engineering**,*
*43*(1), 79-86.

[11] Jasiński, M. (2018). Numerical analysis of soft tissue damage
process caused by laser action. *AIP Conference Proceedings*, 1922, 060002.

[12] 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.

[13] Fasano, A., Hömberg, D., & Naumov, D. (2010). On a mathematical model for laser-induced
thermotherapy. *Applied Mathematical Modelling*, *34*(12), 3831-3840.

[14] Narasimhan, A., & Sadasivam, S. (2013). Non-Fourier
bio heat transfer modelling of thermal damage during retinal laser irradiation.
*International Journal of Heat and Mass Transfer*, *60*, 591-597.

[15] Mochnacki, B., & Majchrzak, E. (2017). Numerical model
of thermal interactions between
cylindrical cryoprobe and biological tissue using the dual-phase lag equation. *International
Journal of Heat and Mass Transfer*, *108*(1-10), 2017.

[16] Jasiński, M.,
Majchrzak, E., & Turchan, L. (2016). Numerical
analysis of the interactions
between laser and soft tissues using generalized dual-phase lag model. *Applied Mathematical Modelling*, *40*(2),
750-762.

[17] Jasiński, M. (2015). Modelling of thermal damage in laser
irradiated tissue. *Journal of Applied Mathematics and Computational
Mechanics*, *14*, 67-78.

[18] Paruch, M. (2017), Identification of the cancer ablation
parameters during RF hyperthermia
using gradient, evolutionary and hybrid algorithms. *International Journal of
Numerical Methods for Heat & Fluid Flow*, *27*, 674-697.

[19] Brebia, C.A., & Dominquez, J. (1992). *Boundary Elements, an Introductory Course *(*Computational Mechanics Publications*). London: McGraw-Hill
Book Company.

[20] Majchrzak, E., & Turchan, L. (2015). The general
boundary element method for 3D dual-phase lag model of bioheat transfer. *Engineering Analysis with Boundary Elements*, *50*,
76-82.

[21] Majchrzak, E. (2010). Numerical solution of dual phase lag
model of bioheat transfer using
the general boundary element method. *CMES*: *Computer **Modeling* *in* *Engineering* *&
Sciences*, *69*(1), 43-60.

[22] Majchrzak, E., & Mochnacki, B. (2016). Dual-phase lag
equation. Stability conditions of a numerical algorithm based on the
explicit scheme of the finite difference method. *Journal of Applied
Mathematics and Computational Mechanics*, *15*, 89-96.

[23] Majchrzak, E.,
Turchan, L., & Dziatkiewicz, J. (2015). Modeling
of skin tissue heating using
the generalized dual-phase lag equation. *Archives of
Mechanics*, *67*(6), 417-437.

[24] Xu, F., Seffen, K.A., & Lu, T.J. (2008). Non-Fourier analysis of skin biothermomechanics.
*International Journal of Heat and Mass Transfer*, *51*(9-10),
2237-2259.