Identification of Laser Intensity Assuring the Destruction of Target Region of Biological Tissue Using the Gradient Method and Generalized Dual-Phase Lag Equation

Temperature distribution in biological tissue is described by dual-phase lag model supplemented by appropriate boundary and initial conditions. Laser–tissue interactions are taken into account in the source function occurring in this model. The problem is solved using the finite difference method. Next, the Arrhenius integral which is a measure of tissue destruction is calculated. The inverse problem formulated here is to estimate the laser intensity assuring the destruction of this tissue region for which the Arrhenius integral exceeds the critical value.


Introduction
The goal of biological tissue artificial heating is to destroy the diseased fragments and to minimize the damage of surrounding healthy tissue. Artificial hyperthermia treatments are performed using, among others, the electromagnetic field or the laser action. In this paper, the influence of lasers heating on biological tissue is analyzed. In the planning of this type of treatment, the mathematical methods are used, which allow the laser intensity to be estimated, assuring the destruction of the target region of biological tissue. To analyze the problem discussed, the mathematical models describing the laser irradiation, the temperature distribution in the tissue domain and the tissue damage degree should be formulated (Jacques and Pogue 2008;Niemz 2007;Tuchin 2007;Welch 1984;Welch and van Gemert 2011). In the case of soft tissues, the scattering generally dominates over the absorption for wavelengths between 650 and 1300 nm and to determine the diffuse fluence rate the steady-state optical diffusion equation can be taken into account (Dombrovsky 2012;Farrel et al. 1992;Fasano et al. 2010;Guo et al. 2003). The temperature field in the irradiated biological tissue can be described by the different mathematical models. Three of them result from the theory of continuum mechanics, namely Pennes (1948), Vernotte (1958) and dual-phase lag (Xu et al. 2008) models. These models are widely used in numerical modeling of thermal processes occurring in the laser-irradiated biological tissues, e.g. Afrin et al. 2012;Fanjul-Vélez et al. 2009;He et al. 2006;Hooshmand et al. 2015;Jaunich et al. 2008;Kim and Guo 2007;Narasimhan and Sadasivam 2013;Zhang et al. 2017;Zhou et al. 2009a, b. The second group of models is based on the theory of porous media (Vafai 2010) and then the biological tissue is divided into two regions: vascular region (blood vessel) and extravascular region (tissue). Using this approach, the generalized dual-phase lag model (GDPL) is formulated (Zhang 2009). This model is starting to be used to determine the tissue and blood temperatures in the heated tissues (Afrin et al. 2012;Jasiński et al. 2016;Kału_ za et al. 2017;Majchrzak et al. 2015). Knowledge of temperature distribution allows one to estimate the degree of tissue destruction using so-called Arrhenius integral (Jasiński 2013;Niemz 2007) or thermal dose parameter (Sapareto and Dewey 1984).
So far, the laser-tissue interactions have been modeled as the direct problems. For the assumed intensity of the & Marek Jasiński marek.jasinski@polsl.pl 1 Silesian University of Technology, Konarskiego 18A, 44-100 Gliwice, Poland laser, the source function related to the laser heating, the temperature distribution and finally the degree of tissue destruction have been determined. In this paper, the inverse problem connected with the estimation of laser intensity assuring the destruction of target region of biological tissue is considered. At first, the direct problem using the finite difference method is solved. To describe the temperature distribution the generalized dual-phase lag model supplemented by appropriate boundary and initial conditions is applied. The source function in the GDPL equation resulting from the laser heating is associated with the total fluence rate, where its diffuse part is determined using the steady-state optical diffusion equation. The degree of tissue destruction is estimated on the basis of Arrhenius integral. Next, taking into account these values, the inverse problem is solved using the gradient method. In the final part, the results of computations are shown.

