Bio-heat response of skin tissue based on three-phase-lag model

In this article, the thermal response of skin tissue is investigated based on three-phase-lag (TPL) model of heat conduction. The governing equation of bio-heat conduction is established by introducing both the TPL model of heat conduction and a modified energy conservation equation. The analytical solution is obtained by adopting the method of separation of variables and a parametric study on temperature responses in TPL model is carried out. It is shown that the TPL model can predict both the diffusion and wave characteristics of bio-heat conduction. Increasing the phase-lag of thermal displacement gradient would result in the rise of thermal propagation speed and decrease the temperature in affected zone. The perfusion rate of arterial blood has no obvious effect on thermal propagation velocity and thermal propagation lagging. Increasing of the rate of blood perfusion contributes to decreasing the temperature of steady state.

The equation indicates that the heat flux and temperature gradient occur at different moments for a specified position. The delay time between the heat flux and temperature gradient is defined as the thermal relaxation time τ . The thermal wave theory ensures that the heat flux is the consequence of the temperature gradient and temperature gradient is the cause so that there is a strong path dependency for the temperature gradient 4,5 .
Later, Tzou 8 proposed a dual phase lag (DPL) model by introducing the phase lag of heat flux which is the time delay due to the fast transient effect of thermal inertia. The conduction relation between heat flux vector and temperature gradient can be expressed as where τ q is the phase lag time for the heat flux vector, and τ T is the phase lag time for the temperature gradient.
(1) � q(� r, t) = −k∇T(� r, t) The above three heat conduction models are widely used by researchers and the effects of various thermophysical properties on thermal responses have been discussed adequately. For example, Lin and Li 9 obtained the analytical solutions of bio-heat transfer for skin tissue with general boundary conditions in the Pennes, C-V and DPL models by using separation variable method. They discussed the effect of phase-lags, boundary conditions and heat conduction models on temperature and thermal damage. Xu et al. 10 presented a comprehensive literature review pertinent to the bio-thermomechanical behavior of skin tissue which reviewed the fundamental concepts, methodologies and various thermophysical properties. Based on the Laplace transform method, Liu 11 theoretically investigated the thermal behavior in a living tissue subjected to constant, sinusoidal, or step surface heating with the thermal wave model of bio-heat transfer. Babaei and Chen 12 utilized the hyperbolic heat conduction model in a functionally radius graded hollow cylinder and analyzed the effects of thermal relaxation time and nonhomogeneous indices by numerical inversion of the Laplace transform. Akbarzadeh and Chen 13 studied heat conduction in one-dimensional functionally graded media based on the DPL model.
Recently, Choudhuri 14 established Three-phase-lag (TPL) constitutive model (4) by introducing the phaselag of heat flux vector, phase-lag of temperature gradient and phase-lag of thermal displacement gradient in heat conduction law based on a model of thermoelasticity developed by Green and Naghdi 15,16 . If the thermal conductivity and phase-lags are assumed as zero, the reduced heat conduction model would be obtained like equation (5) which could be regarded as a special C-V heat conduction model when the thermal relaxation time τ approaches to infinity. Hence the three-phase-lag model based on G-N model could probably predict the non-Fourier phenomenon in medium with large phase lag time, like sand and biological skin. It is necessary to investigate the bio-heat transfer behavior of biological skin tissue induced by external heat source in TPL model.
Besides the temperature gradient which serves as a constitutive variable in C-V and DPL heat conduction models, thermal displacement gradient is considered as a third constitutive variable in TPL model. The introduction of TPL model provides a general theoretical heat conduction model and a new approach to predict thermal response and damage of structures 17 . Akbarzadeh et al. 17 studied heat conduction in a functionally graded, infinitely-long hollow cylinder based on TPL model. Quintanilla and Racke 18 obtained conditions on the phaselag and material parameters to guarantee the exponential stability for two cases in the heat conduction theory with three-phase-lag. And generalized thermoelasticity based on three-phase-lag effect is investigated [19][20][21] . Kar et al. 19 studied the thermoelastic responses in a functionally graded orthotropic hollow sphere due to a sudden temperature change in context of three-phase-lag model of generalized thermoelasticity. Kothari et al. 20 attempted to find the fundamental solutions in the context of generalized thermo-elasticity with three phase-lags. Mukhopadhyay et al. 21 presented a thorough analysis of effects of phase lags on wave propagation and investigated the nature of distributions of different fields in a thick plate based on thermoelasticity with three phase-lags and thermoelasticity of type GN-III (without any phase lag). Biswas and coauthors [22][23][24][25] published a set of researches based on TPL model, such as the propagation of Rayleigh waves in orthotropic thermos-elastic half-space, the solution in closed form for an electro-magneto thermoelastic coupled two-dimensional problem in orthotropic half space, the various heat source responses in a transversely isotropic hollow cylinder, and free vibration of homogenous isotropic cylinder panel with voids. Li et al. 26 investigated the transient thermoelastic responses of bi-layered skin tissue with temperature-dependent thermal material parameters in the context of generalized thermoelasticity without energy dissipation.
To the authors' knowledge, the three-phase-lag effect hasn't been considered in bio-heat transfer model in biological skin tissue. With this motivation in mind the aim of present analysis is to investigate the bio-heat transfer model with three-phase-lag effect in biological tissue. The governing equation is established and an analytical solution is obtained by using the method of separation of variables. A complete and comprehensive analysis on the relative parameters and comparison are presented. Choudhuri 14 introduced a constitutive heat flux relation with three phase lags in which q , T and v are heat flux vector, absolute temperature and thermal displacement which satisfies v = T , respectively. τ q , τ T and τ v are the phase-lags of heat flux vector, temperature gradient and thermal displacement gradient, respectively. t is time and r is the position coordinate vector. k and k * are the thermal conductivity and the rate of thermal conductivity, respectively. ∇ denotes the gradient operator. Equation (6) can be rewritten into the following form like DPL model.

