Open Access
Add to BibSonomy
Abstract
Advancements in medical laser technology have paved the way for its widespread acceptance in a variety of treatments and procedures. Selectively targeting particular tissue structures with minimally invasive procedures limits the damage to surrounding tissue and allows for reduced post-procedural downtime. In many treatments that are hyperthermia-based, the efficiency depends on the achieved temperature within the targeted tissues. Current approaches for monitoring subdermal temperature distributions are either invasive, complex, or offer inadequate spatial resolution. Numerical studies are often therapy-tailored and source tissue parameters from the literature, lacking versatility and a tissue-specific approach. Here, we show a protocol that estimates the temperature distribution within the tissue based on a thermographic recording of its surface temperature evolution. It couples a time-dependent matching algorithm and thermal-diffusion-based model, while recognizing tissue-specific characteristics yielded by a fast calibration process. The protocol was employed during hyperthermic laser treatment performed ex-vivo on a heterogeneous porcine tissue, and in-vivo on a human subject. In both cases the calibrated thermal parameters correlate with the range of values reported by other studies. The matching algorithm sufficiently reproduced the temperature dynamics of heterogeneous tissue. The estimated temperature distributions within ex-vivo tissue were validated by simultaneous reference measurements, and the ones estimated in-vivo reveal a distribution trend that correlates well with similar studies. The presented method is versatile, supported by the protocol for tissue-specific tailoring, and can readily be implemented for temperature monitoring of various hyperthermia-based procedures by means of recording the surface temperature evolution with a miniature thermal camera implemented within a handheld laser scanner or similar.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Pulsed laser light has become widespread in medical applications involving ablative or non-ablative dermal treatments and hyperthermia-based procedures, to cite a few [1]. The approach of selective photothermolysis is often employed to target the preferential chromophore without causing a significant damage to the surrounding tissue [2]. In fact, the effectiveness of such therapies strongly relies on the ability to closely monitor the temperature of the targeted tissues [3,4]. These are usually located in the subsurface regions, thus high-energy visible and infrared lasers are employed to heat targets well below the skin surface [5]. Light absorption by melanin causes epidermal damage, which limits maximum fluence that can be used in these procedures. Consequently, laser irradiation is often used in conjunction with various cooling systems to prevent the surface skin layer from overheating [5–8]. Combination of spatial selectivity of surface cooling and deeper penetration of laser heating (volumetric effect) lead to a colder superficial layer and deeper hot region within the subsurface layer. In fact, higher cooling capacity increases heat diffusion towards the surface, resulting in peak temperature deeper below the surface [7,8]. Temperature peak can be significantly higher than skin temperature (10–20 °C), while its location was reported up to 10–15 mm below the surface during Hyperthermic Laser Treatment (HTLT) [7–10]. In fact, its monitoring becomes challenging with a non-contact techniques that feature small penetration depths (e.g., thermal camera).
Current monitoring methods are either complex and invasive or offer inadequate penetration or spatial resolution. Among them, the most used are: thermocouples or fiber-optic sensors [11], optoacoustic monitoring [12,13], magnetic resonance imaging [14], different optical methods [15], and ultrasound [16,17]. Complex analytical or numerical reconstruction approaches [18–20] are also employed for the prediction of the temperature distribution within the tissue, however, they are computer intensive and can be inaccurate due to the instability of the inverse problem of the heat-diffusion equation. Therefore, alternative approaches relying on numerical simulations of laser-tissue interactions and optimization of laser parameters by ex-vivo studies have been suggested [8,21,22]. These numerical simulations are designed for specific cases [23], in which the model of laser-tissue interaction primarily depends on the laser parameters (wavelength, irradiation time, power, beam width and shape), optical tissue parameters (absorption, scattering, anisotropy) and thermal tissue parameters (thermal conductivity, density, specific heat capacity and blood perfusion) of the irradiated tissue [24].
In fact, the chromophore’s location and the absorbed laser wavelength must be considered when selecting the correct laser source for clinical use. For instance, laser sources such as Er:YAG (2,940 nm), Er:YSGG (2,780) and CO2 (10,600 nm) are primarily absorbed in water, whereas the alexandrite laser (755 nm) primarily targets melanin or hemoglobin [2], thus making numerical models highly therapy-tailored and not suitable for widespread applicability. In addition, such simulations require investigation of the tissue’s optical properties [25], which vary between patients, within the patient and over time [26,27], as a function of temperature [7,28], and skin blood volume fraction [7].
In our previous study, a time-dependent algorithm was developed and conceptually showcased on a homogeneous layer of porcine fatty tissue [9], where the sample’s thermal parameters were sourced from the literature. Here, we present its extended applicability to two-layer tissues (skin plus fat) with higher structural and thermal complexities, including temperature-dependent blood perfusion for a later in-vivo study. The algorithm was employed during a Hyperthermic Laser Treatment (HTLT) performed ex-vivo on a porcine tissue and in-vivo on a human subject, with 755 and 1,064 nm laser wavelength excitations. Examples of clinically performed HTLT are hyperthermic laser lipolysis performed with the 1,064 nm laser [7,8], or “in-motion” laser hair reduction performed with the 755 nm laser [29,30]. Moreover, the algorithm is coupled with a newly introduced tissue parameter estimation (TPE) protocol for fast calibration of the tissue’s thermal parameters, based on Er:YAG laser irradiation. This process allows for the creation of a personalized database prior to the HTLT treatment, and modifies computations for a spatial temperature distributions (STDs) estimation algorithm.
We show the estimation of thermal parameters as well as the estimation of STD within the tissue (ex-vivo and in-vivo) during the HTLT treatments. In both case studies, the thermal parameters correlate with the range of values reported in other research works. For the ex-vivo experiment, the protocol was validated by comparing the estimated STDs to the ones directly measured with a thermal camera. For the in-vivo study we found that the combination of the TPE protocol and the time-dependent algorithm reveals that STDs are in good agreement with similar studies of hyperthermia procedures. The protocol is optimized on the computational level and adaptable for different types of patients.
2. Methods
2.1 Outlining the algorithm
Firstly, we introduce the protocol for thermal tissue parameter estimation (TPE). The process starts with a calibration measurement, and includes both the active period (laser-on), where the surface is irradiated with an Er:YAG laser source, and the relaxation period, where the tissue is left to cool down via natural convection and thermal diffusion (laser-off). This particular laser source is used due to its wavelength (2,940 nm) being primarily absorbed in water as a target chromophore, and making the protocol skin type independent. Moreover, it features a low optical penetration depth, and can be treated as a surface-originating heat source. The average intensity of irradiation required to maintain temperatures below the tissue’s structural change threshold is low. This physical simplicity offers non-complex mathematical implementation and thereby short processing times, enabling the numerical model to replicate the surface temperature (Ts) response to laser irradiation (see sections (2.3) and (4.) for detailed explanation). The tissue-tailored parameters are yielded by fitting the simulated Ts to the measurements with a thermal camera (see Fig. 1).
Fig. 1. The tissue-tailored approach is enabled by the TPE protocol, where tissue thermal data is acquired by fitting a simulated tissue response to the one measured with a thermal camera. The estimation algorithm links the TPE outcome and yields a personalized, time-dependent database by calculating the evolutions of arbitrarily generated STDs. Finally, the STD within the tissue is estimated by correlating the Ts evolution during the relaxation period with the ones clustered in a database.
Secondly, we utilized a time-dependent data matching algorithm to estimate spatial temperature distribution (STD) within the tissue. Uniquely, the algorithm relies on the thermal recording of Ts during the relaxation period of the procedure (laser-off). Therefore, it is not affected by the specific laser source or variation in the tissue’s optical parameters, instead it is primarily dependent on the thermal parameters. The algorithm begins with an arbitrary function to generate STDs and encompasses a bioheat-diffusion-equation (BDE) to calculate field evolutions and yields the time-dependent database. The principle of an algorithm was conceptually presented in our previous study [9]. Here we extend its applicability to tissues with higher complexities. In addition, by introducing the TPE and linking its outcome to the BDE calculations, the algorithm offers a viable approach to personalized hyperthermia-based procedures. The temperature within the tissue is estimated during the interruption of an ongoing procedure, as the data matching process correlates the measured Ts evolution to the database and offers the best fitting dynamic, and, consequently, the best fitting STD.
2.2 Experimental setup and measurement protocol
The TPE experimental setup included a 2,940 nm Er:YAG laser source (Dynamis SP, Fotona, Slovenia), and a thermal camera (ThermaCAM P60, FLIR, 7.5–13 µm). The laser beam illuminated the surface through a diverging lens, placed 15 cm above the sample surface, effectively producing a laser spot approx. 6 cm in diameter. Due to the large laser spot, low average intensities were achieved, maintaining the surface temperatures below the tissue structural-change and damage threshold [31,32]. The active cycle was performed for 30 seconds using VLP mode (“Very Long Pulse”, pulse duration tp=1,000 µs, repetition rate fr=12 Hz). We considered three calibration measurements, where an average intensity of irradiation was systematically increased with each following measurement (from 0.08 to 0.12 W/cm2). Corresponding maximal fluence was Fmax=100 J/m2, which is below the maximum permissible exposure (MPE=998 J/m2) according to the standard of laser safety products (IEC 60825-1). Individual measurements were separated by 30 minutes, allowing for thermal relaxation of the tissue. The thermal camera recordings of the active and relaxation periods were used to estimate the tissue parameters by fitting the TPE simulation of Ts evolution to the measured one. The protocol was used to obtain tissue parameters of an ex-vivo sample and in-vivo human subject. The ex-vivo TPE was performed on a fresh porcine tissue sample (100 × 100 mm2 and 38 mm thick – 2 and 36 mm of skin and fat respectively), and in-vivo TPE was performed on a lumbar and iliac abdomen regions of 28-year-old healthy individual.
The second part of the study involved the STD algorithm that was used to estimate the temperature distribution within the tissue during the HTLT procedure. The HTLT is divided into two periods – the active period (laser radiation on in conjunction with chilled air cooling) and the relaxation period (laser radiation off with natural convection cooling). During the active period, the surface is simultaneously irradiated and cooled, effectively elevating the temperature within the tissue, while maintaining the surface temperature below the damage threshold.
An experimental setup was designed to resemble the realistic procedure, and included two laser systems, a 755 nm Alexandrite laser source, and a 1,064 nm Nd:YAG laser source (both from AvalancheLase LXP, Fotona, Slovenia). For both laser sources, the scanner head (LX-Runner, Fotona, Slovenia) was used for repetitive irradiation of the surface with an area of 5.4 × 5.7 cm2. The scanner produced a laser spot of 11 mm in diameter at the irradiated surface with a top-hat beam profile, and enabled highly homogeneous heating with a repetition rate of 25 Hz. Meanwhile, the surface was actively cooled with chilled air (Cryo 6, Zimmer, Germany). After 120 seconds, the active period was interrupted, flagging the relaxation period. The recorded temporal evolution of the Ts was analyzed (ThermaCam Researcher Pro software 2.10) and used in the data-matching algorithm to estimate the STD, established at the end of the active cycle. Three measurements with each laser system were performed. The average intensity was I=1.2 W/cm2 and was constant for all measurements. In particular, it was found suitable for elevating peak temperatures to 40–47 °C and effectively damaging the adipocytes [6,10]. Colling capacities (CC) were progressed from 50 to 100 W/m2 and were chosen to maintain low surface temperature, which should never exceed 43 °C [33,34]. Measurements were separated by a 30-minute period of tissue relaxation and were repeated once for each experimental condition due to limitations in repeatability of initially distributed temperature field, especially in in-vivo measurements.
The ex-vivo HTLT procedure was performed on a fresh porcine tissue sample (100 × 100 mm2 and 38 mm thick – 2 and 36 mm of skin and fat respectively) as it sufficiently resembles the optical and thermal characteristics of human tissue [25,35]. The porcine sample was thermalized at room temperature and vertically cut in two equal halves prior to measurements. During the active period, the halves were firmly joined. Immediately following the end of the active period, the halves were split, allowing for direct measurement of the reference STD (at the cut cross section). A detailed explanation of the protocol was described in a previous study [9].
Moreover, in-vivo HTLT was performed on the lumbar and iliac abdomen regions of a 28-year-old healthy individual with Fitzpatrik skin type II. Informed consent was obtained from the volunteer and the procedure was performed in the Laser and Health Academy (Ljubljana, Slovenia), by a clinical expert in laser medicine. The study was performed according to the Declaration of Helsinki, and approved by the Slovenian national ethics committee (application number 111/02/12).
2.3 Details on the estimation method
The TPE model and the STD estimation algorithm can be coupled to support an increasing demand of personalized therapeutic procedures. We combined some of the most used laser wavelengths in medicine to estimate the tissue parameters and successively predict the STD within the tissue. Both approaches include the same digital-tissue model, yet the TPE algorithm includes both periods (laser-on and laser-off), whereas the STD algorithm only simulates the relaxation period (laser-off).
Structurally, the tissue model is the same for TPE process and STD algorithm (ex-vivo or in-vivo application), and is described with two planar layers representing skin and fat tissue, whereby a skin layer combines both the dermis and epidermis layers and is 2 mm thick [36]. The thickness of the fat layer was set to 36 mm, simulating an obese person [37] and also matching the thickness of the ex-vivo tissue sample. The photo-thermal interactions diminish with increased tissue thickness and have no major thermal effects at locations deeper than ∼25 mm [38], therefore the muscle layer below the fat tissue was not included in the model.
The general model for temperature field evolution is described with Eq. (1) and has been modified, accounting for the two different sets of measurements (in-vivo, ex-vivo). The ex-vivo model excludes temperature-dependent blood perfusion and metabolic heat generation – the second and third terms of Eq. (1). Due to homogeneous heating of a large surface area, the temperature field evolution can be described with a 1D heat problem, written as follows [7,36,39]:
In the blood perfusion term the cp,b (3800 J/kgK [40]) represents specific heat of blood, ωb blood perfusion rate (kg/m3s) and Tb blood temperature (Tb=37 °C [40]). When computing Eq. (1), different blood perfusion rates were applied to a particular structural layer – ωskin and ωfat for skin and fat layer respectively. The vascular system is temperature-dependent below a certain threshold temperature (TF), and its response to temperature variation can be simulated with a scaling function, introduced by Drizdal et al. [41]:
Where A (10 for skin layer and 1 for fat layer) accounts for temperature dependency of the blood perfusion rate, and B (10 for skin layer and 12 for fat layer) for the steepness of the blood flow change. Once tissue temperature T exceeds critical temperature TF the perfusion rate is equal to the maximum perfusion rate. According to other reports [7,41], the critical temperature for skin was set to TF_skin = 44 °C and for fat tissue to TF_fat = 45 °C. The Qmet represents metabolic heat generation (W/m3).
The interface boundary condition in (Eq. (3)), applied at the air-tissue boundary, governs the heat losses from natural convection and body radiation as follows:
For the TPE model we used an Er:YAG laser source that is implemented into the model by adding a surface-originating heat source qlaser (W/m2) in Eq. (3). This constant irradiation intensity is kept for the entire active period (laser-on). The approach adequately describes laser-tissue interaction due to long illumination times (30 seconds) and when the heat source is localized on the surface (negligible penetration depth of Er:YAG [22], 1/µ ≈ 3µm). The peak temperature surface was kept below 60 °C to avoid structural phase changes or water-complex dissociation [31,32] that would complicate the model, and are not required for the estimation of the tissue thermal parameters.
The TPE algorithm with surface heat source (Eq. (1), (2) and (3)) were applied to each layer of the tissue model. For the skin-fat interface, a boundary condition must be implemented according to the “conductivity matching condition” for finite difference scheme calculations of a heat transfer problem [44]. When computing Eq. (1), the following condition was applied at the skin-fat boundary point:
To run the TPE model, the initial temperature field distribution is required. For the ex-vivo experiment this initial profile was assumed constant along the z-axis and equal to the environmental temperature, due to thermalization of the sample at room temperature prior to measurements. On the other hand, for the in-vivo measurements the initial thermal profile (in equilibrium) was estimated by running a diffusion-based model (Eq. (1)) excluding an external heating source qlaser from the boundary condition in Eq. (3), until the Ts was stabilized. Hence, this initial distribution encompasses the natural gradient in a human subject. For this purpose, the inner temperature of the body was set to 37 °C and metabolic heat generation Qmet to 400 W/m3 [40,45].
Once the tissue parameters were estimated, we ran the custom developed STD algorithm (National Instruments LabVIEW 2020, Austin, TX) to model the relaxation period of the HTLT therapy and estimate the STD within the tissue. The STD algorithm included three main parts: STD generation, calculations of STDs time evolutions and time-dependent data matching process. Similar to a homogeneous fat tissue study [9], the algorithm started with STD generation, whereby an arbitrary function emulated the temperature distributions at the end of the active period. The arbitrary function chosen for this work has been modified to take into account the more complex shape of the profile observed during the validation process of the ex-vivo studies. The first point (at z=0) was constant and was equal to the surface temperature at the end of the active cycle. The final point of interpolation was located at 38 mm from the tissue surface. It was set to room temperature in the ex-vivo experiment, and to the core body temperature of 37 °C for measurements on a human subject. The two points in the middle were used to outline temperature distribution within the tissue. Both points were manipulated within a spatial (zmin, zmax) and temperature (Tmin, Tmax) ranges, whereby the step of each incremental change governed the resolution of the database (in our study 1,835 STDs were created for each measurement). At each iteration the four guiding points were interpolated with a cubic Hermite interpolation.
The generated STDs were clustered in a database and were used in the second part of the algorithm as initial conditions. Here, generated STDs were fed to the physical model (Eq. (1)–(6)), which included higher structural and physiological complexities of the tissue (e.g., blood perfusion, metabolic heat generation and multilayered composition). Presented bioheat-diffusion-equation was solved numerically using the finite difference method implemented in C++ (Microsoft Visual Studio 2019, Redmond, WA) and included in a LabView software (National Instruments LabVIEW 2020, Austin, TX) as a .DLL file. For purpose of temperature evolution calculations, the STDs were discretized with a spatial step of 0.15 mm to satisfy the requirements of numerical stability. The temporal STD evolutions were calculated with a time step of 0.1 s.
The third part of the algorithm included the time-dependent data matching process. Here, the measured Ts evolution was analyzed (ThermaCam Researcher Pro 2.1) and served as an input to the data matching process. The input was correlated with a time-dependent database using residual sum of squares (RSS). The best fit to Ts enabled the extraction of the final STD estimation within the tissue.
3. Results
3.1 Tissue parameters estimation
To estimate tissue parameters (TPE data set), three reference measurements were performed, each with progressively increased laser intensity. The surface temperature response to laser irradiation and the following thermal relaxation are shown in Fig. 2 and reflect the behavior already reported in other studies [23,32,46]. Measured surface temperature Ts dynamics are depicted with colored lines in Fig. 2(a) and 2(b) for ex-vivo and in-vivo, respectively. For both cases, the temporal dynamics of Ts were reproduced well with the TPE model (Er:YAG) presented in the section “Details on the estimation method”. The fit of the digital model is depicted with black dashed lines (Fig. 2(a) and 2(b)).
Fig. 2. Measured Ts evolution (colored lines) during the heating and relaxation cycle of the TPE process: The closest fitting emulations from the digital model are shown with black dashed lines. Three different average intensities were used in the heating phase I1=0.08, I2=0.1 and I3=0.12 W/cm2. Graph a. shows the Ts dynamic on an ex-vivo tissue sample, and b. the Ts dynamic for an in-vivo case. Graphs b. and c. show influence of blood perfusion rates on goodness of fit (RSS) between TPE model and measured Ts dynamic for ex-vivo and in-vivo tissues respectively.
The influence of blood perfusion rates (ω) on goodness of fit between TPE model and measured Ts was analyzed at calibrated tissue’s thermal parameters. The residual sum of squares (RSS) of the fit is shown in Fig. 2(c) and 2(d) for ex-vivo and in-vivo tissues respectively. The ωskin and ωfat were varied within the range of values, while the ratio between skin and fat perfusion rates was maintained (ωskin=3ωfat) [7]. All three TPE calibrations of ex-vivo sample showed best fit (see Fig. 2(c), RSS≈0) with both ω set to zero (Tb was set to bulk temperature of ex-vivo tissue). However, minimal RSS of the fit for in-vivo was obtained with increased ω, further validating the blood perfusion model (Eq. (1) and (2)). In agreement, optimal blood perfusion rates correlate well with the range of values reported in literature (see Table 1).
Table 1. The outcome of the tissue-parameter estimation (TPE) protocol for the calibration of thermal properties of ex-vivo and in-vivo tissues.a
In-vivo Ts (Fig. 2(b)) are generally higher compared to the ex-vivo sample (Fig. 2(a)), due to higher temperatures of the initial STD. Tissue thermal parameters estimated from the TPE fit are listed in Table 1. and correlate well with the range of reported values from literature [8,28,48] and [6,7,21,23,36,41], listed below Table 1. Correlation between ex-vivo and in-vivo thermal parameters suggest the features of porcine and human tissue to be similar, as already recognized by many studies [25,35,47].
3.2 Estimation of spatial temperature distribution
Ex-vivo
The spatial temperature distribution (STD) for Alexandrite and Nd:YAG laser excitation was estimated starting from the relaxation period of the HTLT thermal signal, depicted with colored lines in Fig. 3(a) and 3(c). For STD estimation, we coupled the thermal-diffusion-based model, two-layer tissue structure (see section “Details on the estimation method”) and the data-matching approach reported in our previous publication for a homogeneous tissue [9]. For this analysis we attributed optimized tissue parameters (from TPE–Er:YAG calibration measurements) to the corresponding layer of the digital tissue model (see section “Details on the estimation method”).
Fig. 3. Fitting of Ts evolution and prediction of the STD in porcine tissue: Red color denotes CC1, orange CC2 and blue CC3. Graphs a. and c. show the measured surface temperature (colored lines) and fits (black lines) for Nd:YAG and Alexandrite, respectively. Graphs b. and d. show comparison between the measured STD (colored lines) and estimated distribution (black lines) for Nd:YAG and Alexandrite, respectively.
Black lines in Fig. 3(a) and 3(c) show the best fit to the measured surface temperature dynamic during the relaxation period. The estimated STDs (black lines in Fig. 3(b) and 3(d)) reveal peak temperatures reaching lower values with an increased cooling capacity (CC). The trend was observed in similar studies [8,9] and was further confirmed with simultaneous measurements of STDs at the cross section, depicted as colored lines in Fig. 3(b) and 3(d) for Nd:YAG and Alexandrite, respectively.
In addition, an increasing CC leads to a progressively decreased Ts, while the relaxation peaks become more pronounced (see Fig. 3(a) and 3(c)), indicating a larger temperature amplitude within the tissue (TA = Tmax-Ts). This was corroborated by directly measuring the temperature distribution at the cross-section. The comparison of the reference and estimated STD (Fig. 3(d) and 3(b)) revealed similar peak temperatures and positions, as the average difference was ∼0.35 °C and ∼0.2 mm, respectively. Similarly, the temperature distributions within the tissue were well reproduced, and are comparable with studies on porcine tissue and on human subjects [7–10].
In-vivo
We move the focus of our study to an in-vivo case to investigate the applicability of our approach for the STD estimation when only the Ts evolution is given. Starting from the estimated tissue parameters (TPEs) reported in Table 1., we applied the data-matching algorithm and bioheat transfer equation (see section “Details on the estimation method”) to find the best fitting Ts and the corresponding STD.
The Ts evolutions and corresponding fits (during the relaxation period) obtained with Nd:YAG and Alexandrite laser excitation and different cooling capacities are shown in Fig. 4(a), c. The STDs estimated with an in-vivo model are displayed in Fig. 4(b), d for Nd:YAG and Alexandrite, respectively. The peak’s location was estimated to range 4–6 mm from the irradiated surface. This particular range is consistent to the one obtained in the skin-fat sample and is comparable to a study in human subjects [7], in which the temperature field within the tissue was simulated by a weighted Monte Carlo Multi-Layered (MCML) approach.
Fig. 4. Fitting of Ts evolution and prediction of the STD within a human subject: Red color denotes CC1, orange CC2 and blue CC3. Graphs a. and c. measured surface temperature (colored lines) and fits (black lines) for Nd:YAG and Alexandrite, respectively. Graph b. shows the estimated STDs for Nd:YAG, together with the MCML outcome for the irradiation time of ∼80 s. Graph d. Shows STD estimations for an Alexandrite-based HTLT.
The MCML outcome correlates well with our STD estimations in regard to temperature distributions, temperature peak values and locations. Indeed, our estimations show peak temperatures in a range from ∼40–45 °C after 2 minutes of active period, and supports the estimated range between 43 and 47 °C, reported in [7]. A typical MCML prediction of the temperature field is depicted with a black dotted line in Fig. 4(b). In fact, our estimations of STDs resemble the general shape, however, the MCML outcome predicts the distribution closer to the irradiated surface. This is a consequence of a shorter active period (∼80 s) of the MCML adopted from [7], and thus a shorter time of thermal diffusion within the tissue.
4. Discussion
The STD estimation algorithm relies on the relaxation period of the HTLT procedure, thus the TPE process was focused only on calibration of the tissue’s thermal parameters. In particular, we employed an Er:YAG laser with a wavelength 2,940 nm due to several key factors: (1) at wavelengths shorter than 1,400 nm the scattering becomes dominant over absorption [24], hence a complex light-tissue simulation would be required (e.g., Monte Carlo photon multi-layer procedure) [7,8]. This procedure includes additional temperature-dependent [28] optical parameters that would need to be calibrated (e.g., absorption coefficient, reduced scattering coefficient, anisotropy factor and refraction index) [7,8]. (2) Above 1,400 nm light scattering in tissue becomes insignificant compared to absorption [24]. Moreover, these wavelengths are primarily absorbed in water content instead of melanin, which decreases TPE dependency on melanin volume fraction that varies with skin type (1–10%) [7]. (3) Particular wavelength of Er:YAG (2,940 nm) matches absorption peak of water, hence its absorption is significantly higher compared to Er:YSGG (2,780 nm) [49] and 10-times higher compared to CO2 (10,600 nm) [22,50]. The corresponding optical penetration depth (3 µm [22]) enables implementation of the Er:YAG laser light as a simplified surface originating heat source in a boundary condition of bioheat-diffusion-equation (Eq. (3)), greatly reducing complexity of the TPE process.
Since to optical penetration depth of Er:YAG is smaller than thickness of the skin layer (∼2 mm), the duration of TPE process has to be sufficient to ensure adequate thermal penetration depth. Upon this condition, the predicted thermal penetration depths (TPE model) were 6 mm and 18 mm at the end of the active and relaxation period respectively. As a consequence, both layers significantly affected the temporal evolution of Ts.
The average intensities of TPE process were determined to sufficiently elevate surface temperature within 30 s. However, high temperatures were avoided, due to other effects becoming major when Ts exceeds 65 °C (e.g., excessive water evaporation and structural changes of the tissue [31,32]). In fact, the peak surface temperatures remained below 50 °C for both measured sets (see Fig. 2(a) and 2(b)). Despite progressively higher intensities (see Fig. 2(a) and 2(b)), all measurements were well reproduced with similar thermal parameters. The latter indicates to a non-significant temperature dependance of thermal parameters in particular temperature range. This is in agreement with [24], where a relatively small change of Dfat (-0.15%/°C) was reported. Average standard deviations of estimated thermal parameters were low for both measured sets (∼2–3%, see Table 1), suggesting that single measurement (e.g., I1) sufficiently covers calibration in particular temperature range. In fact, the highest sensitivity to blood perfusion rate was observed at lowest intensity and can thus be estimated with higher accuracy (see Fig. 2(d)).
In ex-vivo experiments a systematic difference in peak temperature location zmax was observed when compared to our previous study performed on a homogeneous layer of fat [9]. In fact, under the same experimental conditions, the zmax was located ∼1–2 mm deeper compared to the sample with an additional layer of skin. This blocking effect of the skin was described by Han et al. [38], and was ascribed to the smaller attenuation length of the skin for IR lasers. In this regard, a thermal attenuation length was estimated for fat-only and skin-fat samples as 5.1 ± 0.4 mm and 4.7 ± 0.7 mm, respectively. Hence, due to a skin layer blocking effect, less photons penetrate deeper within the tissue [38], resulting in a temperature peak closer to the irradiated surface.
Hyperthermia-based treatments of larger thermal mass and smooth thermal gradients require lower spatial (2-3 mm) and temporal resolution [17]. However, the required temperature accuracy is less than 1 °C for regimes below boiling point [17]. Ex-vivo estimations were corroborated by direct measurements of STD (see Fig. 3(b) and 3(d)) and revealed an average error smaller than ±0.5 mm and ±0.5 °C for the peak’s location and value respectively. In-vivo estimations of STDs resemble a general distribution when compared to a similar MCML study. However, in future in-vivo investigations we aim to include multiple subjects and additional validation of the method by 3D temperature monitoring in real-time such as volumetric reconstruction of the optoacoustic (OA) signal (accuracy with a deviation of 10% on the relative temperature [13]) or ultrasound [16,17].
The STD algorithm is potentially limited for real-time temperature monitoring, since the estimation of temperature field is given after the active period. However, this concept could be employed during a long lasting clinical procedures (∼25–30 min) [7,10] that require an incremental heating, interrupted by multiple relaxation periods. In this temporal window the method could provide an estimation of STD, thus the efficiency of the next active cycle would be better addressed. Moreover, the focus on relaxation period enhances personalization and versatility of the method. The latter is achieved with method’s independency of heating or cooling systems. In fact, simulations focused on active period are therapy-tailored and depend on many parameters such as tissue light absorption and scattering, cooling system capacity, and initial STD correctness. Since only thermal parameters are required for STD algorithm, a high and fast personalization could be achieved with a thermal-parameters-oriented calibration approach (TPE model).
The presented approach is beneficial in treatments where a larger areas are processed. Due to its non-contact nature, it avoids the insertion of thermocouples within the skin that have limited spatial resolution and can induce discrepancy between the probed area and the temperature of the surrounding tissue [11,17]. However, 1D approach is only valid if irradiated area is relatively large and if temperature fields are laterally homogeneous. Indeed, this is satisfied within the area of irradiated skin, however it is limited at the boundaries of the area, due to lateral heat diffusion distance (d=(4Dtir)1/2) being 7 and 8 mm for skin and fat layers respectively (at irradiation time 120 s and thermal diffusivities as listed in Table 1.). As currently presented, the STD algorithm assumes heterogeneous structural layers, however, the method is not applicable for treatments of highly localized areas (e.g., tattoo ink, tumor or vessel ablation).
5. Conclusion
We presented a novel approach for estimation of the tissue’s thermal parameters (TPE) and estimation of the spatial temperature distribution (STD) within tissue after laser irradiation. The TPE model accurately reproduces the surface temperature dynamics the tissue thermal parameters with an average relative standard deviation less than 3%. The latter enables the tissue-tailored approach of the STD estimation algorithm that was applied to ex-vivo and in-vivo tissues during a Hyperthermic Laser Treatment. The validity of the protocol was corroborated by direct measurement of the STD across the section of the ex-vivo sample, which also reveals a good correlation of the temperature distributions. In this regard, an average error was found to be smaller than ±0.5 mm and ±0.5 °C for the peak’s location and value respectively. The estimated STD during the in-vivo investigations showed a similar trend to other reported studies. The combination of the TPE calibration protocol and the wavelength-independent data matching algorithm represent a valid tool for the estimation of the STD, offering the possibility to create personalized database and laser treatment such as selective photothermolysis. The presented method can be readily implemented with integrated thermal cameras on medical devices.
Funding
Ministrstvo za Izobraževanje, Znanost in Šport (GOSTOP program, contract number C3330-16-529000); Javna Agencija za Raziskovalno Dejavnost RS (core funding No. P2-0392, project funding No. L2-1833).
Acknowledgments
Authors would like to thank Fotona d.o.o. for providing commercially available medical laser system. Authors thank dr. Bojan Bozic, Institute of Biophysics, University of Ljubljana for help in mathematical modeling, and dr. Jernej Kukovic, M.D., Clinical consultant in Laser Medicine, Laser and Health Academy, Ljubljana for help in experiments performed on human subject.
Disclosures
The authors declare that there are no conflict of interest.
Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
References
1. K. Nouri, Lasers in Dermatology and Medicine (Springer Science + Business Media B.V., 2011).
2. L. Carroll and T. R. Humphreys, “LASER-tissue interactions,” Clin. Dermatol. 24(1), 2–7 (2006). [CrossRef]
3. D. Simhon and A. Katzir, “Soldering System and Its Clinical Application for Bonding Skin Incisions,” J. Biomed. Opt. 20(12), 128002 (2015). [CrossRef]
4. K. F. Chu and D. E. Dupuy, “Thermal ablation of tumours: Biological mechanisms and advances in therapy,” Nat. Rev. Cancer 14(3), 199–208 (2014). [CrossRef]
5. H. H. Zenzie, G. B. Altshuler, M. Z. Smirnov, and R. R. Anderson, “Evaluation of cooling methods for laser dermatology,” Lasers Surg. Med. 26(2), 130–144 (2000). [CrossRef]
6. W. Franco, A. Kothare, S. J. Ronan, R. C. Grekin, and T. H. McCalmont, “Hyperthermic injury to adipocyte cells by selective heating of subcutaneous fat with a novel radiofrequency device: Feasibility studies,” Lasers Surg. Med. 42(5), 361–370 (2010). [CrossRef]
7. M. Milanic, B. T. Muc, N. Lukac, and M. Lukac, “Numerical Study of Hyper-Thermic Laser Lipolysis With 1,064 nm Nd:YAG Laser in Human Subjects,” Lasers Surg. Med. 51(10), 897–909 (2019). [CrossRef]
8. M. Milanic, B. T. Muc, M. Jezersek, and M. Lukac, “Experimental and numerical assessment of hyperthermic laser lipolysis with 1,064 nm Nd:YAG laser on a porcine fatty tissue model,” Lasers Surg. Med. 50(2), 125–136 (2018). [CrossRef]
9. J. Kosir, D. Vella, and M. Jezersek, “Non-contact monitoring of the depth temperature profile for medical laser scanning technologies,” Sci. Rep. 10(1), 20242–10 (2020). [CrossRef]
10. J. W. Decorato, B. Chen, and R. Sierra, “Subcutaneous adipose tissue response to a non-invasive hyperthermic treatment using a 1,060 nm laser,” Lasers Surg. Med. 49(5), 480–489 (2017). [CrossRef]
11. P. Saccomandi, E. Schena, and S. Silvestri, “Techniques for temperature monitoring during laser-induced thermotherapy: An overview,” International Journal of Hyperthermia 29(7), 609–619 (2013). [CrossRef]
12. I. V. Larina, K. V. Larin, and R. O. Esenaliev, “Real-time optoacoustic monitoring of temperature in tissues,” J. Phys. D: Appl. Phys. 38(15), 2633–2639 (2005). [CrossRef]
13. F. J. O. Landa, X. L. Deán-Ben, R. Sroka, and D. Razansky, “Volumetric Optoacoustic Temperature Mapping in Photothermal Therapy,” Sci. Rep. 7(1), 1–8 (2017). [CrossRef]
14. V. Rieke and K. B. Pauly, “MR thermometry,” J. Magn. Reson. Imaging 27(2), 376–390 (2008). [CrossRef]
15. L. Liu, S. Creten, Y. Firdaus, J.J. Agustin Flores Cuautle, M. Kouyaté, M. Van Der Auweraer, and C. Glorieux, “Fluorescence spectra shape based dynamic thermometry,” Appl. Phys. Lett. 104, 031902 (2014). [CrossRef]
16. J. Shah, S. R. Aglyamov, K. Sokolov, T. E. Milner, and S. Y. Emelianov, “Ultrasound imaging to monitor photothermal therapy – Feasibility study,” Opt. Express 16(6), 3776 (2008). [CrossRef]
17. M. A. Lewis, R. M. Staruch, and R. Chopra, “Thermometry and ablation monitoring with ultrasound,” International Journal of Hyperthermia 6736, 163–181 (2015). [CrossRef]
18. L. Yin, M. S. R. Gudur, Y. S. Hsiao, R. E. Kumon, C. X. Deng, and H. Jiang, “Tomographic reconstruction of tissue properties and temperature increase for high-intensity focused ultrasound applications,” Ultrasound Med. Biol. 39(10), 1760–1770 (2013). [CrossRef]
19. T. E. Milner, D. M. Goodman, B. S. Tanenbaum, and J. S. Nelson, “Depth profiling of laser-heated chromophores in biological tissues by pulsed photothermal radiometry,” J. Opt. Soc. Am. A 12(7), 1479 (1995). [CrossRef]
20. L. Jiang, W. Zhan, and M. H. Loew, “A numerical study of the inverse problem of breast infrared thermography modeling,” Med. Imaging 2010 Biomed. Appl. Mol. Struct. Funct. Imaging 7626, 76260O (2010). [CrossRef]
21. S. R. Mordon, B. Wassmer, J. P. Reynaud, and J. Zemmouri, “Mathematical modeling of laser lipolysis,” BioMed Eng OnLine 7(1), 10–13 (2008). [CrossRef]
22. B. Majaron, P. Plestenjak, and M. Lukač, “Thermo-mechanical laser ablation of soft biological tissue: Modeling the micro-explosions,” Applied Physics B: Lasers and Optics 69(1), 71–80 (1999). [CrossRef]
23. P. Barvinschi and O. M. Bunoiu, “A finite element method model to simulate the thermal response of vascularized tissue to laser irradiation,” Rom. Reports Phys.69, (2017).
24. A. J. Welch and M. J. C. Van Gemert, Optical-Thermal Reponse of Laser-Irradiated Tissue, 2nd ed. (Springer Science + Business Media B.V., 2011). [CrossRef]
25. A. N. Bashkatov, E. A. Genina, and V. V Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: A review,” J. Innov. Opt. Health Sci. 04(01), 9–38 (2011). [CrossRef]
26. S.-H. Tseng, P. Bargo, A. Durkin, and N. Kollias, “Chromophore concentrations, absorption and scattering properties of human skin in-vivo,” Opt. Express 17(17), 14599 (2009). [CrossRef]
27. M. M. Paulides, P. R. Stauffer, E. Neufeld, P. F. MacCarini, A. Kyriakou, R. A. M. Canters, C. J. Diederich, J. F. Bakker, and G. C. Van Rhoon, “Simulation techniques in hyperthermia treatment planning,” International Journal of Hyperthermia 29(4), 346–357 (2013). [CrossRef]
28. S. Kim and S. Jeong, “Effects of temperature-dependent optical properties on the fluence rate and temperature of biological tissue during low-level laser therapy,” Lasers Med. Sci. 29(2), 637–644 (2014). [CrossRef]
29. J. Royo, J. Moreno-Moraga, and M. A. Trelles, “Clinical assessment of a new 755 nm diode laser for hair removal: efficacy, safety and practicality in 56 patients,” Lasers Surg. Med. 49(4), 355–360 (2017). [CrossRef]
30. A. Ayatollahi, A. Samadi, A. Rajabi-estarabadi, S. Yadangi, K. Nouri, and A. Firooz, “Comparison of efficacy and safety of a novel 755-nm diode laser with conventional 755-nm alexandrite laser in reduction of axillary hairs,” Lasers Med. Sci. 35, 373 (2019). [CrossRef]
31. J. H. Torres, M. Motamedi, J. A. Pearce, and A. J. Welch, “Experimental evaluation of mathematical models for predicting the thermal response of tissue to laser irradiation,” Appl. Opt. 32(4), 597 (1993). [CrossRef]
32. E. N. Sobol, M. S. Kitai, N. Jones, A. P. Sviridov, T. Milner, and B. J. F. Wong, “Heating and structural alterations in cartilage under laser radiation,” IEEE J. Quantum Electron. 35(4), 532–539 (1999). [CrossRef]
33. D. Yarnitsky, E. Sprecher, R. Zaslansky, and J. A. Hemli, “Heat pain thresholds: normative data and repeatability,” Pain 60(3), 329–332 (1995). [CrossRef]
34. L. Arendt-Nielsen and A. C. N. Chen, “Lasers and other thermal stimulators for activation of skin nociceptors in humans,” Neurophysiologie Clinique/Clinical Neurophysiology 33(6), 259–268 (2003). [CrossRef]
35. N. J. Vardaxis, T. A. Brans, M. E. Boon, R. W. Kreis, and L. M. Marres, “Confocal laser scanning microscopy of porcine skin: implications for human wound healing studies,” J. Anat. 190(4), 601–611 (1997). [CrossRef]
36. P. Wongchadakul, P. Rattanadecho, and T. Wessapan, “Implementation of a thermomechanical model to simulate laser heating in shrinkage tissue (effects of wavelength, laser irradiation intensity, and irradiation beam area),” Int. J. Therm. Sci. 134, 321–336 (2018). [CrossRef]
37. R. A. D. Booth, B. A. Goddard, and A. Paton, “Measurement of fat thickness in man: a comparison of ultrasound, Harpenden calipers and electrical conductivity,” Br. J. Nutr. 20(4), 719–725 (1966). [CrossRef]
38. D. Han, J. Xu, Z. Wang, N. Yang, X. Li, Y. Qian, G. Li, R. Dai, and S. Xu, “Penetrating effect of high-intensity infrared laser pulses through body tissue,” RSC Adv. 8(56), 32344–32357 (2018). [CrossRef]
39. H. H. Pennes, “Analysis of tissue and arterid blood temperatwes,” J. Appl. Phys. 1, 20 (1948). [CrossRef]
40. A. R. Melo, M. M. S. Loureiro, and F. Loureiro, “Blood perfusion parameter estimation in tumors by means of a genetic algorithm,” Procedia Computer Science 108, 1384–1393 (2017). [CrossRef]
41. T. Dřždǎl, P. Togni, L. Víšek, and J. Vrba, “Comparison of constant and temperature dependent blood perfusion in temperature prediction for superficial hyperthermia,” Radioengineering 19, 281–289 (2010).
42. B. W. Scott and A. S. Vance, “A tissue heat transfer model relating dynamic skin temperature changes to physiological parameters,” Phys. Med. Biol. 486, 475–486 (1980). [CrossRef]
43. Z. S. Deng and J. Liu, “Mathematical modeling of temperature mapping over skin surface and its implementation in thermal disease diagnostics,” Comput. Biol. Med. 34(6), 495–521 (2004). [CrossRef]
44. R. I. Hickson, S. I. Barry, G. N. Mercer, and H. S. Sidhu, “Finite difference schemes for multilayer diffusion,” Mathematical and Computer Modelling 54(1-2), 210–220 (2011). [CrossRef]
45. P. Wongchadakul, P. Rattanadecho, and T. Wessapan, “Simulation of temperature distribution in different human skin types exposed to laser irradiation with different wavelengths and laser irradiation intensities,” Songklanakarin J. Sci. Technol. 41, 529–538 (2019).
46. H. Z. Alagha, H. Z. Alagha, and M. Gülsoy, “Photothermal ablation of liver tissue with 1940-nm thulium fiber laser: an ex vivo study on lamb liver thulium fiber laser: an ex vivo study on lamb liver,” (1940).
47. F. A. Duck, Physical Properties of Tissue: A Comprehensive Reference Book (Academic Press, 1990).
48. M. A. El-Brawany, D. K. Nassiri, G. Terhaar, A. Shaw, I. Rivens, and K. Lozhken, “Measurement of thermal and ultrasonic properties of some biological tissues,” J. Med. Eng. Technol. 33(3), 249–256 (2009). [CrossRef]
49. J. Diaci, C. Filipic, T. Perhavec, and M. Lukac, “Influence of water absorption shift on ablation speed of Er:YAG and Er, Cr:YSGG dental lasers,” J. Laser Heal. Acad., 17–22 (2013).
50. M. Lukac, T. Perhavec, K. Nemes, and U. Ahcan, “Ablation and thermal depths in VSP Er:YAG laser skin resurfacing,” J. Laser Heal. Acad. 1, 56–71 (2010).
