Theoretical analysis of thermal response in biological skin tissue subjected to multiple laser beams

Laser has become a well-accepted technique for surgical procedures. To fit the physically irradiated area to the clinically planned target area, a multifibre laser system can be used. This paper theoretically investigated the temperature distribution in a three-dimensional biological tissue when irradiated by multifibre lasers with the aids of the dual-phase-lag (DPL) bio-heat conduction model. First, four laser beams located on the corners of a square are used as the heat sources. The method of separation of variables is adopted to obtain the analytical expression of temperature. The characteristics of DPL bio-heat transfer model and the difference with Pennes model were checked graphically. It is found that the effect of multiple laser beams on temperature response differs from those of single laser beam distinctly. A second heating phenomenon occurs and the contour of a square with curved corners are obtained. Then the laser beams are located to form a triangle and a trapezoid and the corresponding temperature responses are studied. Different irradiated zone can be obtained and accomplished by changing the spot size, the arrangement layout and the interval distance of the laser beams.


Introduction
Nowadays cancer is one of the most common disease with high fatality rate. There are always some drawbacks in traditional treatments such as physical removal, radiotherapy, chemotherapy, and so on. Recurrence of tumor, poor effect on lesion tissue and negative effect on normal tissue are always existing defects which are hard to overcome completely. To cover these shortages, many new therapies have been developed.
The invention of laser has offered a variety of new techniques for surgical procedures. The main modalities for its medical application are photocoagulation, non-thermal ablation (optical knife) and photodynamic therapy [1]. Based on the ability to locally coagulate biological tissue, the laser-induced interstitial thermotherapy (LITT) is proposed [2]. The basic principle of LITT is to position an appropriate laser applicator inside tumor and accomplish irreversible damage to tissue by heating it above 60 • C. Compared with traditional treatments, the LITT possesses prominent advantages, such as minimal invasion, accurate location and less negative effect to normal tissue.
Bown [2] delivered high-intensity light to small well-defined areas under precise control. He pointed out that the possible precision for local treatment of solid tumors with lasers is greater than for almost any other techniques, but careful quantitative studies are necessary to establish the appropriate treatment parameters in any particular situation. After this, the LITT technique is utilized to treat 2. Mathematical model

Governing equation
In 1948, Pennes [7] proposed a conservation equation of thermal energy based on the experiment of forearm ρ t c t ∂T ∂t = − ∇ · q → + w b ρ b c b (T a − T) + Q m + Q l (1) in which T and q → are temperature and heat flux vector that are constitutive variables. T a is the temperature of arterial blood. ρ and c are mass density and specific heat capacity, respectively, with the subscript t and b denoting tissue and blood, respectively. w b represents the perfusion rate of arterial blood which means the volume of blood flowing into unit volume tissue in unit time. Q m is the metabolic heat generation from the tissue cells and Q l results from the heat generation of external laser irradiance. In combination with Fourier's law q → = − k∇T (2) in which k denotes the thermal conductivity of skin tissue, the governing equation of temperature called as Pennes model of bio-heat transfer could be obtained as [7].
in which the ∇ 2 T = ∂ 2 T ∂x 2 + ∂ 2 T ∂y 2 + ∂ 2 T ∂z 2 is the Laplace operator. For eliminating the paradox of infinite speed of thermal propagation, the modified heat conduction model is provided by Cattaneo [8] and Vernotte [9], that is where τ is the thermal relaxation time. Substituting Eq. (4) into Eq. (1), the corresponding governing equation with the aid of conservation equation of thermal energy can be deduced as follows To consider the microstructures interaction effect in small-scale and high-rate heating, the dual-phase-lag model of heat conduction was developed by Tzou [14].
which would turn to the following form with Taylor's expansion with time to the first order [14]. ( Therefore one can obtain the dual-phase-lag model of bio-heat transfer as [16]. The dual-phase-lag (DPL) model of bio-heat transfer is adopted in the following to simulate the thermal behavior with application of multiple laser beams.
The skin tissue applied to multiple laser beams are modeled as a three-dimensional cuboid, which is depicted in Fig. 1. The length, width and depth are L, L and h. Four identical laser beams irradiate on the top surface of the skin and they are placed at the corners of a square. The laser beams are simplified as square column with the side length R 0 . The distance between two neighboring laser beams is a 2 and the distance from one laser beam to the nearest boundary of the skin model is a 1 .