theoretical formulations
in which the term ( T + k * v/k ) can be viewed as an new integral variable and therefore Eq. (7) could be regarded as a kind of extension of DPL heat conduction equation (3). Now Taylor's series expansion of Eq. (7) up to the first-order terms in relative with τ q , τ T and τ v leads to the following generalized heat conduction equation at a point r and time t. www.nature.com/scientificreports/ By using the relation v = T , equation (8) can be rewritten into the following form which is comprised of two variables, namely the temperature T and heat flux q.
For k * = 0, Eq. (9) would reduce to the following equation which can be simplified into the DPL heat conduction equation by integrating with time.
Furtherly, it could reduce to the thermal wave model of Eq. (2) when the phase-lag of temperature gradient is set to zero, namely τ T = 0 . And in the case of k * = 0 and τ q = τ T = 0 , equation (9) would be reduced to the heat conduction equation of Fourier's law. According to Pennes model 3 , the energy conservation equation of classical bio-heat transfer could be written as where the subscript 't' denotes skin tissue and 'b' represents arterial blood flowing into tissue. ρ and c are density and specific heat, respectively. ω b is the perfusion rate of blood which represents the volume of blood flowing into unit volume skin tissue per unit time. T a is the reference temperature of arterial blood entering skin tissue, which keeps the same as constant temperature of inner body, hence the second term on the right-hand side is regarded as heat exchange between arterial blood and the local tissue. Q m and Q laser are heat generation contributions to tissue which are provided by metabolic process and external laser heating, respectively.
It is worthwhile to mention that Eq. (12) cannot be used directly for the TPL model of bio-heat transfer, and a modification is needed to introduce. The necessity and reason of the modification and substitution to energy conservation equation is illustrated in Appendix. According to Eq. (7), the term (T + k * v/k) is viewed as an integral variable, so it takes place of the temperature T in Eq. (12) as well. Then the following energy equation would be adopted in the present study, It is easy to find that Eq. (13) could reduce to the classical energy conservation Eq. (12) in the case of k * = 0. The thermal displacement v can be eliminated from Eq. (13) by using the relation v = T and a partial differential equation of T and q can be obtained as By eliminating heat flux q from Eqs. (9) and (14), the bio-heat transfer equation for temperature based on TPL model in homogeneous skin tissue can be obtained as When k * = 0, equation (15) can also reduce to DPL bio-heat conduction by integral over time, which is By introducing the temperature increment θ = T -T a and considering T a keeps constant, equation (15) can be changed into the following equation solved in the paper.

