Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 57, 1, pp. 249-261, Warsaw 2019 DOI: 10.15632/jtam-pl.57.1.249 MODELLING OF BIOLOGICAL TISSUE DAMAGE PROCESS WITH APPLICATION OF INTERVAL ARITHMETIC Anna Korczak, Marek Jasiński Silesian University of Technology, Institute of Computational Mechanics and Engineering, Gliwice, Poland e-mail: anna.korczak@polsl.pl; marek.jasinski@polsl.pl In the paper, the numerical analysis of thermal processes proceeding in a 2D soft biological tissue subjected to laser irradiation is presented. The transient heat transfer is described by the bioheat transfer equation in Pennes formulation. The internal heat source resulting from the laser-tissue interaction based on the solution of the diffusion equation is taken into account. Thermophysical and optical parameters of the tissue are assumed as directed intervals numbers.At the stage of numerical realization. the interval finite differencemethod has been applied. In the final part of the paper, the results obtained are shown. Keywords: directed interval arithmetic, bioheat transfer, optical diffusion equation, Arrhe- nius scheme 1. Introduction One of the common problems inmathematical modelling of biological systems is the significant variation in the parameters that may result from, among others, individual characteristics. For some parameters, imprecisions in the estimation of their values may be a result of different measurement methods or specific features of tissues. The measurement is not often possible in a living organism. This may cause the values of individual parameters found in the literature to differ significantly from one another. It is also worth noting that tissue parameters may have different values depending on their degree of thermal damage. For example, the value of the scattering coefficient of tissue during the tissue damage process increases several times, what could be observed as “whitening” of the tissue (Jasiński, 2015). One of the most widely used mathematical tools to take into account uncertainties of para- meters are sensitivity analysis methods. By using them, very different problems were analysed, also in the field of heat transfer in the living organism (Jasiński, 2014; Kałuża et al., 2017; Mochnacki and Ciesielski, 2016). The approach using interval or fuzzy numbers is slightly less frequently used, although one can find works in which thermal processes in the human skin or eye are under consideration (Jankowska and Sypniewska-Kaminska, 2012; Mochnacki and Piasecka-Belkhayat, 2013; Piasecka-Belkhayat and Jasiński, 2011). The other initial-boundary value problems for partial differential equations were also considered in many fields of science (Gajda et al., 2000; Di Lizia et al., 2014; Nakao, 2017). One candistinguishbetween two types of interval arithmetic: directedand classical (Dawood, 2011; Hansen andWalster, 2004;Markov, 1995; Popova, 1994). Classical interval arithmetic was proposed in the 60s of the last century by Moore (1966). The two kinds of arithmetic differ mainly in the definition of an interval and the set of arithmetic operations. A classical interval x = [a,b] is defined as a set of real numbers that are between the lower bound a and upper bound b like x = {x ∈ R| a ¬ x ¬ b} while a directed interval is defined by only a pair of real numbers as x= [a,b] where a,b∈ R. Because of that, during successive calculations based on the classical interval arithmetic andmany computational methods, includingmethods based 250 A. Korczak,M. Jasiński on finite differences applied, the widths of intervals increase what leads to the so-called range overestimation problem (Hansen and Walster, 2004; Markov, 1995; Nakao, 2017). Using the directed interval arithmetic, it is possible to obtain the point (degenerated) interval [0,0] = 0 by subtraction of two identical intervals a−a= [0,0] and the point interval [1,1] = 1as the result of the division a/a=1,which is impossiblewhen applying the classical interval arithmetic. From a practical point of view, it is themost important one and causes the directed interval arithmetic to be popular in numerical applications (Dawood, 2011; Piasecka-Belkhayat andKorczak, 2016, 2017). Modelling of laser-biological tissue interactions requires appropriate mathematical descrip- tion. It is known that the scattering dominates over the absorption in soft tissues forwavelengths between 650 and 1,300nm (so-called biological window). Because of this, usually the radiative transport equation is taken into account (Dombrovsky andBaillis, 2010). There are several mo- difications in the discrete ordinates method or statistical Monte Carlo methods which are used to solve this equation (Banerjee and Sharma, 2010;Welch, 2011). In some cases, it is possible to approximate the light transport using the diffusion equation (Dombrovsky et al., 2013; Fasano et al., 2010; Jacques and Pogue, 2008). Modelling of the laser energydeposition is thefirst step in themodelling of physical processes proceeding in biological tissues subjected to a laser beam. Next, the temperature distribution must be calculated by making use of the bioheat transfer equation. The Pennes equation is the earliest one known but is probably still the most popular and widely used (Abraham and Sparrow, 2007; Majchrzak and Mochnacki, 2017; Paruch, 2014). The newest achievements in this field are based on the porous media theory (GDPL equation, generalized dual-phase lag equation) which takes into account the heterogeneous structure of biological tissue (Jasiński et al., 2016; Majchrzak andMochnacki, 2017; Majchrzak et al., 2015). The last step in the analysis of the tissue heating process is to estimate the degree of its destruction (Abraham and Sparrow, 2007; Henriques, 1947; Jasiński, 2018). The Arrhenius in- jury integral is the most frequently applied tool for this purpose, although other models (e.g. thermal dose) are also used (Mochnacki and Piasecka-Belkhayat, 2013). The Arrhenius scheme assumes the exponential dependence between temperature and the degree of tissue destruction. Furthermore, it refers only to the irreversible tissue damage, however, there are models which allow one to take into account the withdrawal of tissue injury in the case of temporary, small local increasing of temperature (Jasiński, 2014, 2018). The purpose of this paper is to analyse the phenomena occurring in the laser-treated soft tissuewherein thermophysical andoptical parameters aredefinedand treatedasdirected interval numbers. The analysis is based on the bioheat transfer equation in the Pennes formulation, whereas, to describe the light distribution in tissue the steady-state diffuse approximation is used. The degree of tissue destruction is also estimated by the use of the Arrhenius scheme. At the stage of numerical realisation, the interval finite difference method is used (Majchrzak and Mochnacki, 2016, 2017; Mochnacki and Suchy, 1995). 2. Formulation of the problem A transient heat transfer in biological tissue is described by the Pennes equation. The interval form of this equation can be expressed in the form (Abraham and Sparrow, 2007; Jankowska and Sypniewska-Kaminska, 2012; Jasiński, 2014; Paruch, 2014) x∈Ω : c ∂T ∂t =λ∇2T +Qperf +Qmet+Qlas (2.1) whereλ [Wm−1K−1] is the interval thermal conductivity, c [Jm−3K−1] is the interval volumetric specific heat, Qperf, Qmet and Qlas [Wm −3] are the interval internal heat sources containing Modelling of biological tissue damage process... 251 information connected with the perfusion, metabolism and laser irradiation, respectively, while T =T(x, t) [K] is the interval temperature. In this paper, the 2D domain of homogeneous biological tissue of a rectangular shape Ω subjected to the laser action is considered (Fig. 1). Fig. 1. The domain considered Equation (2.1) is supplemented by theRobin condition assumed on the external boundary of tissueΓ0, which is subjected to laser irradiationwhile on the remainingparts of the boundaryΓc the no-flux condition is accepted x∈Γ0 : q(x, t) =α(T −Tamb) x∈Γc : q(x, t) = 0 (2.2) where α [Wm−2K−1] is the convective heat transfer coefficient and Tamb is temperature of the surroundings. The initial distribution of temperature is also known x∈Ω, t=0 : T(x, t)=Tp (2.3) In the current work, the metabolic heat sourceQmet [Wm −3] is assumed as a constant interval parameter while the perfusion heat source is described by the formula Qperf(x, t)= cBw[TB −T(x, t)] (2.4) wherew [s−1] is the interval perfusion coefficient, cB [Jm −3K−1] is the volumetric specific heat of blood andTB corresponds to the arterial temperature (Abrahamand Sparrow, 2007;Mochnacki and Piasecka, 2013). The source function Qlas connected with the laser heating is defined as follows (Jasiński et al., 2016) Qlas(x, t) =µaφ(x)p(t) (2.5) where µa [m −1] is the interval absorption coefficient, φ(x) [Wm−2] is the interval total light fluence rate and p(t) is the function equal to 1 when the laser is on and equal to 0 when the laser is off. The total interval light fluence rateφ is the sumof the interval collimated partφc anddiffuse part φd (Banerjee and Sharma, 2010; Dombrovsky et al., 2013) φ(x)=φc(x)+φd(x) (2.6) The collimated fluence rate is given as (Jasiński et al., 2016) φc(x)=φ0exp ( − 2x22 r2 ) exp(−µ′tx1) (2.7) 252 A. Korczak,M. Jasiński whereφ0 [Wm −2] is the surface irradiance of laser, r is the radius of the laser beamandµ′t [m −1] is the interval attenuation coefficient defined as (Banerjee and Sharma, 2010) µ′t =µa+µ ′ s =µa+(1−g)µs (2.8) where µs and µ ′ s [m −1] are the interval scattering coefficient and the effective scattering coeffi- cient, respectively, while g is the anisotropy factor. To determine the interval diffuse fluence rate φd, the steady-state optical diffusion equation should be solved (Dombrovsky et al., 2013; Welch, 2011) x∈Ω : D∇2φd(x)−µaφd(x)+µ ′ sφc(x)= 0 (2.9) where D= 1 3[µa+(1−g)µs] = 1 3µ′t (2.10) is the interval diffusion coefficient. Equation (2.9) is supplemented by the boundary conditions on the boundaries Γ0 and Γc x∈Γ0, Γc : −Dn ·∇φd(x)= φd(x) 2 (2.11) where n is the outward unit normal vector. Damageofbiological tissue resulting fromtemperature elevation ismodelledbytheArrhenius injury integral, and its interval version considered in this paper is defined as (Abraham and Sparrow, 2007; Fasano et al., 2010; Henriques, 1947; Mochnacki and Piasecka-Belkhayat, 2013) Ψ(x, tF)= tF ∫ 0 P exp [ − E RT(x, t) ] dt (2.12) whereR [Jmole−1K−1] is the universal gas constant, E [Jmole−1] is the activation energy and P [s−1] is the pre-exponential factor while [0, tF ] is the considered time interval. The criterion for tissue necrosis is Ψ(x)­ 1. 3. Method of solution In this paper, both analyzed equations: the Pennes equation and the steady-state optical dif- fusion equation have been solved using the interval finite difference method (Mochnacki and Piasecka-Belkhayat, 2013). The information about the directed interval arithmetic are presen- ted in (Dawood, 2011; Piasecka-Belkhayat, 2011; Popova, 2011). In order to determine the function Qlas at the internal node (i,j) (cf. equation (2.5)), steady-state optical diffusion equation (2.9) is solved. The uniformdifferential grid of dimension 2n×2nwith spacing h/2 is used here (Fig. 2). Using such a differential grid, it is easier to take into account boundary conditions (2.11), because a part of the nodes is located exactly on the boundary Γ0 and Γc. In addition, one can distinguish common grid nodes for the temperature and diffuse fluence rate. The following differential quotients are used (∂2φd ∂x21 ) i,j = φdi+1,j −2φdi,j +φdi−1,j (h/2)2 (∂2φd ∂x22 ) i,j = φdi,j+1−2φdi,j +φdi,j−1 (h/2)2 (3.1) Modelling of biological tissue damage process... 253 Fig. 2. Five-point stencil and differential grid and then the approximate form of equation for the internal node (i,j) (i = 1,2, . . . ,2n− 1, j=1,2, . . . ,2n−1) is as follows φdi,j =C1Φdi,j +C2φci,j (3.2) where (superscripts “−” and “+”denote the beginning and the end of the interval, respectively) Φdi,j =φdi−1,j +φdi+1,j +φdi,j+1+φdi,j−1 =φ− di−1,j +φ − di+1,j +φ − di,j+1+φ − di,j−1,φ + di−1,j +φ + di+1,j +φ + di,j+1+φ + di,j−1] (3.3) while C1 = 4D 16D+h2µa = 4[D−,D+] 16[D−,D+]+h2[µ−a ,µ + a ] = [4D−,4D+] [16D−+h2µ−a ,16D++h2µ + a ] = [ 4D− 16D−+h2µ−a , 4D+ 16D++h2µ+a ] C2 = µ′sh 2 16D+µah 2 = h2[µ ′ − s ,µ ′+ s ] 16[D−,D+]+h2[µ−a ,µ + a ] = [h2µ′s,h 2µ ′+ s ] [16D−+h2µ−a ,16D++h2µ + a ] = [ h2µ ′ − s 16D−+h2µ−a , h2µ ′+ s 16D++h2µ+a ] (3.4) while for boundary nodes (cf. equation (2.11)) D φd1,j −φd0,j h/2 = 1 2 φd0,j → φd0,j =C3φd1,j −D φdn,j −φdn−1,j h/2 = 1 2 φdn,j → φdn,j =C3φdn−1,j D φdi,j −φdi,0 h/2 = 1 2 φdi,0 → φdi,0 =C3φdi,1 −D φdi,n−φdi,n−1 h/2 = 1 2 φdi,n → φdi,n =C3φdi,n−1 (3.5) where C3 = 4D 4D+h = 4[D−,D+] 4[D−,D+]+h2 = [4D−,4D+] [4D−+h2,4D++h2] = [ 4D− 4D−+h2 , 4D+ 4D++h2 ] (3.6) 254 A. Korczak,M. Jasiński It should be pointed out that the system of equations (3.2) is solved using the iterative method. The differential quotients approximating the derivatives appearing in the Pennes equation (2.1) are defined in a similar, as previously, manner (∂2T ∂x21 ) i,j = T i+1,j −2T i,j +T i−1,j h2 (∂2T ∂x22 ) i,j = T i,j+1−2T i,j +Ti,j−1 h2 ∂T ∂t = T f i,j −T f−1 i,j ∆t (3.7) where∆t denotes the time step while f−1 and f are the subsequent time levels. By introducing this formulas into (2.1), and after some mathematical operations for the internal nodes (i=1,2, . . . ,n, j=1,2, . . . ,n), one obtains T f i,j =A2T f−1 i,j +A1Υ f−1 i,j +A3 (3.8) where Υ f−1 i,j =T f−1 i+1,j +T f−1 i−1,j +T f−1 i,j+1+T f−1 i,j−1 = [ (T f−1 i+1,j) −+(T f−1 i−1,j) −+(T f−1 i,j+1) −+(T f−1 i,j−1) −, (T f−1 i+1,j) ++(T f−1 i−1,j) ++(T f−1 i,j+1) ++(T f−1 i,j−1) + ] (3.9) while A1 = λ∆t ch2 = ∆t[λ−,λ+] h2[c−,c+] = [∆tλ−,∆tλ+] [h2c−,h2c+] = [∆tλ− h2c− , ∆tλ+ h2c+ ] A2 =1−A1− wcB∆t c = [1−A−1 ,1−A + 1 ]− cB∆t[w −,w+] [c−,c+] = [1−A−1 ,1−A + 1 ] − [cB∆tw − c− , cB∆tw + c+ ] = [ 1−A−1 − cB∆tw − c− ,1−A+1 − cB∆tw + c+ ] A3 = ∆t c (wcBTB +Qmet+Qlas)= ∆t [c−,c+] (cBTB[w −,w+]+ [Q−met,Q − met]+ [Q − las ,Q− las ]) = [∆t(cBTBw −+Q−met+Q − las ) c− , ∆t(cBTBw ++Q+met+Q + las ) c+ ] (3.10) The “boundary” nodes are located at the distance 0.5h from the boundary of the domain. This approach gives a better approximation of the Neumann and Robin boundary conditions, but the final form of the interval FDMequation for the boundary nodes are obtained in similar way. More details of this approach can be found in (Jasiński et al., 2016; Majchrzak andMochnacki, 2016, 2017). On the basis of equation (3.8), the temperature at the node (i,j) for the time level can be foundon the assumption that the stability condition for an explicit differential scheme is fulfilled (Mochnacki and Suchy, 1995). 4. Results of computations At the stage of numerical computations, a 2Dhomogeneous tissue domain of size 4×4cmduring laser irradiation (Fig. 1) has been considered. The uniform differential grid with 40×40 nodes is introduced (Fig. 2). As has been already mentioned, the values of the effective scattering coefficient are different for the native and thermally-damaged (denaturated) tissue, so the two Modelling of biological tissue damage process... 255 simulations have been conducted (for µ′s = µ ′ snat and µ ′ s = µ ′ sden). The optical (µa,µ ′ s) and thermophysical (λ,c,w,Qmet) parameters of the tissue are assumed as the directed interval numbers in form (Dombrovsky et al., 2013; Popova, 2011) p= [p−0.025p,p+0.025p] (4.1) where p denotes the parameter. The following data have been assuming in calculations: λ = 0.609Wm−1K−1, c = 4.18MJm−3K−1, w = 0.00125s−1, µa = 40m −1, µ′snat = 1000m −1, µ′sden = 4000m −1, Qmet = 245Wm −3, cB = 3.9962MJm −3K−1, TB = 37 ◦C, P = 3.1 · 1098 s−1, E = 6.27 · 105Jmole−1, R = 8.314Jmole−1K−1, φ0 = 3 · 10 5Wcm−2, d = 2mm, texp = 25s, α=10Wm−2K−1, Tamb =20 ◦C, Tp =37 ◦C,∆t=1s. Fig. 3. Distribution of the interval collimated fluence rate φ c (x2 =0) Fig. 4. Distribution of the interval diffuse fluence rate φ d (x2 =0) Figures 3-5 are associated with the light fluence distribution in the domain considered. Figure 3 presents the distribution of the interval collimated fluence rate φc calculated on the 256 A. Korczak,M. Jasiński basis of equation (2.7) while in Fig. 4 the distribution of the interval diffuse fluence rate φd resulting from steady-state optical diffusion equation (2.9) is presented. The distribution of the interval diffuse fluence rate φd is also presented in Fig. 5. It is visible that the area of scattering is larger in the case of calculation with µ′snat although the values of φd are lower than in the case of calculation with µ′sden. Fig. 5. Distribution of the interval diffuse fluence rate φ d [kWm−2] The next figure is associated with the tissue temperature. In Fig. 6, the interval tissue temperature history at the nodeN0(0,0) (Fig. 1) obtained for the effective scattering coefficient of native and denaturated tissue is presented. As can be seen, wider temperature intervals are obtained in the case of denaturated tissue. Fig. 6. History of interval tissue temperature at the nodeN0 Modelling of biological tissue damage process... 257 Fig. 7. History of the interval Arrhenius integral at the nodeN0 Fig. 8. Distributions of the interval Arrhenius injury integral for native tissue (µ′ s =µ′ snat ) 258 A. Korczak,M. Jasiński InFig. 7, the history of the interval Arrhenius integral at the nodeN0 is shown.As expected, due to thewider intervals obtained in the calculations withµ′sden, wider intervals for the interval Arrhenius integral are also obtained for this variant of calculations. It is visible that tissue destruction occurs first in the case of native tissue – the injury integral reaches the necrosis criterion (Ψ ­ 1) in the time interval [14,15]s while for calculation with µ′sden this criterion is reached in the time interval [17,21] s. The results associatedwith the tissue destruction are also presented inFigs. 8 and 9. In both figures, the interval injury integral distributions for selected time steps are presented. Thewhite zone in these figures refers to the values of the injury integral below 0.01 (thermally untouched tissue), the grey zone refers to the values 0.01 < Ψ < 1, so it is a partial damage area, and the black zone illustrates the area in which the Arrhenius integral achieved the criterion of tissue necrosis. In both cases, the coagulation zones obtained for the lower and upper bounds of intervals are slightly different, however, in the case of destructed tissue the differences are bigger. Fig. 9. Distributions of the interval Arrhenius injury integral for denaturated tissue (µ′ s =µ′ sden ) 5. Discussion As has beenmentioned before, the reason for the application of the directed interval arithmetic was to take into account the inaccuracies of estimation of biological tissue parameters and their effect on thermal damage to the tissue. It can be seen that the results associated with the tissue Modelling of biological tissue damage process... 259 damage are quite reasonable for the assumed width of interval parameters. The differences between the ranges of coagulation zones (Figs. 8 and 9) obtained for both assumed values of the interval effective scattering coefficient µ′snat and µ ′ sden are clearly visible. On the other hand, the differences between the lower and upper bounds of the interval Arrhenius integral for both cases of µ′s are not very big. Of course, thermal damage can also expand after the impulse causing the temperature ele- vation ceases, what happens also in the analysed cases for t > texp, and the thermal injury is fully formed in less than 40s. Application of the Arrhenius integral also allows one to estimate the time after which the given point of the domain considered is thermally damaged, i.e. the valueΨ reaches the necrosis criterion. The results of this type are shown in Fig. 7. These results allow one to estimate the time of tissue damage in a few second intervals – at the nodeN0 for the calculation with µ ′ snat [14,15] s while for the calculation with µ′sden [17,21] s. Such width of intervals can be accepted as useful information from a practical point of view. It should be pointed out that, in some cases, it is justified to assume a more restrictive necrosis criterion, i.e. Ψ(x) ­ 4.6 which corresponds to 99% of dead cells (e.g. laser cancer therapy). This criterion is also fulfilled at the nodeN0 for the assumed width of the parameter intervals – for calculation with µ′snat [15,17]s, for calculation with µ ′ sden [19,25] s. These width of intervals are a still reasonable outcome. One should also note that an increase in the intervals for optical and thermophysical para- meters leads to an increase in temperature intervals, which are the basic value in calculations of the interval Arrhenius integral. This, in turn, may lead to a situation in which the lower bound of the interval injury integral for a selected point of the domain considered would be below the threshold of necrosis, whereas the upper bound would exceed this threshold. 6. Conclusions In the paper, the analysis of the thermal damage of biological tissue subjected to laser impulse has been presented, whereas, the optical and thermophysical parameters of tissue have been treated as the directed interval numbers. That concerned two optical parameters: absorption coefficient µa, effective scattering coefficient of tissue µ ′ s, and four thermophysical parameters: thermal conductivity λ, volumetric specific heat c, perfusion coefficient w and metabolic heat sourceQmet. The combined effect of those parameters has been considered, although it is known that not all of these parameters (or actual changes in parameters values) have the same influence on the temperature level and, as a result, on the estimated value of tissue damage. It has been described in numerous works related to e.g. sensitivity analysis (Majchrzak and Mochnacki, 2017;Mochnacki andCiesielski, 2016). Of course, interval analysis for the individual parameters of tissue is also possible to perform. The directed interval arithmetic has already been effectively applied to the modelling of bioheat transfer problems.Most of the works were based on the Pennes equation, although it is possible to use a different bioheat transfer equation (Mochnacki and Piasecka-Belkhayat, 2013). Of particular interest here is the use of the GDPL equation based on the theory of porous bodies which binds the process of heat transfer with the inner tissue structure (Majchrzak and Mochnacki, 2017; Majchrzak et al., 2015). The modelling of the tissue thermal damage process also includes a number of other pro- blems, e.g. methods of taking into account changes in parameter values caused by the degree of tissue damage (Jasiński, 2015, 2018). Further work related to the application of directed inter- val arithmetic in this area could concern these issues or the problems related with the use of interval numbers in other algorithms related to thermal damage (Jasiński, 2015;Mochnacki and 260 A. Korczak,M. Jasiński Piasecka, 2013). Because the interval arithmetic is a method developed by mathematicians in order to put bounds on rounding errors andmeasurement errors in mathematical computation, it could also be a useful approach to the modelling of various problems of bioheat transfer. Acknowledgements The research has been funded from the projects of Silesian University of Technology, Faculty of Mechanical Engineering. References 1. Abraham J.P., Sparrow E.M., 2007, A thermal-ablation bioheat model including liquid-to- -vapor phase change, pressure- and necrosis-dependent perfusion, andmoisture-dependent proper- ties, International Journal of Heat and Mass Transfer, 50, 13-14, 2537-2544 2. Banerjee S., Sharma S.K., 2010, Use of Monte Carlo simulations for propagation of light in biomedical tissues,Applied Optics, 49, 4152-4159 3. Dawood H., 2011,Theories of Interval Arithmetic: Mathematical Foundations and Applications, LAP LAMBERTAcademic Publishing GmbH&Co., Germany 4. Dombrovsky L.A., Baillis D., 2010,Thermal Radiation in Disperse Systems: An Engineering Approach, Begell House, NewYork 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. Fasano A., Hömberg D., Naumov D., 2010, On a mathematical model for laser-induced ther- motherapy,Applied Mathematical Modelling, 34, 12, 3831-3840 7. Gajda K., Marciniak A., Szyszka B., 2000, Three- and four-stage implicit interval methods of Runge-Kutta type,Computational Methods in Science and Technology, 6, 41-59 8. Hansen E.,WalsterG.W., 2004,Global Optimization Using Interval Analysis, Marcell Dekker, NewYork 9. Henriques F.C., 1947, Studies of thermal injuries, V. The predictability and the significance of thermally induced rate process leading to irreversible epidermal injury, Journal of Pathology, 23, 489-502 10. Jacques S.L., Pogue B.W., 2008, Tutorial on diffuse light transport, Journal of Biomedical Optics, 13, 4, 1-19 11. Jankowska M.A., Sypniewska-Kaminska G., 2012, Interval finite difference method for bio- heat transfer problem given by the Pennes equation with uncertain parameters, Mechanics and Control, 31, 2, 77-84 12. JasińskiM., 2014,Modelling of tissue thermal injury formation process with application of direct sensitivity method, Journal of Theoretical and Applied Mechanics, 52, 947-957 13. Jasiński M., 2015, Modelling of thermal damage in laser irradiated tissue, Journal of Applied Mathematics and Computational Mechanics, 14, 67-78 14. Jasiński M., 2018, Numerical analysis of soft tissue damage process caused by laser action,AIP Conference Proceedings, 060002, 1922 15. JasińskiM.,MajchrzakE., Turchan L., 2016,Numerical analysis of the interactions between laser and soft tissues using generalized dual-phase lag model, Applied Mathematic Modelling, 40, 2, 750-762 16. Kałuża G., Majchrzak E., Turchan L., 2017, Sensitivity analysis of temperature field in the heated soft tissue with respect to the perturbations of porosity, Applied Mathematical Modelling, 49, 498-513 Modelling of biological tissue damage process... 261 17. Di Lizia P., Armellin R., Bernelli-Zazzera F., Berz M., 2014, High order optimal control of space trajectories with uncertain boundary conditions,Acta Astronautica, 93, 217-229 18. Majchrzak E., Mochnacki B., 2016, Dual-phase lag equation. Stability conditions of a nume- rical algorithm based on the explicit scheme of the finite difference method, Journal of Applied Mathematics and Computational Mechanics, 15, 89-96 19. Majchrzak E., Mochnacki B., 2017, Implicit scheme of the finite difference method for 1D dual-phase lag equation, Journal of Applied Mathematics and Computational Mechanics, 16, 3, 37-46 20. 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 21. Markov S.M., 1995, On directed interval arithmetic and its applications, Journal of Universal Computer Science, 1, 514-526 22. Mochnacki B., Ciesielski M., 2016, Sensitivity of transient temperature field in domain of forearm insulated by protective clothing with respect to perturbations of external boundary heat flux,Bulletin of the Polish Academy of Sciences – Technical Sciences, 64, 3 23. MochnackiB.,Piasecka-BelkhayatA., 2013,Numericalmodeling of skin tissue heating using the interval finite difference method,Molecular and Cellular Biomechanics, 10, 3, 233-244 24. Mochnacki B., Suchy J.S., 1995, Numerical Methods in Computations of Foundry Processes, PFTA, Cracow 25. Moore R.E., 1966, Interval Analysis, Prentice-Hall, EnglewoodCliffs 26. Nakao M., 2017,On the initial-boundary value problem for some quasilinear parabolic equations of divergence form, Journal of Differential Equations, 263, 8565-8580 27. Paruch M., 2014, Hyperthermia process control induced by the electric field in order to destroy cancer,Acta of Bioengineering and Biomechanics, 16, 4, 123-130 28. Piasecka-Belkhayat A., 2011, Interval Boundary Element Method for Imprecisely Defined Unsteady Heat Transfer Problems (in Polish), Silesian University of Technology, Gliwice 29. Piasecka-Belkhayat A., Jasiński M., 2011,Modelling of UV laser irradiation of anterior part of human eye with interval optic, [In:] Evolutionary and Deterministic Methods for Design, Opti- mization and Control. Applications, Burczyński T., Périaux J. (Eds.), CIMNE,Barcelona, 316-321 30. Piasecka-Belkhayat A., Korczak A., 2016, Numerical modelling of the transient heat trans- port in a thin gold film using the fuzzy lattice Boltzmann method with alpha-cuts, Journal of Applied Mathematics and Computational Mechanics, 15, 1, 123-135 31. Piasecka-BelkhayatA., KorczakA., 2017,Modeling of thermal processes proceeding in a 1D domain of crystalline solids using the lattice Boltzmannmethod with an interval source function, Journal of Theoretical and Applied Mechanics, 55, 1, 167-175 32. Popova E.D., 1994, Extended interval arithmetic in IEEE floating-point environment, Interval Computations, 4, 100-129 33. PopovaE.D., 2011,Multiplication distributivity of proper and improper intervals,Reliable Com- puting, 7, 129-140 34. Welch A.J., 2011, Optical-Thermal Response of Laser Irradiated Tissue, M.J.C. van Gemert (Eds.), 2nd edit., Springer Manuscript received May 10, 2018; accepted for print October 25, 2018