Disposal of laser light
Since the model is a part of body, the temperature is equal to the body temperature, before the irradiation of laser beams. That is, the initial condition could be expressed as follows In the unaffected area, the temperature keeps constant and equal to the body temperature. The top surface of the model suffered from the irradiance of laser light directly could be viewed as adiabatic. So the boundary conditions are given by The effect of laser irradiance could be regarded as volumetric heat source [33] with the expression of where μ a , ψ (z, t) and ϕ (x, y) are the absorption coefficient of material and local light distribution along depth direction and transverse direction, respectively. For highly absorbing tissue, the distribution of laser light along depth direction can be obtained according to Beer's law [34] as where I (t), R and μ t , are laser intensity, light reflectance of surface and attenuation coefficient, respectively. From Fig. 1, the distribution of heat energy induced by Q 1 could be expressed as where H (⋅) represents the Heaviside step function. Similarly, the heat energy distributions for Q 2 , Q 3 and Q 4 are given by

Analytical solution
Considering the constant initial temperature T a , the temperature increment (θ = T-T a ) results from the heat generation and is considered as direct response to laser light. Introducing the temperature increment into the temperature governing equation Eq.(8), the following governing equation of temperature rise could be deduced [25].
And the corresponding initial conditions and boundary conditions can be obtained from Eqs. (9) and (10).
The method of separation of variables and series solution are employed to obtain the analytical solution. Based on the boundary condition Eq.(16), the eigenfunctions along various direction can be expressed as Therefore the temperature increment could be assumed as Substituting Eq. (21) into Eq. (15), the governing equation for the time term Θ mns (t) can be obtained as where λ mns 2 = λ m 2 + λ n 2 + λ s 2 is the sum of squares of eigenvalues in x, y and z directions. By multiplying the eigenfunctions  (23) in which Substituting Eq. (21) into Eq. (17), the corresponding initial conditions turn into For general adaption of the solution to various arrangement of laser beams, the nonhomogeneous term on the left-hand side of Eq. (23) could be expanded specifically. If the volumetric heat source Q is expressed as , then the following expression could be obtained The metabolic heat generation could be regarded as volumetric heat source from laser irradiance with in which Therefore the governing equation can be written as the following brief form where Considering the homogeneous initial condition Eq. (25), the solution could be obtained as [35]. Defining , the concrete solution would be determined for the following cases: Case 2: Δ < 0 Case3: Δ = 0 Finally, the temperature response is obtained as

Results and discussions
In this section, the thermal response of the skin tissue under multiple laser beams is discussed by using the DPL model of bio-heat transfer comprehensively and a set of parametric studies are presented. The size of the skin model is taken as length L = 10 mm, width L = 10 mm and depth h = 9 mm. The side length of laser beam is R 0 = 2 mm. The pure interval space between laser beams is taken as a 2 = 2 mm. The thermophysical properties of skin tissue and blood and parameters of laser light are given as [22,36,37]: ρ t = 1190kg/m 3 , c t = 3600J/(kg · K), k t = 0.235W/(m · K), ρ b = 1060kg/m 3 , c b = 3770J/(kg · K), w b = 1.87 × 10 − 3 s − 1 , R = 0.024, μ a = 2000m − 1 , μ t = 2000m − 1 , I 0 = 1 × 10 4 W/m 2 , Q m = 368.1W/m 3 .In the following discussion, the temperatures at several points and a reference line, which are shown in Fig. 1 (b), will be considered. Point 1 is the center of the model, Point 2 is the middle of two neighboring laser beams, and Point 3 is the center of laser Q 3 . Line 1 is a horizontal line passing the centers of two neighboring laser beams.
The convergence of the series solution is checked and a good convergence is reached when M = 50 and S = 50, so 50× 50 items of the series will be selected in the following calculations. Due to the limitation of the paper, the detailed results are not shown here.

