Thermal damage in three-dimensional vivo bio-tissues induced by moving heat sources in laser therapy

The thermal damage of a three-dimensional bio-tissue model irradiated by a movable laser beam was studied in this work. By employing the DPL biological heat conduction model and Henriques’ thermal damage assessment model, the distribution of burn damage of vivo human tissue during laser therapy was analytically obtained. The influences of laser moving velocity, laser spot size, phase lags of heat flux and temperature gradient were discussed. It was found that the laser moving speed and the laser spot size greatly influence the thermal damage degree by affecting the energy concentration degree. The increases of the laser moving speed and laser spot size can enlarge the irradiated region and reduce the burn degree. A greater phase lag of temperature gradient led to lower accumulation of thermal energy and lower burn degree. However, the increment of heat flux phase lag leads to the thermal energy accumulation and more serious burn degree in the irradiated region.

The thermotherapy is a useful approach in medical field, such as hyperthermia, laser soldering, laser ablation, laser surgery and other thermal treatment methods 1 . Many researchers have explored the pathological mechanism and thermal behavior in different laser surgeries. For example, Tuncer et al. 2 declared the advantages of the CO 2 laser surgery in oral soft tissue pathologies, such as minimal side effects, good pain control, minimally invasive procedure and favorable prognosis. Semenyuk 3 explored the laser-induced thermal behavior of retinal pigmented epithelium and the adjacent layers in a human eye, in order to specify laser optimal parameters and protect the normal tissue in an ophthalmic surgery.
In order to maximize the therapeutic efficiency while ensuring the patients' safety, the thermal damage induced by the laser beam need to be predicted precisely. The accurate description of the bio-heat conduction process in living tissues is primary to predict the thermal behavior. Pennes 4 was the first to develop a biological heat conduction model to describe the thermal transfer process in the vivo tissue. The Pennes' model has been employed as an efficient model in a lot of researches on bio-heat transfer. Among them, Ma et al. 5 investigated the thermal response in the laminated bio-tissue during the nanoshell assisted laser hyperthermia for subcutaneous tumor. Shih et al. 6 presented a theoretical analysis of the Pennes bioheat transfer equation with sinusoidal heat flux condition on skin surface.
However, the Fourier's law assumed that the heat propagation speed is infinite, which was not precise in transient thermal conduction. In order to describe the thermal transfer process in living bio-tissues more accurately, the non-Fourier effects should be considered. So many researchers made modifications on the Pennes' model by introducing the non-Fourier effects. Among these, two models are commonly used at present. The first one is a hyperbolic heat conduction model, which is developed by Cattnaeo 7 and Vernotte 8 . This model introduced a heat flux relaxation time to describe the thermal transfer process with a finite velocity, and is referred to as C-V model hereinafter. Zhang et al. 9 proposed a thermomechanical model for multilayered structure subjected to heat deposition with the thermal barrier under consideration. Tzou 10 improved the C-V model by considering the phase lags of both heat flux and temperature gradient and proposed a dual-phase-lag heat conduction model (DPL model). Ma et al. 11 discussed the combined effects of these phase lags and moving thermal source on a square plate. Kumar et al. 12 theoretically investigates the thermal behavior in a living biological tissue under various coordinate systems and different non-Fourier boundary conditions with the dual-phase-lag bioheat transfer

