Surrogate-based optimization of thermal damage to living biological tissues by laser irradiation

The surrogate-based analysis and optimization of thermal damage in living biological tissue by laser irradiation are discussed in this paper. Latin Hypercube Sampling (LHS) and Response Surface Model (RSM) are applied to study surrogate-based optimization of thermal damage in tissue using a generalized dual phase lag model. Response value of high temperature as a function of input variables and relationship of maximum temperature and thermal damage as a function of input variables are investigated. Comparison of SBO model and simulation results for different sample sizes are examined. The results show that every input variable individually has quadratic response of maximum temperature and maximum thermal damage in highly absorbing tissues.


Introduction
Computational optimization is an important paradigm itself with a wide range of applications. In every sector in engineering, there are always needs to optimize something, whether to computational time, cost, efficiency, energy consumption, maximize the profit, output, or performance. In many cases of engineering practice, optimization of objective functions comes from measurement of physical system or from computer simulations in a straightforward way. One reason behind it is that simulation based objective function are often analytically interactable and another reason is the high computational cost of simulations. Long simulation times can be feasibly handled by using surrogate model which refers the optimization of the original objective replaced by iterative re-optimization and updating of analytical computationally cheap surrogate models.
The unique characteristics of laser makes its application dramatically increased in every sector of science and engineering. In biomedical science, for example, laser is used in many treatments such as photo dynamics therapy, cosmetic dermatology, laser mammography, and plastic surgery. Thermal damage is one of the major concern of laser application in biomedical science. Pennes bioheat equation [1] is the most widely used model to obtain temperature distribution in the living biological tissue. Welch's [2] three-step model for laser induced thermal damage in biological tissue attracted many researchers' attentions. The thermal damage of the tissue was determined based on protein denaturation evaluated by a chemical rate process equation. Fourier's law, thermal wave model [3,4] were used to model the thermal effects on living biological tissues. Tzou's [5] DPL model introduced two different time delays between the temperature gradient and the heat flux which remove the precedence assumption of the thermal wave model. It allows either the temperature gradient precede the heat flux or the heat flux to precede the temperature gradient in a transient process. The general bioheat equation can be express as The main cause of the dual phase lag phenomena in the living biological tissue is nonequilibrium between the blood and the surrounding tissue. Zhang [6] derived a generalized DPL model based on nonequilibrium heat transfer [7] in living biological tissue. It was demonstrated that the phase lag times depended on intrinsic properties of blood and tissue, blood perfusion rate, and convection heat transfer. The values of phase lag times might vary from place to place in human body. The following equation with the tissue temperature as sole unknown was derived: Afrin et al. [8] investigated temperature response and thermal damage induced by laser irradiation based on the non-equilibrium heat transfer in living biological tissues using a generalized dual phase lag (DPL) bioheat model. It was showed that the generalized DPL model predicted significantly different temperature and thermal damage compared with classical DPL and Pennes bioheat model. They also studied the effects of laser parameters such as laser exposure time, laser irradiance, and coupling factor on thermal damage in living tissues. They found that the phase lag time for heat flux had more impact on the temperature earlier, while the phase lag time for temperature gradient had more impact on the temperature later. In another research, Afrin et al. [9] studied the effects of uncertainties of laser exposure time, phase lag times, blood perfusion coefficient, scattering coefficient and diffuse reflectance of light on the thermal damage of living biological tissues by laser irradiation using a sample-based stochastic model. The variabilities of input and output materials were quantified using the coefficient of variance (COV) and interquartile range (IQR), respectively. The IQR analysis concluded that phase lag times for temperature gradient and heat flux, laser exposure time and blood perfusion rate had more significant influences on the maximum temperature and maximum thermal damage of the living tissues than the diffuse reluctance of light and scattering coefficient.
SBO techniques are concerned with accelerating the optimization of simulation problems that are expensive and time consuming. In many engineering and scientific disciplines where complex numerical simulations or physical experiments need more data by additional experiments, surrogate modeling is relatively easier and cheaper to carry out [10]. It represents a class of optimization methodologies to rapidly find the local and global optima with surrogate modeling.
To build a surrogate model, Design of Experiments (DoE) methods are used to determine the location and distribution of the sample points in the design space. The goal is to gather maximum amount of information from a limited number of sample points. There are two categories of DoE methods in the literature such as classical and modern DoE methods [11]. Modern DoE methods such as Latin Hyperbolic Sampling (LHS), Orthogonal Array Design (OAD) and Uniform Design (UD) have great advantages for deterministic computer experiments without random error as arises in laboratory experiments. Kriging model and sometimes polynomial chaos expansion are combined with LHS are often used for surrogate modeling optimization [12,13]. Uncertainty and sensitivity can also be analyzed with Monte Carlo analysis in conjunction with Latin hypercube sampling [14].
In this paper, one of the most popular DoE for uniform sampling distribution LHS is used to choose random sampling in the design scheme. LHS is a method of sampling that can be used to produce input values for estimation of expectations of functions of output variable. The surrogate-based optimization is a feasible solution where the optimization of the high-fidelity model does not work or impractical [15]. The surrogate model optimization provides an approximation of the minimizer associated to the high-fidelity model. The surrogate-based optimization process can be summarized as follows: The objective of this article is to investigate optimization of thermal damage in living biological tissue based on a generalized dual phase lag model. Surrogate-based optimization is used optimize thermal damage in highly absorbing biological tissue. In the following sections, the physical model and modeling with surrogate-based optimization are briefly summarize and results are discussed. (1 ) (1 )