Solution to the problem
In this section, the method of separation of variables is taken to solve the problem and obtain analytical solution for the three-phase-lag bio-heat transfer model. Two kinds of boundary conditions are under consideration to stimulate two load cases in reality.
case i: constant temperature on the surface of skin. This case stimulates the bio-heat transfer of skin tissue when the external surface of skin is abruptly attached to a high temperature object and the temperature of outer surface is viewed as the temperature of the object's temperature. A step heating on the outer surface is adopted, which means the outer surface is heated to reach the specific temperature of T 0 and then it keeps constant temperature for simplicity. The inner surface contacts with the arterial blood of the body, so it is assumed to keep the same constant temperature as the body temperature T a . The thickness of the skin is L. The boundary conditions are defined as where H(t) is the Heaviside function and θ 0 = T 0 -T a .
Case II: Constant heat flux exerted on the surface of skin. This case predicts the bio-heat transfer behavior when the highly absorbing skin suffers from a laser light irradiance. In such a case, the absorbed laser energy could be viewed as boundary heat flux 27,28 as where R d is the diffuse reflectance of light at the irradiated surface. φ in is the intensity of laser irradiance varying with time as follows.
where Q 0 is the maximum intensity of laser. The internal surface is still assumed as constant temperature T a like the second term of Eq. (18). Therefore the boundary condition for the case of constant heat flux could be written as Solution of the problem. In practice, a one-dimensional model of the tissue is sufficient as the size of laser beam is much larger than the thickness of tissue. The x-axis points inwards the skin tissue with the original point positioned on the outer surface of tissue. The temperature increment governing equation (17) turns to the following form considered in this paper.
The initial conditions are set as To deal with the nonhomogeneous boundary condition and find out the steady state of eigen mode of temperature increment, the following equation is deduced by setting all the time derivative terms in Eq. (22) to be zero, The general solution of above equation is θ = c 1 e x + c 2 e − x . Firstly with the boundary condition (18) of case I, the solution of Eq. (24) can be expressed as The steady state solution that satisfies boundary conditions equation (18) can be expressed as Therefore temperature increment can be expressed as the following form.
By substituting Eq. (27) into Eqs. (22), (18) and (23), the following equation for the transient temperature  where V 0 (t), V 1 (t), V 2 (t) are effect functions of initial conditions on the temperature response, which has four different formulations decided by the coefficients of Eq. (36). By defining the four cases are shown in the following: α, β are the real part and imaginary part of the complex number 3 − q 2 + i √ − respectively. Since the solution of θ 1 (x, t) is obtained, the veritable temperature increment θ(x, t) is expressed as As for the case II of constant heat flux, only some modifications needs to perform to obtain the analytical solution by the method of separation of variables. With boundary condition equation (21), the solution of equation (24) can expressed as Then the transient problem θ 1 (x, t) cooperate with the following homogeneous boundary condition and nonhomogeneous initial condition equation (30).

Results and discussions
In this section, the analytical results would be presented and illustrated. For case I, the outer skin surface is specified to be heated to a temperature of T 0 = 80 °C, and the temperature of arterial blood keeps constant as T a = 37 °C. For case II, the maximum intensity of laser irradiance is considered as Q 0 = 2 W/cm 2 and the diffuse reflectance is taken as R d = 0.05 28 . The thickness of the tissue is deemed as L = 9 mm. The thermophysical parameters of the tissue and blood are 18,[29][30][31] : ρ t = 1190kg/m 3 , c t = 3600J/ kg · K , ρ b = 1060kg/m 3 ,c b = 3770J/ kg · K ,ω b = 1.87 × 10 −3 s −1 ,k = 0.235W/(m · K),k * = 0.1W/(m · K · s) . The three phase lags τ q , τ T and τ v may be assumed as different values in the following discussions, and their values will be given together with the discussions. The rate of thermal conductivity hasn't been studied in bio-transfer of biological tissue before, so the constant value of the constant rate of conductivity is assumed based on request of solution 18 and relative value to thermal conductivity 23,25 . convergence of series solutions. Figure 1 depicts the time-history of temperature with different number of terms in series solutions in TPL model. The phase lags are τ q = 16 s, τ T = 6 s 32 and τ v = 2 s 18 , which satisfy the inequality 0 ≤ τ v < τ T < τ q 14 . It is necessary to mention that τ v is assumed as 2 s based on the stable request of solution 18 and relative size of numerical value between the phase lags of heat flux, temperature gradient and thermal displacement gradient 14,22,24 . It is shown in the figure that the results converge quickly as the number increases. The curves for the cases of n = 300 and n = 400 are almost coincided together, so 300 terms in series solution could be adopted in the present calculation and analysis.