Mathematical Models
Heat conduction model. The DPL model is employed to describe the thermal transfer procedure, which is modified from the Fourier's law with two relaxation times under consideration. The DPL model can be described by the following formulation 21 , where t, T, q, → r and k represent the time, the temperature of the tissue, the heat flux, the position vector and the heat conductivity, respectively. τ q is the heat flux phase lag and T τ represents the temperature gradient phase lag. According to the DPL model, the heat conduction equation can be expressed as: where α is the thermal diffusivity.
Bio-heat transfer model of living tissue. The heat conduction process in living tissue is significantly complicated because of the involvement of the heat conduction between blood and tissues, blood perfusion in vascular beds and metabolic heat generation. With the non-Fourier effects under consideration, an effective bio-heat transfer model that includes effects of the capillary vessel systems, the metabolism and the relaxation time of temperature gradient and heat flux can be expressed as 22 : ( ) where, ρ and ρ b represent the density of the human tissue and the blood. C and C b are specific heat of the vivo tissue and blood. w b is the blood perfusion in the human tissue. T a represents the human core temperature. Q m represents heat production rate caused by metabolism and Q is energy power of external source. energy absorption model. In the above equation, the external heat source Q can be written as 23 : t where, μ t and r t ( , ) ϕ → are the absorption coefficient of material and the local energy distribution, respectively. For the three-dimensional vivo human tissues, ϕ → r t ( , ) can be expressed in the following form: where, R represents the optical reflection ratio of the human tissue. I(t) is the power density of laser source. μ t represents the optical attenuation coefficient, and φ x y ( , ) is the energy distribution on the surface.
thermal damage model. The burn damage caused by a thermal source on a living tissue is complex and multidisciplinary, which depends on the power of the heat source and its duration. Henriques and Moritz 24 proposed an expression for the denaturation process based on the first order approximation of the Arrhenius equation. The denaturation rate K can be written as: where, A is the tissue frequency factor, R the universal gas constant and E a the activation energy of the denaturation reaction. The evaluation parameter of burn damage can be expressed as: Governing equation of the problem. The human tissue is modelled as a cuboid as shown in Fig. 1. The parameters L, b and h represent the length, width and height of the cuboid tissue, respectively. The top surface of the human tissue is irradiated by a square laser spot, which moves along the middle axis with velocity v. The governing equation based on DPL model can be written as: where Q l represents the laser heating source, θ = − T T a the temperature rise, =°T 37 C a the temperature of the human body. The heat energy of the laser beam in this study is treated as a body heat source in the tissue model and the expression can be written as: represents the distribution in laser irradiation direction for the high-absorption tissues. Q 1 (x) and Q 2 (y) are the spatial distribution functions of the laser power density which can be expressed as: www.nature.com/scientificreports www.nature.com/scientificreports/ where, H(x) is the Heaviside function, and R 0 is the laser spot measure.
Considering this living tissue model as a part of human body, the temperatures on the side-faces and bottom surface are set as T a and the top surface is adiabatic. So the boundary conditions can be written as: The initial conditions can be expressed as:

Solution for Governing equation
The governing equation (8) can be solved by using the variables separation method and the mode superposition method. The solution of equation (8) can be written as following 25 : where, X m (x), Y n (y) and Z s (z) are the mode functions, which are expressed as and β m , γ n and μ s are eigenvalues, which are given by m n s Substitute equation (14) into equation (8), and T t ( ) mns can be obtained as: By introducing the following parameters where, mns m n s can be derived as the following: www.nature.com/scientificreports www.nature.com/scientificreports/  While ⁎ T t ( ) mns can be obtained as the following:   Finally, substituting the temperature distribution formulations into equation (7) yields the burn damage in the human tissue. Table 1 gives the parameters of the laser source used in this problem. Table 2 shows the thermo-physical parameters of the human tissue 26 .

physical parameters
The parameters for burn evaluation, which are introduced in equations (6) and (7), are shown in Table 3 26 . It has been accepted that: Ω = .

Results and Discussions
The empyrosis distribution in the bio-tissue due to a moving laser is derived on the basis of the analytical solution of temperature response, which is introduced in Section 3. The thermal conductivity of bio-tissue is small, so the burn mainly occurs in a small vicinity region of the irradiated area where the laser energy is concentrated. Firstly, the empyrosis distribution in an appropriate part of the target tissue (the grey part in Fig. 2(a)) is exhibited in Fig. 2(b) in order for a clarified description of the thermal damage. The numbers in the nephograms in Fig. 2(b) show the coordinates of the locations. And the color level shows the value of thermal damage.
It is obvious in Fig. 2 that the empyrosis occurs in the area subjected to the laser irradiation. The most severe burn appears along the center line of the affected area on the top surface, which agrees well with Fig. 3(b,c). The irradiation duration at the starting point of the path is shorter than the part behind, so the burn of the starting point is at a slight degree, which is also shown in Fig. 3(a,d). Figure 3 shows the influences of moving velocity of heat source on the burn. The movement of the heating source brings about energy dispersion in the irradiated area 28,29 . The power input during a specific period is constant, so that a greater moving speed makes the laser energy spread in a broader area. Figure 3(a) shows the empyrosis variation on the laser moving path. The laser beam induces serious burn in the irradiated region. And the area of irradiated region is proportional to the moving velocity, with the double speed leading to the double area. However, the degree of thermal damage decreases because of the less energy concentration caused by the higher laser moving speed. Figure 3(b,c) show the thermal damage variation along the y-and z-coordinates, respectively. The heat-transfer capability of bio-tissue is weak, so that the width of the burnt region is just a little greater than the laser spot size R 0 . The laser moving speed has a great influence on the thermal damage in the irradiated region. The laser energy concentrates in a thin layer below the top surface, so the thermal damage  Table 2. Thermo-physical parameters of the human tissue.

