Assessing temperature increase during photodynamic therapy: a simulation model

Photodynamic therapy (PDT) is a non-invasive procedure mainly used for treatment of malignancies. It is based on cytotoxic products generation after the excitation of absorbed photosensitizing drugs by non-ionizing radiation. The wavelength chosen is usually in the 620 - 700 nm region, resulting in limited tissue penetration depth. The relatively low irradiance values and the treatment time that, in most cases, lasts for a few minutes result in limited photothermal effects, which are usually overlooked. Therefore, in this study we computationally assess the temperature distribution during PDT in a cancer bearing mouse model. A user-friendly application is created that could be used in the treatment planning step of PDT. It receives as input various parameters, such as laser power, beam radius, irradiation time, body temperature and returns the maximum tissue temperature and irradiance values. Furthermore, it offers visualization of the generated effects through the spatial distribution of temperature, irradiance, acute necrosis and damaged tissue percentage that are presented in the form of interactive 3D plots. The conducted simulations reveal that 43 °C, which are in the hyperthermia range, are difficult to be excessed in PDT clinical practice, although topical thermal effects are observed.


Introduction
Photodynamic therapy (PDT) is a clinical procedure used for cancer treatment, as well as for skin and eye disorders (e.g. acne [1] and macular degeneration [2]). It is based on three, simultaneously acting fundamental agents, i.e. photosensitizing drug (PS), molecular oxygen ( 3 O2) and light. The PS molecules are excited by irradiation with light of the appropriate wavelength (λ) and the cytotoxic products generated kill the cancer cells [3]. Among the main advantages of PDT one could find the non-invasive nature of the procedure, as the excitation light irradiates the tissue surface and the PS is topically or intravenously administered. Moreover, the non-ionizing part of the electromagnetic spectrum used for the illumination allows the treatment repetition and the sparing of the adjacent tissue. The latter is further supported by the high localization of the PS molecules in the cancer cells and the selectivity of drug activation only under illumination by light of certain wavelength.
The procedure itself is rather simple, compared to other therapies, since there are no strict, globally accepted guidelines and protocols for the physicians to follow, especially regarding personalized treatment planning. In many cases the protocols followed are the ones proposed by the PS suppliers and the physicians seem to overlook dosimetry guidelines as the Report No. 88 of American Association of  [4,5]. Furthermore, to the best of our knowledge, medical physicists do not always participate in the treatment process or in the pre-treatment simulation and dosimetry.
The current trend in PS development is a swift towards absorption coefficients (μa) in red and infrared (IR) wavelengths, due to their increased penetration depth [6]. For the most PSs the irradiation wavelength is in the 630 to 700 nm region, resulting in superficial effects; while the beam is absorbed mainly by the drug, melanin and haemoglobin molecules and practically not at all by water [6,7]. The relatively low irradiance (E, expressed in mW/cm 2 ) values and the treatment times used in clinical practice, which in most cases last for a few minutes, result in limited photothermal effects that are usually not taken into account.
Hence, in this computational study the temperature distribution during PDT in a cancer bearing mouse is assessed. A user-friendly application is created that receives as input various parameters, such as laser power, irradiation time and body temperature, and returns the maximum tissue temperature and irradiance values at user specified time points. Moreover, the spatial distribution of temperature, irradiance, acute necrosis and damaged tissue percentage are presented in the form of interactive 3D plots. Therefore, through this application the physician and/or the medical physicist can adjust the treatment parameters, in order to predict, avoid or induce hyperthermic effects.

Materials and methods
A tumor bearing mouse is created using the COMSOL modeling software (COMSOL Multiphysics®, v.5.5, Stockholm, Sweden), which consists of muscle tissue, an elementary skeleton and a 220-μm skin layer. A breast cancer tumor, corresponding to the 4T1 cancer cell line, is designed subcutaneously in the left front flank of the mouse (see figure 1 (a)). The mouse is subjected to PDT with Metvix, a PS that is metabolized to intracellular porphyrins, including Protoporphyrin IX (PpIX), which is highly photosensitive and mainly accumulates in cancer cells [8]. The thermal effects of the PpIX molecules are assumed to be negligible due to their low μa value at the wavelength used (625 nm), compared to the one of the cancer tissue [9]. The corresponding values are based on the available literature [10][11][12] and are presented in table 1. The light beam is of gaussian profile and is produced by a 625 nm laser. The COMSOL Bioheat Transfer module is used for the photothermal study and the parameters values used are presented in table 1. The simulated tumor and irradiation setup correspond to the ones used in a preliminary study of ours [10]. The software application builder is used to create a program with userfriendly interface (UI), executable for both COMSOL and non-COMSOL users. In the relevant fields (see figure 1) the following parameters are defined: tumor position (in Cartesian coordinates) and radius (rt), light source power (P) and radius (rs), treatment time (t), body temperature (Tb) and simulation accuracy via mesh size. Furthermore, the user is able to set irradiation time and temperature threshold values (tth and Tth, respectively) for cancer cell death via apoptotic procedures as well as for acute necrosis [13]. The application calculates the maximum temperature, irradiance values and damage percent values; while spatial distribution of temperature, irradiance, damaged tissue percentage and acute necrosis are presented in the form of interactive 3D plots.  (c), although the maximum and minimum values are displayed in yellow font (in the online, colored version), they are not fully visible in this snapshot due to the selected, for clarity reasons, mouse orientation. As it can be seen from the graphs, the user may also choose to project the desired treatment time point (in this case 0 s in (b) and 500 s in the rest).
In order to assess the effect that PDT has on the mouse temperature fifteen (15) case scenarios are studied. Each one of them is based on a fixed irradiation time, while irradiance values range from 50 to 250 mW/cm 2 , with a 50 mW/cm 2 step. For the calculations presented the mesh size is set to "fine".  Figures 2 (a) and (c) are indicative of the results returned by the application, regarding temperature and damage fraction spatial distribution. For demonstration purposes, t and E values are higher than the ones used in clinical practice. Specifically, t is set at 45 min and E at 350 mW/cm 2 . The maximum recorded temperature is 43.9 o C, while less than half of the tumor volume has been subjected to irreversible damage.

Results and concluding remarks
The case scenarios studied reveal that, for a fixed irradiance, there is no significant difference regarding the maximum temperature values observed at the end of the treatment, regardless of the irradiation time chosen (see figure 2 (b)). Moreover, the temperature values recorded, consistent with literature [14], are not sufficient to induce hyperthermic effects but high enough to affect PpIX photobleaching rate and PDT efficacy [15,16]. The temperature increase between the 5 and 15 min scenarios for each irradiance value ranges from 0.30 to 1.25 % (from 50 to 250 mW/cm 2 ), while the corresponding increase between the 15 and 30 min scenarios ranges from 0.06 to 0.20 %. In order to assess the temperature behavior during PDT treatment, a 5 min -200 mW/cm 2 scenario is presented in figure 2 (d). The temperature values present a logarithmic rise, denoting the significance of the first treatment seconds.