Validation of the solving method.
To check the validation of the solution procedure and method used in this paper, the same problem in Liu's work 11 is calculated with the same parameters being used. According to reference 11 , the following parameters are applied: ρ t = 1000 kg/m 3 , c t = 4200 J/(kg·K), ρ b = 1000 kg/m 3 , c b = 4200 J/ (kg·K), k = 0.2 W/(m·K), ω b = 0.5 × 10 -3 s -1 , and L = 0.0012 m. These parameters are just used in this subsection. The comparison between the present results and Liu's results is depicted in Fig. 2.
In Liu's work 11 the Pennes model and C-V model are adopted. In this present model, with assumption of τ v = 0 s and τ T = 0 s, the model can reduce to C-V bio-heat transfer model. Furtherly, it reduces to Pennes model when τ q = 0 s. The time history of temperature at the location of x = 0.002 m for the Pennes model and C-V model of present model and Liu's model are shown in Fig. 2. Temperature in the skin increases immediately when the outer surface is heated. The C-V model is adopted with phase lag of heat flux taken as τ q = 20 s. It is shown that the temperature doesn't change at once, while it jumps to a high value at the time of t = 41 s. The thermal wave speed in C-V model is given by V C−V = k τ q ρ t c t , which gives the value of 4.88 × 10 −5 m/s under present configuration. At the time of t = 41 s, the thermal wave arrives at the location of x = V·t = 0.002 m, which is just the location under study in Fig. 2. Good agreement is obtained between this work and Liu's results in Fig. 2. Hence the method used in present study is efficient and will be adopted in the following calculations.
Parametric study of influences on temperature response in TPL model for the case of constant surface temperature. Comparison between TPL model and former models. Figure 3 shows the compari-  Figure 3a shows the time-history of temperature at the location x = 0.5L for the bio-heat models under discussion. It is found that the thermal propagation speed is the largest in Pennes model and C-V model predicts the lowest one. It is apparent that the curve in TPL model is between the curves of C-V model and DPL model, which means that the TPL model can show both the diffusion and wave characteristics of bio-heat transfer. The thermal propagation speed in TPL model is smaller than that in DPL model so that the thermal response in TPL model occurs later than that in DPL model. Figure 3b depicts the temperature distribution for the above models at the time t = 50 s. It can be seen that there is a drop of temperature in C-V model, which shows the location where the thermal wave arrives. Before this location, the temperature in C-V model is larger than those in the other three models, while the temperature predicted in Pennes model is the largest after this location. These could be reasonable because in C-V model the lowest thermal propagation speed would lead to more heat concentrated in the affected area by surface heating and on the contrary, the more thermal energy would transfer to the unaffected zone quickly due to the infinite wave propagation speed in Pennes model.
Effect of τ v on temperature response. Due to the uninvestigated problem in biological tissue and lack of specific value for the phase lag of thermal displacement gradient, the influence of the phase-lag of thermal displacement gradient τ v needs to be investigated and the results are shown in Fig. 4. In present paper, all the values of phase lag times attempt to be determined by the stable request of temperature solution 18 , the fundamental inequality .   Figure 4b depicts the temperature distribution at the time t = 50 s for different phase-lags of thermal displacement gradient. It can be found that the temperature response behavior would get closer to the C-V model of wave characteristic when the phase-lag of thermal displacement gradient τ v is getting smaller. It is apparent that the decreasing of τ v requires more time to motivate thermal response at the same location from Fig. 4a and the larger τ v leads to the longer thermal propagation distance at the same time from Fig. 4b. That indicates the increasing of τ v would increase the thermal propagation speed and accelerate thermal propagation process so that there would be more thermal energy transfer to the unaffected area from the affected zone. Therefore it can be seen from Fig. 4b that the temperature with large τ v is lower than that with small τ v at the location near to the surface while temperature is lower with small τ v away from the surface.
It is worthwhile noticing that small oscillation occurs in temperature response when τ v is zero, which hasn't been shown in former bio-heat transfer model before. This phenomenon is more obvious at the location nearer to the heating surface in Fig. 4c,d which results from the thermal propagation wave by instantaneous heating. The wave characteristic in space is similar to the vibration phenomenon of a string motivated by constant transverse displacement at left end abruptly in Fig. 4d. Therefore the wave characteristic would be more obvious for heating way of instant temperature rise when phase lag of thermal displacement gradient is taken as zero (τ v = 0 s). This indicates the phase lag of thermal displacement gradient will be small if the wave characteristic of bio-heat transfer in some biological tissue and the phase lag τ v is significant in bring about the oscillation characteristic.
Effects of τ q and τ T on temperature response. The effect of τ q on temperature response is illustrated in Fig. 5a,b.
As mentioned in C-V model, τ q indicates a lag phenomenon of heat flux which means thermal energy would not transfer instantly. The larger τ q is, the more obvious the lagging characteristic of thermal behavior, which results in higher temperature in affected area (Fig. 5b) and steeper rise of temperature (Fig. 5a). Comparing with τ q , the influence of τ T is smaller on lag characteristic and thermal propagation speed, which is shown in Fig. 5c,d.   Figure 6a and has no obvious effect on Non-Fourier effect at time t = 50 s from Fig. 6b. It could be reasonable and acceptable that heat transfer quantity increases with the increasing of perfusion rate of blood because the flowing blood would take more thermal energy away from subjected skin tissue which leads to the lower steady-state temperature in local tissue.
Parametric study of influences on thermal behavior in TPL model for case of constant surface heat flux. Due to limitation of the paper length, only the influences of τ v and ω b on temperature response are discussed for this case. Figure 7a depicts the variation of temperature at x = L/2 with different phase lag times of thermal displacement gradient. Reduction of τ v makes the variation of temperature closer to C-V model but in a smoother way, which indicates the increasing of τ v could enhance the speed of thermal propagation a little. This could be concluded in Fig. 7b which plots the temperature distribution at time t = 50 s. However the discrepancy with Fig. 4 is the oscillation phenomenon doesn't arise under the present configuration, which is attributed to the difference of boundary conditions. A specific temperature is applied on the surface of skin tissue in Fig. 4, while a specific heat flux is loaded on the surface in Fig. 7. One case is field variable increment, while the other is gradient rise of the corresponding field variable which leads to gentle increment of the field variable rather than abrupt increment or oscillation increment.
The effect of ω b on temperature response is shown in Fig. 8. Blood perfusion in energy conservation equation means the heat transfer to blood from affected skin tissue. Hence enhancement of ω b represents more blood flowing and participating heat convection, which leads to lower temperature in the same location at the same time, as is shown in Fig. 8a. In Fig. 8b no difference between different blood perfusion rate at time t = 50 s is attributed to fact that the blood perfusion primarily affects the stable solution of temperature according to equation (25) and hardly shows discrepancy in short time transient process.

conclusions
The thermal response behavior in a one-dimensional skin tissue model which is subjected to surface heating and heat source of laser irradiance is investigated based on the TPL bio-heat conduction model. By combining the TPL heat conduction equation and the modified energy conservation equation, a new temperature governing equation has been obtained in the heat transfer process of skin tissue. The influences of the phase-lags of heat flux, temperature gradient and thermal displacement gradient and the rate of perfusion of arterial blood on temperature responses are investigated emphatically. The comparison between TPL, DPL and C-V and Pennes models has been made to find more attractive phenomenon about TPL bio-heat conduction theory. The following conclusions can be drawn from the results and analysis: (a) The new bio-heat conduction model with three phase-lags predicts both diffusion and wave characteristics of bio-heat transfer. Therefore the behavior of temperature obtained from TPL model is between that of the C-V model and DPL model.  www.nature.com/scientificreports/ (c) The phase lag of thermal displacement gradient τ v is helpful to decrease local tissue temperature. And temperature increment oscillates for a while with constant temperature of surface boundary condition when the phase lag of thermal displacement gradient is assumed as zero (τ v = 0). (d) The increase of the rate of perfusion of arterial blood ω b leads to more heat taken away from local affected tissues resulting in lower temperature of steady state.

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