Effect of optical energy modulation on the thermal response of biological tissue: computational and experimental validations.

This study develops an energy modulation technique to attain a constant interstitial tissue temperature and to induce the predetermined thermal coagulation without carbonization in tissue. An optical diffuser was employed to deliver 1064 nm light to the biological tissue. The combined mode maintained the interstitial temperature at 70 °C for longer durations compared to the continuous wave mode. Coagulation volumes increased linearly with the time and met the predetermined treatment volume range (0.32-0.52 cm3) after the combined treatment for 100 s. The combined modulation can be a feasible modality to induce the predetermined extent of thermal coagulation for treating papillary thyroid microcarcinoma.


Introduction
Papillary thyroid carcinoma with a diameter less than 8.5 mm (volume ≤ 0.32 cm 3 ) is referred to as papillary thyroid microcarcinoma (PTMC), which accounts for approximately 50% of all papillary thyroid carcinoma cases [1,2]. Almost PTMCs are either clinically silent or sometimes found in autopsies of non-thyroid-related patients [3]. Typically, PTMCs can be treated via open surgery. The remaining cases include patients that are either unwilling to undergo invasive resection or are at a high risk for open surgery [2]. Moreover, surgical complications such as damage to laryngeal nerve and parathyroid gland may be more dangerous than the disease itself [4]. For such patients, several non-surgical thermal treatments, including radiofrequency ablation (RFA), microwave ablation (MWA), and photothermal coagulation (PTC), have been proposed [4,5].
Among various thermal treatments, ultrasound-guided PTC is a promising option yielding high precision, minimal invasiveness, and fast recovery; it is also relatively less painful than traditional open surgery [4,5]. PTC uses an optical fiber to induce heat to damage neoplastic tissue. Currently, PTC has been widely deployed to treat liver, brain, pancreas, prostate, and lung tumors and cancers [5]. The efficacy of PTC depends on tissue properties (absorption, scattering coefficients, and anisotropy factor) and laser light characteristics (wavelength, spot size, power, and irradiation time) [6]. In principle, tissues are irreversibly coagulated at approximately 60 ℃, and heat diffusion around the irradiated tissue causes supraphysiological hyperthermia at 42-60 ℃, inducing delayed thermal damage (24-72 h) [7,8]. Temperature increases above 100 ℃ induce carbonization and charring effects in the regions surrounding the fiber tip [9,10], resulting in unfavorable complications and a delayed healing response after the treatment [11]. Moreover, the other drawbacks of PTC include difficult fiber placement, non-uniform injury, and mismatching with the target tumor shape, resulting in unexpected side-effects after the PTC treatment [5]. Two groups of therapeutic monitoring techniques can be combined with PTC: invasive methods (i.e., thermocouple, thermistor, and fiber optic sensors) and non-invasive methods (i.e., computed topography, magnetic resonance imaging, and ultrasound imaging) [5,12]. Recently, photoacoustic imaging is considered as a hybrid imaging technology for clinical treatment [13]. This modality uses both optical and acoustic energy to produce high optical absorption contrast while maintaining rich ultrasonic spatial resolution in biological tissues [14]. Previous studies demonstrated a capacity of the photoacoustic imaging for photothemal therapeutic monitoring in various tissues [14].
To manage the thermal damage within target tumors, recent studies have proposed several feasible approaches to monitor and constrain the temperature within a predetermined range. A controllable shutter integrated with both a 1064 nm laser and a thermocouple was used to maintain the temperature at a desired point between 42°C and 44°C by switching the power on and off [15]. However, this study was unable to predict the period of time for which the temperature must be maintained at the desired point to achieve a coagulation zone of the desired size, as no computational modeling was performed. Similarly, a thermistor-based, controlled power modulation system was developed to modulate the energy of an 805 nm laser to maintain the temperature at 50 ℃ using a stepwise power technique [16]. Nevertheless, this Monte Carlo simulation-based model merely predicts the temperature distribution without simulating thermal damage. Subsequently, a model based on a proportional-integrative-derivative (PID) controller was used to maintain a temperature of 80 ℃ at a strategic position within the tissue [11]. However, a manual tuning algorithm was used for the simulation models in this study, whereas the auto-tuning algorithm was used for experiments, leading to inherent inconsistencies between the results obtained using the two approaches. Although feedback-controlled temperature may be useful for surgeons to adjust laser light dosages, PTC is usually implemented without monitoring the temperature in real-time [17]. Recently, we proposed a linear energy modulation model to maintain the temperature in the preset range (65-75 ℃); however, the model did not consider the optical properties of both native and coagulated tissues, resulting in unexpected errors in simulated and experimental outcomes [5]. Therefore, the development of additional models is desirable to simplify the screening step of dosimetry for PTC; this will ultimately enhance treatment outcomes, such as the size of coagulation zones, and ensure that there is no carbonation. The current study modulated 1064 nm light energy in a continuous wave (CW) and combined manner. This was done to maintain the interstitial temperature at 70 ℃, with a deviation of ∼ 3%, for a considerably long duration and to induce the predetermined volume of thermal damage (0.32 ÷ 0.52 cm 3 ) for the treatment of tumors with a diameter of ∼ 8.5 mm (volume ∼ 0.32 cm 3 ) without carbonization. We hypothesized that performing energy modulation could achieve a constant interstitial temperature for a prolonged duration and induce the extent of thermal coagulation to the predefined tissue volume in a controlled manner.