Temperature rage ( o C) E α /R (K)
A, Frequency factor (s −1 ) T ≤ 55 7.5 × 10 4 3.1 × 10 98 T > 55 3.54 × 10 4 5.0 × 10 45 www.nature.com/scientificreports www.nature.com/scientificreports/ maximizes on the top surface, namely z 0 = , and decreases rapidly with the increment of z. It can be found from Fig. 3(a-c) that the magnitude of empyrosis increases by hundreds of times with a thirty percent decrement of laser moving speed (from 0.006 m/s to 0.004 m/s). In order to visually demonstrate the effects of moving velocity on the degree of burn, the variations of the depth of second-degree ( 1 Ω = ) and third-degree (Ω = 10 4 ) thermal damage along the x-axis are shown in Fig. 3(d). The dashed lines show the second-degree burn while the solid lines for the third-degree burn. Figure 3(e) shows the depths of the three levels of burn vs the laser speed. It can be found from Fig. 3(d,e) that the minimum speed value causes the deepest burn wound. And the depth of thermal damage in the target tissue decreases with the increment of laser moving speed. The third-degree burn does no t occur when the moving speed increases to v = 0.008 m/s. Since the heating power density is constant, the energy absorbed in the specified region increases as the moving speed of laser beam decreases. Consequently, the burn wound in this region will be more serious and deeper.
The laser spot size will affect the energy input per unit area. With a larger spot size, a wider region is irradiated by laser beam. Since the total laser energy power is constant, a specified position of the tissue absorbs less energy in such a case. That is to say, the energy dispersion caused by a bigger spot size will lead to a lower empyrosis degree. Figure 4(a) compares the variation of thermal damage along the laser moving path, which is induced by the laser beams with different spot sizes. It can be found that the spot size has significant influence on the burn injury. When the spot size decreases from = . R 0 0015 m 0 to = . R 0 001 m 0 , the spotted area falls by 55%, while the value of Ω increases by one thousand times. Figure 4(b) shows the variation of burn injury along the y-axis, and it is shown that the maximum value of thermal damage shows up on the path of laser spot center. The decrement of the spot size leads to a significant increment of the peak value of Ω and reduces the burnt area. The concentration of laser energy caused by the decreasing of the spot size also leads to higher-level burn in the depths of the target tissue, as is shown in Fig. 4(c). Figure 4(d,e) show the influences of laser spot size on the depth of different degrees of burn injury. It is shown that the depths of the three levels of burn decrease with the increment of laser spot size and the third-degree burn disappears when = .
R 0 002 m 0 . www.nature.com/scientificreports www.nature.com/scientificreports/ It is known that the non-Fourier effect has significant influence on the empyrosis in the target bio-tissue. In the following, the influences of the two relaxation time q τ and τ T on the burn damage will be discussed. The heat flux relaxation τ q means the heat flux lags behind the temperature variation. Consequently, the thermal energy propagation from the irradiated area to the surrounding tissue is arrested. A greater τ q leads to more energy accumulated in the exposed region and a more remarkable vibration characteristics in thermal response 30 . Figure 5(a) shows the distribution of thermal damage along x-axis under different values of q τ . The fluctuation appears in the empyrosis distribution when q τ increases from 0.05 s to 0.5 s. When τ q is great enough, the fluctuant characteristics would appear in the temperature response, which has been reported in the existing literatures 22,31,32 . Since the thermal damage parameter Ω is determined by the integration of temperature function with regard to time, the fluctuation occurs as a result. The thermal damage in the affected region becomes more and y 0 009 m, z 0 = ); (b) Thermal damage variation with ycoordinate ( = .
x 0 0045 m). www.nature.com/scientificreports www.nature.com/scientificreports/ more severe with the increasing heat accumulation caused by the growing of τ q . Moreover, the increment of τ q also brings a smaller burn area, which is shown in Fig. 5(d). The variations of thermal damage along y-and z-directions are shown in Figs. 5(b,c). The thermal accumulation due to the delay of thermal energy propagation is remarkable. The increment of τ q leads to an obvious rise of thermal damage in the irradiated region. Figure 5 , the retardation of the heat flux is not distinct, and the third-degree burn does not occur in the target tissue. On the other hand, when q τ increases to 0.5 s, the hysteretic heat flux can not transfer the thermal energy on time. So more energy accumulates in the target region and the third-degree burn occurs in the depths of the tissue, as is shown by the solid lines in Fig. 5(d). The dashed lines show the depth of the second-degree burn. It is shown that the largest τ q ( 0 5 s q τ = . ) leads to the greatest depth of third-degree burn and the smallest depth of second-degree burn. However, the smallest τ q (τ = .

Scientific RepoRtS
0 05 s q ) causes the smallest depth of third-degree burn and the greatest depth of second-degree ).
burn. This phenomenon agrees well with the conclusion that the phase lag of heat flux can arrest the energy transmission and cause the thermal energy concentration in the irradiated region 30 . Figure 5(e) shows the variation of the depth of thermal damage with q τ . It can be found that the depths of first-and second-degree burn are not sensitive to the value of τ q . However, the depth of the third-degree burn increases obviously with the increment of τ q , especially in the range of 0 ~ 0.5 s.
The phase lag of temperature gradient, τ T , means that the heat flux precedes the temperature variation. The thermal transfer, which precedes the temperature increment, causes less energy accumulation and lower peak temperature in the irradiated area. The vibration characteristic in the thermal response recedes due to the increasing of T τ 30 . Figure 6(a,b) show the burn injury distribution along the laser moving path with different T τ . The increment of τ T enhances the dispersion of the thermal energy, which means that the thermal damage is milder. Moreover, the thermal response area expands as τ T increases. As a result, the phase lag of temperature gradient leads to milder burn in a wider affected region. Figure 6(c) shows the distribution of empyrosis in z-direction. The , z 0 = ); (b) Thermal damage variation with ycoordinate ( = .
x 0 0045 m www.nature.com/scientificreports www.nature.com/scientificreports/ temperature retardation means that heat transfer occurs before the temperature response, so the heat accumulation in the surface layer of the target tissue is weakened and the peak value of thermal damage decreases. When τ = .
0 05 s T , the thermal damage, Ω, decreases rapidly along the thickness direction. However, when τ T increases to 1 s, Ω decreases more gradually. The obstruction of thermal accumulation leads to lower degree of burn and wider heat-affected zone, as is shown in Figs. 6(d,e). As T τ increases, the depth of third-degree burn decreases obviously until it disappears when τ T reaches 0.5 s. However, the regions of first-and second-degree burn are not sensitive to the variation of the temperature relaxation time τ T . www.nature.com/scientificreports www.nature.com/scientificreports/ conclusions Combined action of the non-Fourier effects, the metabolic thermal production and the heat diffusion by blood circulation makes the thermal transfer procedure in human tissue very perplexing. The present study explored the thermal damage distribution in a three-dimensional model of human tissue which is irradiated by a movable laser source. DPL biological thermal transfer model and Henriques' thermal damage model were employed. The effects of the laser moving velocity, the laser spot size and the phase lags of heat flux and temperature gradient were discussed. The following conclusions can be drawn from the research: The source moving velocity and the magnitude of laser spot influence the thermal damage by affecting the energy concentration degree. The laser moving speed mainly influences the thermal damage along the moving direction, namely the x-direction, however the effects of laser spot size show up on both x-and y-directions. Increment of the source moving velocity and the magnitude of laser spot makes the laser energy distributed over a much larger area and reduces the burn intensity level. With the comparison with spot magnitude, the source moving speed has greater influence on the burn depth.
The two thermal relaxation parameters, q τ and T τ , also show distinct influences on the thermal damage. The phase lag of heat flux q τ delays the thermal energy propagating from the irradiated region to the vicinity. Consequently, the thermal accumulation and burn degree increase with the increment of q τ . The depths of firstand second-degree burn reduce gently with τ q , but the depth of third-degree burn is more sensitive. On the contrary, the phase lag of temperature gradient T τ delays the temperature response in the target tissue, implying the energy transmission occurs before the temperature varying. τ T impedes the thermal accumulation and reduces the burn degree in the irradiated region. The depth of third-degree burn decreases obviously with the increasing of T τ .

Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.