Temporal and spatial temperature distributions on glabrous skin irradiated by a 1940 nm continuous-wave laser stimulator

: For predicting pain stimulation effects and avoiding damage in 1940nm laser evoked potentials (LEPs) experiments, a 2-layer finite element model (FEM-2) was constructed. A series of experiments were conducted on ex-vivo pig skin pieces to verify temperature distribution predicted by this model. Various laser powers and beam radii were employed. Experimental data of time-dependent temperature responses in different sub-skin depths and space-dependent surface temperature was recorded by thermocouple instrument. By comparing with the experimental data and model results, FEM-2 model was proved to predict temperature distributions accurately. A logarithmic relationship between laser power density and temperature increment was revealed by the results. It is concluded that power density is an effective parameter to estimate pain and damage effect. The obtained results also indicated that the proposed FEM-2 model can be extended to predict pain and damage thresholds of human skin samples and thus contribute to LEPs study. mm waist radius and 4 different power levels (0.25, 0.5, 0.75, and 1 W), as well as a laser pulse with 1 W power and 3 different waist radii (2, 3, and 3.5 mm), were used to irradiate the ex-vivo pig skin. With each combination of power and waist radius, 5 samples were tested, and each of them was irradiated for 30 s.


Introduction
The development of yttrium aluminum garnet (YAG) laser technology, (e.g. neodymium-YAG [1064 nm, Nd: YAG] and thulium-YAG [1940 nm, Tm: YAG]) has recently impacted the clinical diagnosis of pain [1]. Benefits by the effective penetration to dermis and high absorption in skin, 1940 nm Tm: YAG laser has been recognized as a novel laser source to stimulate nociceptors directly [2]. Thus, in some cases, it can effectively replace CO 2 lasers [3]. However, because of different penetration and absorption characteristics of this laser, its temperature distribution is obviously different from that of the traditional CO 2 laser. Furthermore, as 1940 nm irradiating temperature distribution for laser-evoked potentials (LEPs) was less investigated, this laser is slightly difficult to control and predict when it penetrates relatively deeply into the skin, which may cause skin damage or invalid stimulation in actual studies [1].
Pain threshold in LEPs study fully depends on the temperature rise caused by the thermal effect of laser in the cutaneous tissue [4].This begins with the absorption of the laser energy and activation of the transient receptor protein (TRP) channel family members (vanilloid receptors [VR1] and vanilloid receptor-like proteins [VRL-1]) [1,5]. The VR1 in the receptor membrane can be activated by temperatures above 43 °C, leading to a noxious stimulation of unmyelinated C fibers. Similarly, the VRL-1 can be activated by temperatures above 52 °C, and ultimately activate thin myelinated Aδ fibers [1,[5][6][7]. Damage threshold, in another aspect, relies on the temperature distribution depended on laser power, laser beam radius and exposure duration [3,8]. Considering that both tonic pain stimulation effects and skin damage depends on the time-temperature response of cutaneous tissues, a numerical model which can predict temperature distribution accurately is a feasible method to guild LEPs study by 1940 nm continuous-wave laser. It cannot only explain the stimulation effect by 1940 nm laser, but also provide guidance and prediction for the development of this laser stimulator and more effective stimulating mode.
Several authors have contributed to numerical simulations of the optical-thermal response in tissues for LEPs study. Most of them are models of CO 2 lasers. These are largely built on the simple one-layer (1-D) heat transfer models and the original works of Treede et al. [9], with many authors contributing through the addition of increased the fidelity of models [10,11]. The models that represent the skin as a two-layer (2-D) or three-layer (3-D) construct (with a Beer's law absorption term for the laser energy deposition, and with the Pennes bioheat equation used to predict the change of temperature distributions) have been shown to predict the temporal and spatial temperature distribution of skin more accurately, as well as pain stimulation effects and damage risks [12][13][14]. Moreover, the accuracy of these models has been shown to increase by opting for the 3-D model in short pulse laser studies [13]. When the parameter of intra-epidermal nerve fiber (IENF) density is applied to the model, the number of activated Aδ fiber endings can be forecast accurately [14].Although these models have been successfully applied to CO 2 lasers, there has been little success in producing models for 1940 nm continuous-wave lasers in persistent pain stimulation study. Stolarski et al. built a model of 1940 nm laser for damage threshold analysis [3]. The output of a numerical model based on a 2-D (cylindrical coordinate), finite difference solution was presented. Furthermore, the results of in-vivo experiments and tests on the skin of mini-pigs were also presented for various exposure durations and beam radii applied to skin. However, in that study, the experimental conditions were constrained due to the nature of IR surface temperature measurement and the in-vivo experimental method [3].Thus temperature in different depths of epidermis and dermis could not be measured. Moreover, the main purpose of that study was to examine laser damage thresholds, not for the pain study.
The main purpose of this study was to clarify the cutaneous temporal and spatial temperature distribution caused by 1940 nm laser irradiation. Thus, the effects of laser thermal pain stimulation and safe laser doses in LEPs study and clinical application can be predicted. To achieve this goal, a 2-layer finite element model (FEM-2) was built and a series of experiments were carried out. The FEM-2 for continuous-wave 1940 nm laser irradiation, based on Lambert-Beer's Law and Pennes bioheat equation, was constructed according to the parameters of pig skin and human skin. Temperature distributions for different spatial position of ex-vivo pig skin pieces were measured using a fast-response thermocouple in order to verify the reliability of this theoretical model. The result is beneficial to demarcate safe laser doses and exposure durations for volunteers in future pain related brain study or neuropathy patients in clinical diagnosis.

Laser stimulator and skin temperature measurement method
In this study, a continuous output Tm: YAG laser stimulator (maximum power: 5 W; wavelength: 1940 nm) was employed as the source for the persistent ex-vivo pig skin irradiation, as shown in Fig. 1. A thulium laser model (Fibolaser Ltd. Hangzhou, China) was integrated into the optical system. The output beam from the armored fiber exhibited a quasi-Gaussian distribution. The size of laser beam delivered to the target was controlled via a beam-expanding and collimating telescope (BET). Beam radii from 2 to 3.5 mm (e −2 -width) with quasi-Gaussian profiles were employed in this study. The output from the telescopes was collimated to maintain a spot size over the expected range of the target distance, which was determined from the actual permanent position of the BET frame (15 cm, verified by a vertical vernier caliper).
An armored needle wire thermocouple was used for temperature measurement. The size of the temperature measuring point was as small as 0.25 mm 2 . The thermocouple was able to penetrate deeply into the ex-vivo tissue and thereby measure the laser-induced temperature rise at different tissue depths. The signal processing circuits read the temperature signal and then transferred it to the PC in the form of RS485 serial communication. The highest sampling frequency of the temperature data was 60 Hz. A trigger module was incorporated in this measurement module to ensure that the signal sent from the laser stimulator could trigger the temperature measurement. This established synchronization between the temperature measurement and the laser irradiation, as illustrated in Fig. 1.

Preparation of ex-vivo skin and experimental tools
All experimental procedures were approved by the Institutional Animal Care and Use Committees at Tianjin. Fresh pig skin from the back of a pig which slaughtered 4 h afterward was cut into rectangular cubes measuring 4 cm long and 3 cm thick. The bottoms of these cubes were soaked in a 37 °C physiological saline solution pool to simulate an in-vivo environment. The BET was fixed on a bracket to ensure irradiating pig skin surface vertically. The solution pool was placed on a platform with a guide rail which could ensure 0.5 mm relative displacement accuracy. Accordingly, the position of this saline solution pool could be precisely adjusted, giving the pig skin relative movement to the laser spot (See Fig. 2).

Temperature measurements in experimental design
The thermocouple temperature measuring point is illustrated in Fig. 3(a). At the beginning of surface temperature measurement experiment, the thermocouple was placed in the center of laser indicator spot. It was fixed with respect to the solution pool so that the measurement point under the skin surface could not be changed. Before each irradiation, the solution pool was moved on the platform to a specified distance. Thus, the temperature distribution in a variety of positions could be obtained, starting at the central spot and then moving out toward the edge of the laser irradiation. In the internal temperature measurement experiment, the pig back skin cubes were frozen for 1 h. Then they were drilled with a hollow needle to obtain an oblique line from the surface irradiation point to 3mm subcutaneous depth. The distance from the needle's entry point on the surface to the edge of the skin cube was recorded in order to estimate the depth. The thermocouple which fixed by the platform was placed through the drilled hole. Thus, by controlling the movement distance of solution pool, we were able to adjust the depth of temperature measurement point. The experiment began after the pig skin thawed and returned to room temperature.
In order to study the correlation between laser power density and temperature distribution of the irradiated skin, a laser pulse with a 3.5 mm waist radius and 4 different power levels (0.25, 0.5, 0.75, and 1 W), as well as a laser pulse with 1 W power and 3 different waist radii (2, 3, and 3.5 mm), were used to irradiate the ex-vivo pig skin. With each combination of power and waist radius, 5 samples were tested, and each of them was irradiated for 30 s.

Finite element model for the estimation of skin subsurface temperature
A numerical model was developed to simulate laser energy deposition, heat transfer, and thereby to predict the nociceptive stimulating effect. Considering that the irradiated region of laser beam in skin was approximately a cylinder, two dimensional cylindrical coordinate was applied to this model. Since the previous research concerning CO 2 lasers showed that FEM-2 and FEM-3 predicted identical effects in long-term pain stimulation, this model was deemed sufficient for predicting the temporal and spatial distributions of continuous laser stimulation [12,13].
This model takes into account both the epidermis and dermis layers. For glabrous skin, the former is 133 μm thick and the latter is 1 mm thick. The correlation between the thermal parameters and moisture content of the skin is subject to the following equations [10,15]: where ρ is the density of skin tissue, c p is the specific heat of the tissue, k is heat transfer coefficient of the tissue, and ω c is the moisture content of the tissue. In epidermis layer, the average moisture content of the tissue is about 40%; in dermis layer, the average moisture content can reach 65% [14]. According to these equations, the thermal parameters of the model were calculated and listed in Fig. 3(b).
where P is the laser output power, R is skin reflectivity, μ α is the absorption coefficient in skin, z is the depth of the laser's cutaneous penetration, r is the distance from laser irradiation center, and w 0 is the beam waist spot size of the Gaussian distribution. S is the laser heat source. ρ b , c b , and w b represent the blood density, specific heat capacity of blood, and the blood perfusion of dermis, respectively. For expression convenience, we definedα as the product of ρ b , c b , and w b . T 0 represents the base temperature of skin before irradiation. Thus represents the influence of blood perfusion and heat dissipation to bio-heating effectiveness.
Since the output laser is continuous, the impact of blood perfusion and heat dissipation must be considered in this study [16]. When ρ b , c b , and w b denote pig skin parameters, α is equal to 0.582[J/g °Cs], and the spatio-temporal distribution of temperature output by this model can be verified by experimental data; when the above parameters are values related to human skin, α is equal to 0.272 [J/g °Cs], then the temperature field model obtained can be used to guide the process of clinical pain production [17][18][19].
The laser wavelength used in the present study was 1940 nm, which corresponding to the absorption peak of water, with the absorption coefficient at 10 4 m −1 and the penetration depth of 74 μm in water [20]. According to the pioneers' absorption coefficient data of the epidermis and dermis, the absorption coefficient of epidermis at 1940 nm is 3.31 × 10 3 m −1 , while that of the dermis is 8.83 × 10 3 m −1 [3]. The epidermal laser penetration depth was 302 μm. Therefore, the light distribution of the 1940 nm laser reached the dermis since the thickness of the epidermis in this model was 133 μm. Consequently, the heat source S should be divided into the epidermal heat source S 1 and dermal heat source S 2 . According to Eqs. (4), (5), S 1 can be written as: Considering that the bottom of epidermal was the boundary of epidermal and dermal heat source, the luminous fluxes of laser out from epidermal and into dermal were equal. Thus, S 2 can be written as: where the epidermal thickness, D 1 , is equal to 133 μm. The calculation of above all the equations was performed using ANSYS software. In particular, 27,440 hexahedrons were used in building this model. The lower boundary of the skin sample in the model was kept as a constant temperature. For the surface of the model, the convection coefficient was set as 10 and the ambient temperature was 25 °C. For the inner surface of the model, the convection coefficient was set at 50 and the ambient temperature was 37 °C. The outputs of this model were temperature-space distributions along the depth direction after 30 s of irradiation, as well as temperature-temporal distributions during the whole irradiation. These data were compared against the experimental results to verify the model's validity. Lastly, by utilizing the parameters of in-vivo human skin, this model output the temporal and spatial distributions of the human glabrous skin temperature distribution.

Data analyses and statistics
For the convenience of data analysis, the following parameters were defined: Highest ΔTemp was defined as the temperature rise from room temperature to the highest temperature during irradiation; warm feeling time was defined as the time from the beginning of irradiating to the time when the temperature reached up to 48 °C (i.e. stimulating threshold of effective first pain); effective pain radius was defined as the radius of the region for which the temperature exceeded 43 °C. The temperature curves of the actual temperature distribution were acquired by a thermocouple module with a 60 Hz sampling rate. The mean value of each combination of laser power and waist radius is presented in Fig. 4 and Fig. 5, while the standard deviation with 5 s intervals was also used in Fig. 4 in order to demonstrate the repeatability of this experiment. The temperature rise of different depths at 30 s of irradiation was defined as the highest temperature rise after 30 s laser irradiation in different subcutaneous depths, which is presented in Fig. 6. The temporal temperature data of FEM-2 were output with 4 Hz sampling rate. The equation of linear regression and the Pearson correlation coefficient (PCC) test were used to evaluate the correlation of the results of model with the actual experimental results.

Correlation of average laser power density and pain-producing area
All the experimental data in this study are presented in the form of mean value ± standard deviation. In the study of all 30 skin samples, the baseline temperature of the skin surface was 29.58 ± 1.95 °C. After the elimination of the impact of baseline temperatures on final results, the spatio-temporal distribution of the skin surface temperature field can be obtained, as shown in Fig. 4. Figure 4(a) shows the temperature variation trend on the central irradiated point of the skin surface with different powers and spot sizes. Figure 4(b) shows the spatial distribution of the skin surface temperature field under the laser beam irradiation with 1 W output power and various waist radii. The data of stimulating effects, which were defined in Sec. 2.5, are presented in Table 1. Since the distribution depth of nociceptors in human skin is 20-570 μm, the temperature rise between 500 and 600 μm was also considered in this Table [21].

Verification of FEM-2 ex-vivo pig skin model
Temperature distributions for irradiations of 0.25 W with a 3.5 mm waist radius, 1 W with a 3.5 mm waist radius, and 1 W with a 2 mm waist radius were simulated by FEM-2 in this study. The comparison between the FEM-2 model results and the actual measured values are presented in Fig. 5 and Fig. 6. A comparison of the Pearson correlation coefficient shows that in the estimation of subcutaneous temperature distribution, the correlation coefficient predicted by the model and the actual measurement value were 0.976, 0.976 and 0.981 (P<0.0001) respectively. Thus, the model shows a strong agreement with the experimental results. In addition, Fig. 6 demonstrates that the 3.5mm waist radius laser beam had a slower radial diffusion of energy than 2mm waist radius. Thus uniform temperature distribution in the target region could be acquired to avoid damage and ensure stimulation effects.   Figure 7 shows the difference of the temperature field distribution between human skin and ex-vivo pig skin under the irradiation of a 1 W laser with a 3 mm waist radius. The model shows that after 30 s of irradiation with the same laser power, the overall temperature rise in the human skin temperature field is about 4 °C higher than that in pig skin.  Figure 8 shows the temporal and spatial temperature distribution of human skin as determined by the FEM-2. The irradiating durations were 0.25, 10, 20 and 30 s with a laser power of 1 W and a waist radius of 2 mm (see Fig. 8(a)). The differences in the temperature distributions from irradiation by a 1W laser beam with 3.5 mm waist radius and 2 mm waist radius at 10 s are also presented in Fig. 8(b). The x-axis represents the horizontal distance from irradiation center, and the y-axis represents the depth of skin. The damage region and the region which can stimulate nociceptors were noted by different colors.

Discussion
The primary goal of this study was to find a reliable method to evaluate tonic pain stimulation effect and thermal damage by 1940 nm continuous waveform laser. We achieved this goal by building a mathematical model which based on experimental findings, allowing for parametric analysis and prediction of above effects relative to temporal and spatial temperature distribution. The ex-vivo pig skin experiments verified the accuracy and reliability of this model. In the reason of that pig skin has been widely used in the study of infrared laser surgery and its thermal properties have been verified, a more precise prediction of temperature distributions thereby can be utilized as a reference for selecting laser powers in the study of LEPs [17,22]. Moreover, there are some new findings in this study.
The first one is the relationship between power density and temperature rise. Transient laser stimulating studies have indicated that there is a positive correlation between average laser power density and temperature rise in the irradiated skin. More specifically, a higher average laser power density is correlated with a quicker temperature rise [9].From the parameters of Fig. 4 and Table 1, it can be ensured that this conclusion is also applicable to the continuous laser pain stimulus in the present study. From the data in Table 1 and Fig. 9, it can be seen that when under a laser irradiation of 30 s, there is a logarithmic relationship between the average laser power density and maximum temperature rise at the irradiated central point. Moreover, the temperature rise at depths of 500-600μm also has a similar relationship, as shown in Fig. 9. Consequently, a higher laser power density can cause a quicker temperature rise and a better synchronicity between continuous pain stimulation and evoked potentials. However, it can be seen in Fig. 4 that the higher power density may imply a higher possibility of overheating, which can cause irreversible damage on the skin surface. Therefore, in the practical control of continuous pain production, a laser output with a high average power density should be used to make skin surface temperature rise quickly to the pain threshold; and then, a lower power density irradiation should be applied to maintain a constant thermal stimulation, thereby generating a tonic pain. This working mode cannot only fulfill the requirements of precise control of pain-producing time, good synchronism between pain simulation and LEPs, but also prevents skin damage caused by excessive laser energy.
In previous studies of the damage threshold, beam exposure size larger than 3.5mm diameters was recognized as "safe factor" in damage threshold experiments [3,8]. In this study, the result partly supports this conclusion. In Fig. 6 and Fig. 8, it can be seen that the effective pain-stimulating radius is positively correlated to the laser waist radius. The result of 10 s laser irradiation showed that when the laser beam radius was 2 mm, damage depth (temperature above 60 °C, may cause protein denaturation) was more than 500 μm, while the pain stimulating region was limited with nearly 2mm radius. In contrary, the cloud chart of 3.5 mm radius demonstrated uniform temperature distribution beneath damage threshold, which can both elicit tonic pain and avoid irreversible damage (see Fig. 8(b)). Considering that 200 to 250 μm damage depth was regarded as a threshold to reach minimally visible lesion endpoint, larger beam irradiation was an effective method for long duration stimulation in tonic pain study [3]. Thus, the laser waist radius can be used to control the size of the painproducing area. As mentioned above, the distribution of cutaneous nociceptors in the skin is uniform. Thus the number of activated nociceptors can be controlled more precisely by controlling the laser waist radius. However, the higher laser power density in 2mm radius power irradiation caused visible damage and 3.5 mm diameter is not suitable as the safety factor anymore. But the different average power densities of these two irradiation situations can still be utilized to predict pain and damage effects (see Fig. 9). In conclusion, laser power density is a more precise parameter for pain and damage prediction. The second finding is the possibility of building a precise numerical model for 1940 nm laser in LEPs study. Some studies indicated that the energy transmission of near-infrared and mid-infrared lasers in LEPs studies is not easily estimated because lasers in this wavelength have high reflectivity and diffusivity, as well as low absorptivity in cutaneous irradiation [2]. However, the results of low power density irradiation in this study showed that this model fits the real temperature distribution quite well, and is promising in predicting the real in-vivo temperature distribution. The reason is that the 1940 nm wavelength laser, as well as the CO 2 laser, both have high absorptivity in skin. The reflectivity and diffusivity of them can also be neglected. Thus the model in this study shows that the temperature distribution is predictable. The difference between this model and the previous CO 2 laser pain stimulus model resides in the penetrating depth used. From the temperature rise of skin under laser irradiation to the depths shown in Fig. 6, and the cloud chart shown in Fig. 8, it can be concluded that compared with transient CO 2 laser stimulation, the temperature rise in skin from the laser used in this study causes smaller changes at greater depths. Consequently, the continuous pain stimulating working mode of a 1940 nm wavelength laser can stimulate nociceptive receptors directly, which is consistent with the results of past research [21].Since the 1940 nm wavelength laser can penetrate into the dermis of the skin, the setting of the heat source during irradiation is a combination of the heat source in both the epidermis and dermis. Furthermore, since the irradiation time is relatively long, the heat dissipation should also be accounted for in this model. On the other hand, a notable discrepancy between the model output temperature and actual measurement at the center of laser plot with 1W and a 2 mm waist radius could be observed in Fig. 5. This is due to the damage of skin under high laser dose irradiation. Our model assumed that the structure of skin and the moist content were constant. However, when the laser energy beyond damage threshold, skin damage and the variation of moist content would cause this discrepancy. Although the model also showed the trend of temperature variation reliably even in this situation, in further study, we will consider these factors mentioned above and get more precise model.
The third finding is the difference between pig skin and human skin. Since the temperature measuring method used in this study requires drilling cutaneous holes, the invasive method cannot be applied to volunteers and experimental animals. Consequently, the FEM-2 model is an useful tool for predicting the pain-eliciting effect and the risk of injury. Because the FEM-2 model output results of the in-vivo human skin temperature field are relatively higher than those of ex-vivo pig skin, the output laser energy density should be lower in the case of in-vivo human skin to prevent irreversible skin damage in volunteer experiments (see Fig. 7).
In summary, laser power density is a reliable parameter to evaluate pain stimulation effect and damage situation. Additionally, an effective FEM-2 module of 1940 nm laser cutaneous irradiation is useful in guiding tonic laser-evoked pain irradiation for clinical and research applications. By the results of our model, although there are some unconsidered factors such as variation of moist content and thermal damage, it can still be concluded that the laser power applied in volunteer experiments should be slightly lower to avoid damage. Our results are beneficial in interpreting and predicting the physiological responses induced by a continuous 1940 nm laser, as well as for assessing the thermal effects on tissue.