Mathematical algorithm of proposed model
A mathematical algorithm was employed to investigate the optimal mode ( Fig. 1(A)). First, three criteria, namely, the temperature range, treatment volume, and carbonization effect were set. The preset temperature range was 70 ± 2 ℃, which is approximately 10 ℃ higher than the coagulation temperature threshold (∼ 60 ℃) in tissue, to achieve uniform coagulation based upon the previous study [5]. The treatment volume was set within the range of 0.32 cm 3 (i.e., volume of an 8.5 mm PTMC) to 0.52 cm 3 (i.e., diameter ∼ 10 mm). This value of the maximum volume was chosen to ensure the complete coagulation of the predetermined PTMC [2]. To avoid follow-up complications in practical treatment, no visual carbonization was assumed. The mathematical algorithm can be described in several following steps. At first, the CW mode (constant power and build-up time) was employed to achieve the preset treatment volume. Then, the temperature range, treatment volume, and carbonization effect were compared with preset criteria. In case of unsatisfied results, the combined modes with time and power changes were employed to select the optimal mode and to the end of the optimization.

Theoretical model
A theoretical model was developed to predict the thermal responses (temperature characteristics and thermal injury) of the tissue under laser treatment ( Fig. 1(B)). A physical model comprising two components, namely, porcine liver tissue (dimensions = 6 × 6 × 6 cm 3 ; volume = 216 cm 3 ) and an optical diffuser (core diameter = 600 µm; glass cap diameter = 1.2 mm; active length = 10 mm) was considered. A 1064 nm laser source with a long optical penetration depth (∼ 4.8 mm for native tissue and ∼ 1.8 mm for coagulated tissue) was employed to irradiate the tissue [18,19]. Computational analysis was conducted using COMSOL Multiphysics software 5.3 (COMSOL Inc., Burlington, MA, USA) because of its intrinsic multi-physics modules and ease of use [4,5].

Heat transfer in tissue
A simple diffusion approximation approach was employed to mathematically explain the transport of the irradiated light in the tissue using Beer's law [11,20]: where E 0 (W/cm 2 ) is the irradiance at the glass cap and tissue interface, E (W/cm 2 ) is the fluence rate at the tissue depth z (cm), and µ t (cm −1 ) is the total attenuation coefficient. The incident light was absorbed by the tissue chromophores and converted into heat, which then diffused into the surrounding area. These thermal phenomena can be described using Pennes' bio-heat transfer equation [21]: where ρ (kg/cm 3 ), c (J/kg·°C), k (W/m·°C), T (°C), and t (s) are the density, specific heat, thermal conductivity, temperature of the tissue, and irradiation time, respectively (Table 1). Blood parameters of ρ b (kg/m 3 ), c b (J/kg·°C), ω b (m 3 /m 3 ·s), and T b (37℃) were employed as the density, specific heat, perfusion rate per unit volume, and arterial blood temperature, respectively. The terms Q r (W/m 3 ) and Q met (W/m 3 ) denote the heat generated owing to light absorption and metabolic heat generation per unit volume, respectively. Because all experiments were conducted ex vivo, the metabolic heat source (Q met ) and blood perfusion in the tissue (ω b ) were neglected in the numerical simulation [11,22]. The thermal properties of the liver tissue assumed in this model were retrieved from previous studies (Table 1) [11,19]. The initial temperature was assumed to be uniform at 25 ℃.