Formulation of the Direct Problem
An axially symmetrical domain of biological tissue exposed to the laser beam, as shown in Fig. 1, is considered. Thermal processes can be described by generalized dual-phase lag model (Majchrzak et al. 2015;Zhang 2009 and T b ðr; z; tÞ ¼ Tðr; z; tÞ À eq b c b G where e is the porosity defined as the ratio of blood volume to the total volume, C ¼ eq b c b þ 1Àe ð Þqc is the effective heat capacity (q, q b , c, c b are the densities and specific heats of tissue and blood, respectively), ð Þk is the effective thermal conductivity (k, k b are the thermal conductivities of tissue and blood, respectively), G is the coupling factor, s q is the relaxation time, s T is the thermalization time, Q m and Q mb are the constant metabolic heat sources, Q ex (r, z, t) is the source function connected with the laser irradiation and T(r, z, t) and T b (r, z, t) are the tissue and blood temperatures.
The generalized dual-phase lag Eq. (1) is supplemented by boundary condition (Mochnacki and Majchrzak 2017): where C 0 corresponds to the axis of symmetry, C is the boundary of the cylinder, n is the normal outward vector, rT(r, z, t) is the temperature gradient and q b (r, z, t) is known boundary heat flux. Here, the no-flux boundary condition q b (r, z, t) = 0 is assumed. It should be noted that for strongly scattering tissues (like the soft tissues analyzed here), laser heating is considered as a body heat source but the irradiated surface is thermally insulated (Afrin et al. 2012). Thus, on the upper surface of the domain considered q b (r, z, t) = 0 is also assumed. The initial conditions are also given: where T p is known temperature and w(r, z) is the initial heating rate. The ordinary differential Eq.
(2) is supplemented by condition Source function Q ex (r, z, t) connected with the laser heating can be defined as follows (Kim et al. 1996): where l a [1/m] is the absorption coefficient, /(r, z) [W/ m 2 ] is the 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 light fluence rate /(r, z) is the sum of collimated part / c and diffuse part / d (Gardner et al. 1996;Guo et al. 2003;Kim et al. 1996) In the case of soft tissues, to determine the diffuse fluence rate the steady-state optical diffusion equation (Farrel et al. 1992;Fasano et al. 2010;Welch and van Gemert 2011) should be solved and l 0 s ¼ ð1 À gÞl s [1/m] is the effective scattering coefficient (l s is the scattering coefficient, g is the anisotropy factor).
Equation (8) is supplemented by boundary conditions The collimated fluence rate is given as (Fasano et al. 2010) where I [W/m 2 ] is the surface irradiance of laser, r D is the radius of laser beam and l 0 t [1/m] is the attenuation coefficient defined as Solution of the above-presented problems allows one to determine the Arrhenius integral (Jasiński 2013;Niemz 2007) where P[1/s] is the pre-exponential factor, E[J/mol] is the activation energy, R[J/(mol K)] is the universal gas constant and [0, t f ] is the time interval under consideration. A value of damage integral A(r, z, t f ) = 1 corresponds to a 63% probability of cell death at a specific point (r, z), while A(r, z, t f ) = 4.6 corresponds to 99% probability of cell death at this point (Chang and Nguyen 2004).

Inverse Problem
The inverse problem formulated here concerns the estimate of laser intensity I which ensures the destruction of target region of biological tissue. Thus, the following criterion can be formulated: where A m ðr i ; z i ; t f Þ is the 'measured' Arrhenius integral. Aðr i ; z i ; t f ; IÞ is the calculated Arrhenius integral obtained from the direct problem solution with current estimate of the unknown parameter I, while M is the number of points and F is the number of time steps. This inverse problem is solved using the gradient method (Kurpisz and Nowak 1995). Thus, using the necessary condition of optimum one has where (c.f. Eq. 13) while are the sensitivity coefficients and A f i ¼ Aðr i ; z i ; t f ; IÞ. Function A f i is expanded into a Taylor series about known value of I k , this means where k is the number of iteration and I k for k = 0 is the arbitrary assumed value of I, while I k for k [ 0 results from the previous iteration. Introducing formula (18) to Eq. (15) one obtains

Iran J Sci Technol Trans Mech Eng
After the mathematical manipulations one has where K is the assumed number of iterations.
To determine the sensitivity coefficients (17), the governing equations should be differentiated with respect to the unknown parameter I (Kleiber 1997;Majchrzak and Mochnacki 2014). At first, the source function Q ex (r, z, t) connected with the laser irradiation is differentiated Taking into account the dependence (11) one has W c ðr; zÞ ¼ expðÀl 0 t zÞ exp À while the function W d (r, z) results from the differentiation of Eq. (8) and boundary conditions (10). Thus r; z ð Þ 2 X : Dr 2 W d ðr; zÞ À l a W d ðr; zÞ þ l 0 s W c ðr; zÞ ¼ 0 ð25Þ and r; z ð Þ 2 C : ÀD n Á rW d ðr; zÞ ¼ W d ðr; zÞ 2 ð26Þ Next, Eqs. (1) and (2) The boundary condition (3) and the initial conditions (4) (5) are also differentiated ðr; zÞ 2 C : In Fig. 2, the summary of the identification algorithm is presented.

Method of Solution
The direct problem and additional ones connected with the sensitivity functions are solved using the finite difference method. For this purpose, the uniform spatial grid of dimensions n 9 n (h is the grid spacing) shown in Fig. 3 is introduced.

It means
Because the explicit scheme of finite difference method is considered, the stability criterion should be determined (Majchrzak and Mochnacki 2016).
Finally, the Arrhenius integral (13) and the sensitivity function (16) are calculated and At first, the direct problem is solved under the assumption that the radius of laser beam is equal to r D-= 0.001 m and laser intensity equals I = 0.3 MW/m 2 . The initial tissue and blood temperatures are equal to T p-= 37°C and initial heating rate w(r, z) = 0. In the Arrhenius integral (13): P = 3.1 9 10 98 1/s, E = 6.28 9 10 5 J/mol, R = 8.314 J/(mol K). The problem is solved using explicit scheme of the finite difference method (100 9 100 nodes, time step Dt = 0.02 s).
In Figs. 4 and 5 the distributions of collimated fluence rate / c (r, z) and diffuse fluence rate / d (r, z) are shown, Iran J Sci Technol Trans Mech Eng while Fig. 6 illustrates the distribution of source function Q ex (r, z) related to the laser irradiation. As can be seen, only in the small part of the tissue (0 B r B 0.004 m, 0 B z B 0.006 m) the source function associated with the laser irradiation is non-zero. Figures 7 and 8 illustrate the courses of the same functions, as previously in the z direction for r = 0. Maximal values of these functions are clearly visible.
The tissue temperature distributions after the time 60 and 120 s are shown in Figs. 9 and 10, while in Figs. 11 and 12 the distributions of Arrhenius integral for the same moments of time are presented. The subdomain marked in black (Figs. 11,12) corresponds to this part of the tissue which is destroyed (A(r, z, t f ) C 1).
It should be noted that to better illustrate the results of computations, in the Figs. 4, 5, 6, 11 and 12 only the fragment of the domain corresponding to 0 B r B 0.006 m, 0 B z B 0.006 m is shown.
Next, the inverse problem is solved. Thus, on the basis of the values of Arrhenius integrals, at the nodes located in the subdomain 0 B r B 0.001 m and 0 B z B 0.001 m, Iran J Sci Technol Trans Mech Eng obtained from the direct problem solution the laser intensity is identified. In Fig. 13 the results for two initial values I 0 = 0.25 MW/m 2 and I 0 = 0.5 MW/m 2 are shown. For these values, the iteration process is convergent and the final value of I is obtained after 20 and 25 iterations, respectively.
In conclusion, it is possible to estimate the intensity of the laser which will ensure the appropriate distribution of Arrhenius integral-it means the assumed tissue destruction.

Conclusions
The solution of the inverse problem consisting in estimating the laser intensity which ensures the destruction of assumed tissue subdomain is presented. The temperature field in the biological tissue is described by the generalized dual-phase lag model, while the degree of tissue destruction is determined using the Arrhenius integral. The direct problem and additional one related to the sensitivity analysis are solved using the finite difference method, while to solve the inverse problem the gradient method is applied. However, it should be noted that the gradient algorithm is not always convergent. Thus, appropriate selection of the starting point I 0 , which ensures convergence, is, therefore, extremely important.
Presented approach can be used, among others, in the case of electromagnetic field interaction on biological tissue, e.g. (Paruch 2017). For the postulated degree of tissue destruction it is possible to identify, for example, the voltage on the electrodes.