The characteristic of DPL bio-heat transfer model
The isotherms at different times are obtained to observe the temperature response of the irradiated area, which is shown in Fig. 2. In the calculation, the phase lags used are τ q = 16s and τ T = 8s. The laser intensity is taken as I 0 = 1 × 10 4 W/m 2 and the size of laser beam is assumed as 2 mm. The exposure time of laser irradiance is 100s. At the beginning of heating, the temperatures of the regions around the fibers increase quickly. Then the heat propagates outside and the temperature of the region among the fibers is higher than that outside of the rectangle. Since the four fibers form a rectangle, the irradiated area nearly forms a rectangle with rounded corners and curved edges. In addition, the isotherms around the laser fibers deform towards the center of the rectangle due to the interaction of the four heat sources and change into oral shapes.
It is known that the thermal propagation speed is assumed as infinite in Pennes model, while it is finite in DPL model. Under the irradiation of four laser beams located at different locations, it is clear that the heat energies from the four beams reach a specified point in different times. The points 1, 2 and 3 are three typical locations to show the propagation of heat from the four heat sources. Fig. 3 (a) depicts the temperature developments at the three points with two different phase lags in DPL model and Pennes model.
The phase lags in DPL model are taken as τ q = 16s and τ T = 0.05s. The laser light keeps on after being turned on. The difference of heat transfer behavior between DPL model and Pennes model could be observed in Fig. 3 (a). In initial stage, the temperature increases immediately and quickly in Pennes model, while the speed of temperature rise is relatively small in DPL model. However in later time, either DPL model or Pennes model can predict the same temperature increment after a relatively long time or in the stable stage, which could be obtained in theory deduction. From Fig. 3 (a), an interesting phenomenon occurs at Point 3 that a second temperature rise or second heating appears in temperature variation in DPL model due to the bio-heat transfer process. To give a good explanation to this phenomenon, the temperature responses of the skin under several cases of laser beams are shown in Fig. 3 (b), including single laser (Q 1 and Q 3 ), combination of several lasers (Q 2 +Q 4 , Q 1 +Q 2 +Q 4 ) and all the lasers. In the subfigure (b), obviously different contributions to temperature increment from Q 3 laser and other lasers can be found. Since Point 3 is located at the center of laser beam Q 3 , it will absorb the heat energy from laser beam Q 3 as soon as the lasers are switched on and the first temperature increment is obtained. After 50s approximately, the temperature rises again remarkably because of the heat transfer from the other laser beams (Q 2 +Q 4 laser beams). The laser beam Q 1 is far from Point 3 that its influence is so small that the third heating is not distinct. Obviously this time period 50s is the propagation time of thermal wave from the laser beams' locations (Q 2 and Q 4 ) to Point 3. The propagation speed can be calculated if the distance is known between the other laser source and the reference point. It is also shown in Fig. 3 (a) that there are distinct reacting times to laser heating for various locations in DPL model due to finite thermal propagation speed. Point 1 is the center of the whole model, which has the same distance away from all laser beams, so it is reasonable that the second heating phenomenon doesn't occur. By contrast the phenomenon shows up at points 3 and 2 because the heating sources have different distances from the points so that there is a sequence of response times for all the laser heating sources.