Table 1. Optical and thermal parameters for modeling of laser-tissue interaction (1064 nm) in porcine liver
Parameter Value References Native Coagulated

Laser light heat source
For simplicity, the proposed model assumed that the laser light diffuses radially along the optical fiber without a forward light component (i.e., cylindrical propagation). In addition, optical scattering at 1064 nm in the liver tissue predominates over absorption [23]; therefore, the effective attenuation coefficient (µ eff ) was used based on a diffusion approximation [9]. The heat source term Q r from light absorption was calculated according to the following equations [11]: where µ a (1/mm), µ s (1/mm), µ eff (1/mm), g, and r (mm) denote the absorption coefficient, scattering coefficient, effective attenuation coefficient of the tissue, anisotropy factor, and radial distance from the glass cap, respectively (Table 1).

Laser power modulation for simulation
During PTC, the power and the therapeutic duration (i.e., the total applied energy) are crucial to managing the degree of irreversible thermal coagulation within the tissue. The constant power required to induce deep necrosis may involve carbonization within the tissue due to overheating. Therefore, to attain a preset temperature (70 ± 2°C) at the point of interest (POI: 2 mm from the optical diffuser axis based on optical penetration depth ∼ 4.8 mm for native tissue and ∼ 1.8 mm for coagulated tissue; Fig. 1(B)) [19], the laser power was modulated with time. The concept of combined energy modulation involves switching off the power supply once the temperature reaches 70 ℃. Subsequently, the temperature will decrease exponentially to the room temperature [5]. To indirectly compensate for this decrease, the power is exponentially modulated to maintain a stable temperature at the predetermined level. The power modulation P e (t) in the exponential period can be calculated as follows: where P eo , a 1 , a 2 , b 1 , b 2 , t, t 1 , and t 2 are the power constant, amplitudes, time constants of the exponential function, exposure time, starting time point, and ending time point for exponential energy modulation, respectively. Based on our previous study [11], 10 W in the CW mode was selected as a pilot power for the current study. After multiple iterative simulations for optimization, the CW mode (10 W for 35 s) and combined mode (constant 10 W for 20 s and exponential decrease from 3.6 W to 2.1 W for following 80 s) were chosen for a quantitative comparison (Figs. 1(A) and 3(A)).

Thermal injury in tissue
Thermal injuries in tissues are dictated by the Arrhenius equation [11]: where , A (1/s), E a (J/mol), and R (J/mol·K) are the dimensionless damage index, preexponential factor, activation energy, and universal gas constant, respectively ( Table 1). The tissue was assumed to experience irreversible damage when = 1 (i.e., log 10 = 0 and 63% of tissue undergoes complete cellular denaturation as the temperature rises to 60 ℃) [11]. Both undamaged (f u ) and damaged (f d ) tissue fractions can be expressed with , as follows [11]: However, the absorption coefficient, scattering coefficient of the tissue, and anisotropy factor in Eqs. (3) and (4) are temperature-dependent, and the thermal damage index is given as follows [11]: where µ a,na , µ a,co , µ s,na , and µ s,co (1/mm), g na , and g co are the native absorption coefficient, coagulative absorption coefficient, native scattering coefficient, coagulative scattering coefficient, native anisotropy factor, and coagulative anisotropy factor of the tissue, respectively. The optical parameters of the porcine liver tissue in this model were retrieved from the previous study (Table 1) [19].