Physical Model
When considering highly absorbed tissues, the laser heating is approximated as a boundary condition of the second kind and the laser volumetric heat source or laser irradiance, QL, is considered as zero. The boundary conditions are given below (1 ) where τL is the laser exposure time, in is the incident laser irradiance and Rd is the diffuse reflectance of light at the irradiated surface. The initial condition is: The damage parameter is evaluated according to the Arrenius equation [1]: where A is the frequency factor, 3.1 10 98 s -1 [2]; E is the energy of activation of denaturation reaction, 6.28 10 5 J/mol [2]; R is the universal gas constant, 8.314 J/ (mol. K); T is the absolute temperature of the tissue at the location where thermal damage is evaluated; t0 is the time at onset of laser exposure; and tf is the time of thermal damage evaluation. First performing integration of Eq. (3) over the control volume of grid point P and over the time step from t to t+Δt, then applying backward difference in time and piecewise-linear profile in space, the equation leads to the following form After replacing the values of the temperature-involved terms into Eq. (12) for the source term b, the discretization Eq. (8) becomes a linear system of algebraic equations and can be solved by TDMA (Tri-diagonal matrix algorithm). The temperature can be computed from the discretization form of the bioheat transfer model [8].

Modeling with surrogate-based optimization
There are several surrogate models available in the literatures, such as Response Surface Model (RSM), Kriging model, and Radial Basis Function (RBFs) etc. [16][17][18]. RSM consists of a group of mathematical and statistical techniques used in the development of an adequate functional relationship between a response of interest, y, and many associated inputs variables denoted by x1, x2, ….xk. A relationship between them can be approximated by a degree polynomial model of the form (13) where x= (x1, x2, ….. xk)', f(x) is a vector function of p elements, β is a vector of p unknown constant coefficients, ε is a random experimental error assumed to have a zero-mean value. The Latin Hyperbolic Sampling (LHS) is a type of stratified Monte Carlo sampling. The basic idea is to make sample point distribution close to probability density function (PDF). The coding structure of LHS is given in Fig. 2. For producing a LHS of size N [19], define P=pjk to a N×K matrix, where each column of P is an independent random permutation of (1,2…., N}. Xjk is defined by where is the NK independent and identical distributed random variables independent of P [20]. The nominal mean values of τL, wb, μs, Rd, τT, and τq are chosen as 5 s, 1.87×10 -3 m 3 / (m 3 tissue s), 120 cm -1 , 0.05, 0.05 s, and 16 s, respectively. Standard deviation is taken to be 0.5 for all the input parameters. Defined MATLAB function is used to calculate LHS for the observed domain.
A series of simulation is run using a factorial design or fractional factorial design to identify the most significant explanatory variables. A focused design is used to determine the response to the explanatory variables within a range of interest. Using least squares regression to a polynomial, RSM is produced to predict responses of values to the explanatory variables and parameters that are not simulated in the original model including an optimum (maximum or minimum) value. The inputs from LHS are x=x1, x2, x3…xk and S =[x(1), …., x(n)] T and outputs ys=[ ys (1), …., ys (n)] T . The true quadratic RSM can be written by Eq. (13). The pair (S, ys) denotes the sample data sets in the vector space. The second -degree model (d=2) can be defined as follows y= β0+ βi xi + ∑∑βij xi xj + ∑βij xi 2 +ε (15) where β0, βi, βii and βij are unknown coefficients to be determined. Next targets are first to establish a relationship between y and x1, x2, x3…xk, then determine optimum setting of x1, x2, x3…xk that result in the maximum response over a certain region of interest. The least square estimator of β is A response surface model can be generated using the least squares regression procedure. After determining all the coefficient of response surface, the approximated response can be obtained from untried point by Eq. (15). The model is specified for the problem by the following equation [21]: y=β 0 +β 1 x 1 +β 2 x 2 +β 3 x 3 +β 4 x 4 +β 5 x 5 +β 6 x 6 +β 12 x 1 x 2 +β 13 x 1 x 3 + β 14 x 1 x 4 +β 15 x 1 x 5 + β 16 x 1 x 6 +β 23 x 2 x 3 + β 24 x 2 x 4 + β 25 x 2 x 5 +β 26 x 2 x 6 + β 34 x 3 x 4 +β 35 x 3 x 5 + β 36 x 3 x 6 + β 45 x 4 x 5 +β 46 x 4 x 6 + β 56 x 5 x 6 + β 11 x 1 2 + β 22 x 2 2 + β 33 x 3 2 +β 44 x 4 2 + β 55 x 5 2 + β 66 x 6 2 (17)