Effect of size of laser beam
To achieve more appropriate damage shape for the tumor block and better therapeutic effect, the influences of different size of laser beam and interval space between laser beams need to be studied. Firstly the interval distance between laser beams' center keeps constant as d = a 2 + R 0 = 4 mm. The size of laser beam is taken as R 0 = 1 mm, 2 mm and 3 mm, respectively. To guarantee the equal output of laser irradiance, the intensities of laser light are adopted as I 0 = 4 × 10 4 W/m 2 , 1 × 10 4 W/m 2 and 0.44 × 10 4 W/m 2 , respectively. Fig. 4 (a) depicts the temperature variation at Points 1 and 3. It can be seen that as the laser beam size increases, the steady  temperature at Point 3 decreases. At the same time, the second heating phenomenon occurs earlier due to the shortening of the pure interval distance between lasers, a 2. Different from Point 3, the temperature at Point 1 increases as the laser beam size increases, which is because the distance from the edge of laser beam to Point 1 decreases. Similarly, it takes a shorter time for temperature rising to occur at Point 1. Fig. 4 (b) shows the temperature distribution along Line 1, with the shape of a saddle. The two peaks is located at the center of the two laser beams, Q 1 and Q 4 , while the valley occurs at the middle of the two laser beams. It is acceptable that increment of laser beam size makes the temperature distribution more gentle and mild. It is shown that as the laser beam size increases, the peak temperature drops obviously while the width of the peak increases. This is because the increase of laser beam size leads to weakness of heat energy accumulation, and less energy is absorbed at the center of the laser beams but more regions are illustrated by the laser beams.
The isotherms at the surface z = 0 are depicted in Fig. 5 for various size of laser beam. The moment is t = 200s. When R 0 is small, for example, 1 mm, the isotherms around the laser beams are nearly circular, while they change into the oval shapes for the cases of larger R 0 , like 2 and 3 mm. This is because as R 0 increases, the interaction between the four laser beams becomes stronger. When R 0 is small, the heat energy is absorbed by a small region of skin and accumulate heavily, so the temperature of the irradiated region is high. However, the larger distance for heat energy to propagate can isolate all laser beams and hyperthermia zone just shows up around the laser beams, which can't live up to the treatment effect of multiple laser beams.

Effect of interval distance between laser beams d
To investigate the influence of distance between laser beams of thermal response, the size of laser beams keeps constant as R 0 = 2 mm and the distance (d = a 2 + R 0 ) is taken as 3 mm, 4 mm and 5 mm, respectively. The time histories of temperatures at points 1 and 3 are revealed in Fig. 6 (a). At Point 3, the interval distance d mainly influences the second heating stage, while the temperature starts increasing immediately after the lasers are turned on, because the heat energy of laser Q 3 is absorbed at the point. While at Point 1, the temperature starts increasing after a while. With the rise of d, the needed time increases, and at the same time the steady temperature decreases. Fig. 6 (b) describes the spatial distribution of temperature on line 1 for various interval distance at time t = 200s. Apparently the augment of interval distance would decrease the temperature between laser beams remarkably and reduce the temperature at laser beams a little bit due to the longer distance of thermal wave propagation. The temperature at Point 1 depends on the interval distance prominently, which could be concluded from Fig. 6.

Temperature distribution with triangle and trapezoid arrangement of laser beams
To fit well to the different shape of the tumor, triangle and trapezoid arrangement of laser beams are also studied, which are illustrated in Fig. 7 (a) and Fig. 8 (a). Fig. 7 depicts the isotherms at z = 0 with different internal distance d at time t = 200s for the triangle arrangement of laser beams. Like the case of rectangle arrangement, the isotherm shape is like a triangle with round corner because of the triangle distribution of laser beams. The same conclusion could be obtained that enhancement of distance would make the same isotherm more concave inwards the gap between laser beams. And the maximum temperature at the whole model lessens with increasing of the distance between laser beams. The isotherms at z = 0 in trapezoid arrangement of laser beams are plotted in Fig. 8 with different internal distances between laser beams. There is a concave dent of the irradiated region in this case. The size of laser beams is taken as R 0 = 1

Conclusions
Bio-heat transfer process of biological skin tissue applied to multiple laser beams is investigated comprehensively by applying the DPL bio-heat conduction model. Some key factors are explored to affect the thermal response, such as the size of laser beam, interval distance between laser beams and arrangement of laser beams. The following conclusions could be drawn from the above discussions: (1) Due to the finite propagation speed of heat energy, the energy of the multiple laser beams arrive at a specified point at different times, which leads to a second heating phenomenon. (2) The size of laser beam mainly affects the temperature at the laser beam and its maximum temperature. Increasing of size of laser beam would decrease the temperature at the location irradiated directly by laser beam and smooth the distribution of temperature. (3) Enhancement of the interval distance diminishes the temperature at the zone between laser beams and isolates the laser beams, which can't live up to the effect of multiple laser beams. An appropriate interval space between laser beams would maximize the therapeutic effect. (4) The arrangement of laser beams is the advantage of multifibre system compared with single laser beam system. According to the outline shape of tumor, multiple laser beams are arranged to result in a heated zone like tumor's focal zone, which can kill the pathological tissue but have little impact on normal cells.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.