Experimental study
Porcine liver tissue was procured from a local slaughterhouse and stored at -15 ℃. The tissue was sliced, allocated into two sample groups (CW and combined groups; dimensions = 6 × 6 × 6 cm 3 ; volume = 216 cm 3 ; N = 5; Fig. 1(C)), and defrosted to room temperature (25 ℃) before the tests. To validate the simulation results, a laser diode system with a wavelength of 1064 nm (FC-W-1064B-30W, CNI, Changchun, China) was deployed to irradiate the light through the customized optical diffuser (core diameter = 600 µm, active length = 10 mm; TeCure, Inc., Busan, Korea) to treat the porcine liver tissue [24]. The light intensity distribution of the active length in the air was measured as previously described [24]. The delivered energy was modulated by temporally altering the power supply using a customized power modulator, as described in a recent report [5]. During treatment, the temperature was monitored using a thermocouple (type K, Omega Engineering Inc., Norwalk, USA), and the data were acquired at the rate of 10 Hz and sent to a central controller. A customized holder was designed to integrate the optical diffuser with the thermocouple for convenient insertion into the tissue. The tip of the thermocouple was firmly positioned in the middle of the active part of the diffuser. To minimize the direct absorption of the light from the diffuser, the tip of the thermocouple was covered with aluminum foil and placed 2 mm away from the diffuser axis, based on the optical penetration depth of the irradiated light (∼ 4.8 mm for native tissue and ∼ 1.8 mm for coagulated tissue) [19]. A portable ultrasound machine (UGEO HM70A, SAMSUNG Co., South Korea) was used to confirm the position of the optical diffuser and the thermocouple as well as to monitor thermal damage in the tissue in both cross-sectional and longitudinal views (Figs. 1(D) and 1(E)). The ultrasound head of the machine, with a linear transducer (5-13 MHz), acquires images of objects located deep in the tissue.
After the experiments, all treated samples were kept at -15 ℃ for two hours and cross-sectioned in half; then, images of the coagulated areas were captured using a digital camera (D5100, Nikon Corporation, Tokyo, Japan). The coagulated areas were identified as discolored regions (i.e., tan color) in the tissue after treatment and measured using Image J (National Institute of Health, Bethesda, MD) for quantitative comparisons. The coagulation shape was assumed to be an ellipsoid, and its volume (V C ) was calculated according to the following equation [11]: where x, y, and z are the coagulation thicknesses measured in three dimensions.
To address the changes in the longitudinal coagulation shape with time, the ellipticity (E L ) was calculated as follows: Moreover, the eccentricity (e) of the cross-sectional coagulation rim was also evaluated to assess the distribution of thermal coagulation around the optical diffuser. To investigate the degree of the coagulation eccentricity (E C ), the ratio between the eccentricity (e) and the diameter of the glass cap (D C ) was calculated as follows [24]: The Mann-Whitney U test was used to assess the statistical differences between the control and treated samples. The non-parametric method was employed because of the small sample size and degree freedom of the data distribution [25,26]. SPSS (SPSS, Inc., Chicago, USA) was used to analyze the data, and p < 0.05 was considered statistically significant.