Results and Discussions
Response volume of maximum temperature as function of (τL, Rd, μs) and (wb, τq, τT) are shown in Figs. 3 and 4, respectively. The figures show a general increase of temperature from the point (0,0,0) to the maximum values of the predictor variables (τL, Rd, μs , wb, τq, τT ). Afrin et al. [9] concluded that phase lag times for temperature gradient and heat flux, laser exposure time, and blood perfusion rate had more significant influences on the maximum temperature and maximum thermal damage of the living biological tissues. Therefore, the laser exposure time has chosen in this work to calculate the differences between the correct and estimated mean and variance. Table  1 shows the comparison of different of corrected and estimated mean of laser exposure time for Sample Random Sampling (SRS) and Latin Hypercube Sampling (LHS) for different sample size (N). It is shown from the table is that with the increase of sample size (N), Δμ decreases significantly for LHS than SRS which implies that LHS is better sampling procedure than SRS.  The relationship of maximum temperature to each of the response variables two at a time instead of all are showing in Fig. 6. In this case, the value of Rd at 0.05 and plotted the response of τL and μs (Fig. 5). Figure 7 (a) and (b) show the response of the thermal damage as a function of (a) τL, Rd, μs, and (b) wb, τT, and τq. The plots in figure 7 (a) show a quadratic response of thermal damage parameter as a function of μs and linear increase of damage parameter as a function of Rd and τL. Figure 7 (b) shows a quadratic response of damage parameter as a function of wb, τT and τq. The two dashed curves represent 95% simultaneous confidence band for the fitted response. The relationship of thermal damage to each of the response variables two at a time instead of all are showing in Fig 8. In this case, the value of Rd at 0.05 and plotted the response of τL and μs.

Conclusion
Surrogate based optimization has been used to optimize thermal damage in living biological tissues. Generalized DPL bioheat model based on the nonequilibrium heat transfer is applied to determine temperature response and thermal damage in highly absorbing tissues. Laser exposure time (τL), blood perfusion (wb), scattering coefficient (μs), diffuse reflectance of light (Rd), phase lag times for temperature gradient (τT), and for heat flux (τq) are chosen as input variables. It is concluded that every input variables individually have quadratic response of maximum temperature and maximum thermal damage in highly absorbing tissues. Comparison of difference between corrected and estimated mean of laser exposure time for Sample Random Sampling (SRS) and Latin Hypercube Sampling (LHS) for different sample size (N) shows that with the increase of sample size (N), Δμ decreases significantly for LHS than SRS.