Figure 2(A)
illustrates the cross-sectional and longitudinal light distributions arising from the optical diffuser. The cross-sectional distribution indicates a nearly isotropic emission, with less than 5% of the light intensity deviation above 2π. In the longitudinal direction, the digital and goniometric measurements are in good agreement and present an almost cylindrical light distribution along the diffuser axis with less than approximately 20% of the light intensity deviation over the active length. Figure 2(B) and 2(C) show cross-sectional and longitudinal temperature maps, respectively, in the CW and combined modes, along with their temperature distribution profiles (far left). In summary, the CW mode exhibited a sharper temperature distribution gradient arising from the optical diffuser than the combined mode. Across the active length (10 mm) of the optical diffuser, the temperature difference was 39°C in the CW mode, which was 27°C higher than that in the combined mode (∼ 12°C). Figure 3 shows a comparison of the numerical and experimental results obtained for the power modulation, temporal temperature development, laser energy comparison, and treatment time, maintained within preset temperature limits, in the CW and combined modes. Figure 3(A) describes the modulated power used in the current simulation model, whereas Fig. 3(B) depicts the power controlled by the customized modulator in the experiments. In the CW mode, a constant power of 10 W was delivered for 35 s during the treatment. In contrast, in the combined mode, 10 W of power was supplied to rapidly increase the tissue temperature to 70 ℃ for 20 s by irradiation, after which the power was instantly dropped to 3.6 W. Subsequently, the power was exponentially decreased to 2.1 W within 80 s to maintain the desired temperature towards the end of the treatment. Figure 3(C) and 3(D) show the temporal development of the temperature at the POI in the CW and combined modes for both numerical simulations and experimental results. In the simulated CW mode, the temperature drastically increased from room temperature (∼25 ℃) to 91 ℃ (transient increase rate = 1.9 ℃/s), which was equivalent to the experimental result (90 ± 5 ℃). In contrast, the simulated combined mode raised the temperature to 70 ℃ after 20 s of laser exposure. For the remaining 80 s of the treatment duration, the temperature remained constant at 70 ℃ (transient increase rate = 0.5 ℃/s). This was in agreement with the experimental observations. However, after irradiation at 10 W of laser power for 20 s, the combined mode increased the tissue temperature to 67.3 ℃ (transient increase rate = 2.1 ℃/s); this was followed by a slight increase to 72.1 ℃ after further laser application for 80 s. Figures 3(E) and 3(F) compare the treatment conditions (CW and combined modes) based on the treatment time, depending on whether the temperature is maintained within the preset limits, and laser energy delivery.
Overall, the combined mode maintained the predetermined temperature (70 ± 2℃) for a longer period of time at the POI than the CW condition. The temperature was only maintained for approximately 15 s (increased for 3 s and decreased for 12 s, respectively), within the predetermined temperature range, in the CW condition. On the other hand, the combined mode maintained the appropriate temperature for 84 s and 88 s in the simulation and experiment, respectively (p < 0.01). Although the combined mode maintained the preset temperature for a longer duration, its total irradiated energy (395 ± 15 J) was only 45 J higher than that of the CW mode (350 ± 10 J; p < 0.05). Figure 4 shows the longitudinal assessment of thermal coagulation, via simulation and experiments, after laser irradiation under the two irradiation conditions (CW and combined modes), along with the corresponding ultrasound images. The CW mode exhibits an elliptical rim in the coagulative region, whereas a circular shape is observed in the combined mode (Figs. 4(A)  and 4(B)). Figure 4(C) depicts the progressive occurrence of a hyperechoic (bright) region with micro water bubbles (temperature ∼ 100℃) around the diffuser tip in the CW mode, whereas a homogeneous echo region (temperature ∼ 70℃) was observed in the combined mode. The volumes of the thermally injured tissue were quantitatively measured in the CW mode, yielding 0.41 cm 3 and 0.44 cm 3 for the simulation and experiment, respectively (Fig. 4(E)). However, carbonization was observed in the CW mode at a volume of 0.04 cm 3 (i.e., ∼ 9.5% of the total injury volume; Fig. 4(D)). For the combined mode, a coagulated tissue volume of 0.39 cm 3 was predicted via the simulation, and 0.42 cm 3 of the tissue was coagulated in the experimental study without carbonization (p = 0.55). Figure 5 presents the increase in the coagulated volumes and shape factors with time in the combined mode, obtained from simulations and experiments ( Table 2). The simulation and experiments yielded comparable coagulation volumes regardless of the treatment time. The ellipticity (E L ) of the longitudinal coagulation profiles decreased from 87% to 54%, 38%, 21%, and 15% in the simulation and from 82% to 46%, 42%, 25%, and 20% in the experiments (i.e., the profile became almost circular). The eccentricity (E C ) of the cross-sectional coagulation profiles was almost below 11%. The coagulation volumes increased linearly with the treatment time. After a 20 s treatment, the simulated coagulation volume was 0.03 cm 3 (coagulation increase rate = 0.0015 cm 3 /s); the volume changed to 0. 16

Discussion
The significance of the current study lies in the advancement of optical energy modulation technique in order to maintain the constant interstitial tissue temperature with less fluctuation and to achieve the predictable and precise photothermal treatment. Our previous study demonstrated a linear energy modulation model to attain the tissue temperature in the predetermined range. However, the linear energy modulation still suffered from a wide deviation of ∼ 10 ℃ in the temperature. In addition, the previous study merely reflected the optical properties of the native tissue in the simulation, leading to a substantial difference between numerical and experimental results [5]. On the other hand, the current study proposed a non-linear energy modulation to achieve more constant tissue temperature with a deviation of ∼ 3%. The optical properties of both native and coagulated tissues were reflected in the current simulation model, which showed a good agreement with the ex vivo experimental results. The current study was performed because PTC usually poses the risk of carbonization and charring effects in the regions surrounding the fiber tip due to overheating within the tissue, resulting in further complications during the follow-up and healing processes [11,27]. In addition, the other limitations of PTC include difficult fiber placement, non-uniform injury, and mismatching with the target tumor shape, leading to unexpected side-effects after the PTC treatment. Moreover, the safety and efficacy of thermal treatments, including MWA, RFA, and PTC, [5] also rely on the number, size, and location of the tumors and on the techniques employed [28,29]. MWA, which uses electromagnetic (EM) signals (∼ 915 MHz-2.45 GHz) to heat tissues, could offer several benefits, such as the treatment of larger necrosis volumes (i.e., larger than 2-3 cm) and a lower sensitivity to heat-sink effects in comparison with RFA [30]. Heat-sink effects usually occur when tumors are located proximal to large blood vessels, leading to incomplete treatment. Furthermore, MWA might encounter some critical issues during the treatment of small-sized tumors (diameter < 2 cm) because of cumbersome antennas and long exposure times [5]. Conversely, RFA, which uses a high-frequency heating source (∼375-500 KHz), may effectively treat small-sized tumors (< 3 cm) [30]. However, RFA usually employs a standard catheter with a 4 mm tip for clinical applications. This could result in difficulties in insertion and operation for treating tumors that are ∼ 8.5 mm diameter, as in this case, the catheter tip diameter is a half of the tumor diameter. Mellal et al. agreed that RFA remains a controversial technique, as it suffers from a lack of accuracy and precision [31]. In contrast, Tombesi et al. suggested that PTC, which uses a fine optical fiber (core diameter = 300-600 µm) to achieve a low light absorption into the tissue via water, precisely induces tissue coagulation and treats small-sized tumors (< 2 cm) in a rapid manner [5,30]. Therefore, PTC is considered a good choice for PTMC treatment when the tumor diameter is smaller than 8.5 mm [29,32]. The current experimental study deployed a fine optical diffuser to obtain the moderate coagulation volume of 0.42 ± 0.08 cm 3 after combined modulation for 100 s; this coagulation volume was larger than the tumor volume (∼0.32 cm 3 ). According to the Arrhenius equation, thermal injuries are primarily a function of the irradiation time and the local temperature in the tissue. However, the current research maintained the temperature at 70 ℃ in the combined mode; therefore, the degree of the thermal injuries in the tissue was dependent on the irradiation time only (Fig. 5). The thermal damage induced can be increased by prolonging the irradiation time in a controllable and predictable manner to treat PTMC ranging from 3 to 10 mm in diameter. Therefore, the combined modulation technique, which maintains a constant temperature, may offer other potential benefits, such as the treatment of lesions in high-risk (e.g., near the gallbladder or main biliary duct) and hard-to-reach locations (e.g., the caudate lobe and hepatic dome) without carbonization issues. Figure 1(D) depicts the ultrasound-assisted percutaneous process of insertion of the customized device (including the optical diffuser and thermocouple) into the tissue and the feasible management of the thermal damage under 1064 nm laser light irradiation. There was a visibly safe route for the insertion of the optical diffuser and thermocouple and for the determination of their advancement in the region of interest (i.e., the target lesion) [33]. No significant puncture-related difficulties occurred. Without ultrasound imaging, the treatment was considered a blind procedure in this study. However, with ultrasound assistance, treatment images were acquired in the three-dimensional perspective and displayed on the computer screen for further processing. In comparison with the MRI and OCT techniques, ultrasound imaging has a higher portability and is more cost-effective [28,34]. Under the ultrasound guidance, the gradual coagulation of the tissue could be observed with the temperature increase during 1064 nm laser light irradiation. The tissue structure was altered, and echoes were displayed as bright white areas (Fig. 4(C)) [35]. Furthermore, hypoechoic regions were observed in the CW mode owing to water vaporization (temperature ∼ 100 ℃), whereas homogeneous echo areas occurred in the combined mode owing to the relatively low temperature (∼ 70 ℃). Nevertheless, hypoechoic rims were hardly detectable in the lesion images. The acoustic shadowing associated with early ultrasound images of acute thermal damages makes a full "rim" difficult to detect and evaluate [36]. For clinical applications, the sound attenuation caused by air gaps and the deep location of the organs in the human body should be considered by using a gel as the interface between the transducer and the skin [37]. Figure 3 shows a slight difference in temperature between the simulated and experimental results. This was caused because in this simulation, the tissues were assumed to be either completely native or completely coagulated; thus, the optical properties of tissues in the transition states were not considered because of lack of data. Regarding the combined mode, the temperature (70 ℃) in the simulation remained unchanged after 20 s of treatment. However, the temperature increased slightly in the experiment. The margin of the coagulation regions was rather unclear in the experiment, resulting in measurement overestimation. Thus, an approximately 10% difference in thermal damage volumes was commonly observed between the simulation and experimental results (Figs. 4 and 5) [11].
Although the present study revealed some of the advantages of the combined modulation, several limitations must be overcome to validate the feasibility of this modulation technique into the clinical applications. A new optical diffuser should be designed and machined to emit a flat-top light profile to be close to the simulated assumptions. In this study, liver tissue was used instead of thyroid tissue owing to the ease procurement of liver tissue and the lack of PTMC data. Interestingly, absorption coefficients of native thyroid and native liver tissues at 1064 nm are comparable (0.04 mm −1 vs. 0.05 mm −1 ) [19]. In addition, Sekar et al. confirmed that the optical scattering in thyroid at 1064 nm predominates over absorption (0.75 mm −1 vs. 0.04 mm −1 ) [38], which is similar to the liver tissue. However, the slightly lower optical scattering coefficient of the thyroid tissue (0.75 mm −1 for thyroid vs. 8 mm −1 for liver) is expected to accompany a narrower temperature distribution with a higher maximum temperature on account of less lateral photon distribution [38]. Furthermore, as the current experiments were conducted under ex vivo conditions, blood perfusion effect were assumed to be insignificant for the numerical simulation due to lack of pulmonary circulation [22]. However, under in vivo and clinical situations, many blood vessels surrounding the papillary thyroid carcinoma lesion can play a main role in cooling the induced temperature via convective heat transfer, determining the extent of the irreversible coagulation volume in the target tissue. In fact, additional simulation was conducted to verify the potential effect of the blood vessels on the interstitial thermal responses by incorporating the blood perfusion term into the Eq. (2). Blood parameters including the density (ρ b = 1050 kg/m 3 ), the specific heat (c b = 3640 J/kg·°C), and perfusion rate per unit volume (ω b = 0.0182 m 3 /m 3 ·s) were derived from previous reports [21,39]. The simulation results confirmed that for both irradiation modes, the blood perfusion effect reduced the maximum tissue temperature and the coagulation volume by around 10% and 30%, respectively, compared to no blood perfusion. The maximum tissue temperatures were 80 ℃ and 66 ℃ while the coagulation volumes were 0.32 and 0.27 cm 3 for CW and combined modes, respectively. These findings are comparable with the previous reports [21,39,40]. The differences in the temperature and the coagulation volume thus implicate that the real tumor may need higher irradiation energy to achieve the predetermined temperature (70 ℃) and the predetermined coagulation volume (0.32-0.52 cm 3 ) during the treatment in spite of a good correlation between the current simulation and ex vivo experimental results. For clinical translations, further in vivo studies should be performed in tumor-bearing animal models to validate the current findings along with the effect of multiple or larger blood vessels on interstitial temperature developments, the extent of thermal injuries, and wound healing after photothermal treatment. In addition, owing to the advantages of large depth (ultrasonic imaging) and high resolution (optical imaging), an ultrasonic-integrated photoacoustic imaging can be employed to acquire the blood perfusion rate prior to the treatment and to incorporate the acquired data into the simulation of the accurate treatment planning [41]. Moreover, photoacoustic imaging can provide the real-time feedback on optical property changes (absorption and scattering coefficients) as photoacoustic signals are dependent on the tissue condition [42]. Hester et al. reported that an ultrasound imaging system integrated with a pulsed laser acquired PA images of the target tissue lesion in real-time after acoustic wave generation through thermoelastic expansion upon light absorption [43]. Thus, instead of invasive thermocouple measurements, the proposed method can be employed concurrently with the photoacoustic imaging modality to externally measure the tissue parameters from photoacoustic signals and eventually to monitor the interstitial temperature development in a non-invasive manner [12].

Conclusions
Numerical and experimental validations of the thermal responses of irradiated tissue demonstrate the feasibility of the combined energy modulation for the potential treatment of PTC. Maintenance of the interstitial tissue temperature entailed a linear correlation between tissue coagulation volumes and irradiation time once the predefined temperature was set at 70 ℃. The predetermined volume of thermal coagulation can be optimized for the treatment of tumors with a diameter of ∼ 8.5 mm (volume ∼ 0.32 cm 3 ) without carbonization. The proposed energy modulation method is a feasible modality for performing predicable thermal coagulation to entirely treat the predetermined volume of the tumor and must be further validated for the treatment of PTMCs.

Funding
Ministry of Health and Welfare (HI16C1017).