Abstract

Photodynamic Therapy (PDT) modeling allows the prediction of the treatment results depending on the lesion properties, the photosensitizer distribution, or the optical source characteristics. We employ a predictive PDT model and apply it to different skin tumors. It takes into account optical radiation distribution, a nonhomogeneous topical photosensitizer spatial temporal distribution, and the time-dependent photochemical interaction. The predicted singlet oxygen molecular concentrations with varying optical irradiance are compared and could be directly related with the necrosis area. The results show a strong dependence on the particular lesion. This suggests the need to design optimal PDT treatment protocols adapted to the specific patient and lesion.

1. Introduction

Photodynamic therapy is a cancer treatment modality based on various types of photochemical interactions involving an extrinsic chemical agent known as photosensitizer, oxygen, and light [1]. The biological effects of localized oxidative damage caused by PDT provoke tumor selective destruction. PDT provides advantages over traditional techniques such as radiation therapy or chemotherapy, like noninvasivity, selectivity and the lack of severe collateral effects. PDT is currently applied to several tumoral tissues and it is particularly appropriate in dermatology, due to the noncontact character and the good cosmetic results [2]. Nonmelanoma skin cancer is one of the most frequent malignant neoplasia in humans, which includes basal cell carcinomas and squamous cells carcinomas or precursor cancer diseases such as actinic keratosis. These skin disorders have grown during the last years due to an increase of solar exposition of the population, the increase of the longevity, or the lack of ozone in the atmosphere. Conventional techniques against these pathologies like surgery, electrocoagulation, cryosurgery, radiotherapy, or even pharmacological solutions present some limitations that are successfully solved by PDT [3].

PDT in clinical praxis employs usually a fixed protocol that is independent of the lesion and of the patient. As a consequence some pathologies present recurrence some time after the application of the treatment. Clinical praxis in the Dermatology Department of Marqués de Valdecilla University Hospital shows that the treatment response varies with the type of pathology treated [4]. Regarding the response of the patients to the treatment, quite good results have been obtained without relevant side effects. These effects are mostly heat and slight pain and they appear in the great majority of the cases.

An optimal adjustment of the treatment is convenient, depending on the light source and the photosensitizer parameters, and also on the type of skin disease and the patient. In recent years there have been different works regarding the improvement of light delivery, the strategy for the photosensitizer delivery, the development of new photosensitizers, or a better oxygenation of the treated tissue [5, 6]. The treatment adjustment requires a PDT predictive modeling approach. There are a limited number of works related to the development of models that characterize the photodynamic process in order to optimize the treatment [7]. In these studies the focus is on the development of photochemical [8, 9] and optical propagation models [10]. Most PDT models that take into account the whole photodynamic process usually consider a systemic application of the photosensitizer [11], with a uniform spatial distribution in the target tissue. However, there are several studies that show that in the case of a topical photosensitizer the distribution of the drug through the skin results in an inhomogeneous spatial pattern [12, 13]. Significant advances have been made in the study of the photosensitizer photobleaching and oxygen diffusion and perfusion in the tissue [14, 15]. Liu et al. [16] have recently studied singlet oxygen direct dosimetry and photosensitizer fluorescence photobleaching as dosimetry tools during 5-aminolevulinic acid (ALA) induced Protophorphyrin IX (PpIX). A one-dimensional model was proposed to simulate the dynamic process of ALA-PDT of normal human skin in order to investigate the time-resolved evolution of PpIX, ground state oxygen, and singlet oxygen distributions. Although this work presents one of the most complex models used up to now, the tissue model employed is quite distant from a realistic clinical application due to the fact that it is applied to normal human skin simplified as a three-layer semi-infinite tissue. The influence of the type of pathology in the result of PDT treatment has been explored in several clinical studies, in which patients exposed to PDT show different treatment responses under the same dosimetry protocol. However, as far as we know, treatment outcome prediction by means of a PDT photochemical model for different skin tumors has not been realized.

The present work is devoted to the application of a complex PDT model to different types of human skin tumors. The prediction of the PDT treatment response depending on the particular lesion would allow the design of optimal PDT treatment protocols adapted to the specific patient and lesion. The PDT model used takes into account optical radiation distribution by means of a Monte Carlo approach, a nonhomogeneous topical photosensitizer spatial temporal distribution, the photoactive compound synthesis and the time-dependent photochemical interaction. The results of the different elements’ concentrations in the photochemical model, specially singlet oxygen production, could be directly related with the necrosis area. The model is applied to skin tumors, such as nodular and infiltrative basal cell carcinomas (BCCs), or squamous cells carcinomas (SCCs). The predicted singlet oxygen molecular concentrations of the different tumoral skin tissues are compared with varying optical irradiance.

The manuscript is organized as follows. Section 2 describes the PDT model employed in detail, including the optical propagation, photosensitizer distribution, and photochemical analysis. In Section 3 the model is applied to several skin tumors, the treatment parameters are varied, and the predictive results are obtained. These results are compared and discussed. The conclusions of this study are contained in Section 4.

2. Photodynamic Therapy Model

The predictive PDT model takes into account optical radiation distribution, a nonhomogeneous topical photosensitizer spatial temporal distribution, and the time-dependent photochemical interaction. Each part is described in the following subsections.

2.1. Optical Radiation Distribution

There are different optical models to be used in biological tissues. Modeling a biological tissue implies dealing with an heterogeneous medium, which does not allow a analytic exact approach of the radiation pattern with Maxwell equations. For the problem we are dealing with, the distribution of light in a three-dimensional tissue must be obtained. This objective is reached by means of the Radiation Transport Theory (RTT) [1]. The model assumes that the scattering events are sufficiently numerous as to the light to be considered incoherent, in such a way that polarization or interference effects can be neglected. As a consequence, the basic parameter of light is the specific intensity , that is, the light power per unit area per unit solid angle. The radiation is expected to be at point , and to follow the direction . The scattering events are treated according to the scattering phase function . Optical radiation comes from direction and is redirected to . The basic idea in order to write the differential radiation transport equation is that radiation from a particle attenuates due to absorption and scattering and also gains power because another particle can scatter light in the direction of the particle of interest. If there are no sources inside the tissue and a steady-state situation, this can be written as Numerical analysis has been widely applied to a great variety of problems governed by differential equations. In the particular topic of the radiation transport equation, the Monte Carlo method has demonstrated its applicability and accuracy, compared with exact solutions. Perhaps the most used implementation of the Monte Carlo method applied to the RTT model is the one by Wang et al. [10, 17]. This implementation of the Monte Carlo model is multilayered, with their borders always perpendicular to the laser beam. This is very useful because tissues usually can be divided in different strata. For the appropriate definition of the model, the characteristics and dimensions of each layer are required. The optical parameters needed are the index of refraction , the absorption coefficient , the scattering coefficient , and the anisotropy of scattering .

2.2. Topical Photosensitizer Distribution

Methyl Aminolevulinate (MAL) is a prodrug that is endogenously metabolized to the photoactive element Protophorphyrin IX (PpIX). After its application on the skin surface, a diffusion process through the different skin layers occurs and it is converted to PpIX. Several studies have evaluated the PpIX content of the skin using fluorescence techniques [12] and the influence of the stratum corneum as the main barrier to the diffusion of the photosensitizer to deeper layers of skin [13].

The inhomogeneous distribution of a topical photosensitizer precursor through the skin and the photosensitizer endogenously produced play an important role to determine the concentration of photoactive substance to be accumulated during the incubation period. For this reason we used Fick’s law to characterize the inhomogeneous photosensitizer precursor distribution and to calculate the concentration reached at each point of the tissue during the incubation period. According to Fick’s law, if there are differences of concentration of a substance, its molecules move from higher to lower concentration regions, and so the flow of substance goes in the opposite direction of the concentration gradient: where is the flux vector indicating the direction and magnitude of substance, is the diffusion coefficient, is the prodrug concentration, and is the depth in the tissue.

The distribution of the photosensitizer in the skin is limited by several factors, including the stratum corneum which acts as a diffusion barrier and is characterized by the permeability, , the diffusion coefficient through the epidermis and dermis, , the relaxation time of the precursor as a consequence of the generation of the photosensitizer and other processes (lymphatic flow and blood perfusion), , and the conversion rate of photosensitizer precursor in its photoactive compound [18]. The temporal evolution of the photosensitizer precursor concentration for each depth in the tissue sample can be calculated as follows: where is the concentration of photosensitizer precursor in the skin surface at and is the distance from the corneal layer located at . Once the concentration of MAL is known at each point, the accumulated concentration of active substance in the tissue during an incubation period of 3 hours is calculated. This lets us know the amount of photosensitive substance at every point of the cancerous tissue at the beginning of the radiation interval. It is assumed that the PpIX relaxation time is fast compared to the MAL diffusion time,, and therefore the concentration of PpIX is proportional with the instantaneous value of MAL concentration. This value can be calculated by means of equation (4), where is the yield of the conversion process and the relaxation time of the photosensitizer precursor due to the generation of PpIX:

2.3. Photochemical Interaction

Interaction of light with a photosensitizer at an appropriate wavelength produces an excited triplet state photosensitizer that interacts with ground state oxygen via two types of reactions, known as Type I and Type II. The Type II reaction is believed to be predominant and responsible for singlet oxygen production, which is considered as the cytotoxic element in charge of killing carcinogenic cells.

The photochemical reactions are characterized by means of a photochemical model [8, 9], which takes into account the transitions between states of the particles involved, such as the photosensitizer and oxygen. The solutions of the stiff differential equations system employed, (5), are obtained by means of a differential equation solver within the Matlab platform. In order to obtain coherent results, we had to adjust relative and absolute error tolerances. These solutions provide the temporal evolution of the different molecular concentrations of the compounds involved in a Type II reaction everywhere in the tissue sample: In these equations, is the concentration of the photosensitizer in ground state, is the concentration of the photosensitizer in singlet excited state; is the concentration of photosensitizer in triplet excited state; is the concentration of oxygen in ground state; is the concentration of singlet oxygen; is the concentration of singlet oxygen receptors; is the scavengers concentration; is the relaxation time from state to ; is the relaxation time from state to ; is the relaxation time from state to ; is the quantum yield of the transition from state to ; is the quantum yield of the transition from to ; is the quantum yield of transition to ; is the quantum yield of transition to ; is the efficiency factor for energy transfer from to ; stands for the biomolecular photobleaching rate; is the biomolecular cytotoxicity rate; is the rate of reaction of with various oxygen scavengers; is light speed in tissue; is the photon density present at a point; is the absorption cross-section of molecules; is the rate of oxygen diffusion and perfusion; is the cell damage repair rate.

3. Results and Discussion

The optical radiation distribution was obtained by means of the Monte Carlo method presented in the previous section. A cylindrical laser beam with a radius of 0.3 cm perpendicular to the tissue sample was used to deliver different irradiances (100, 150, and 200 mW/cm2). The radiation time was 10 minutes. The tissue model is composed of two layers, the upper one was corresponding to the tumor zone and the other below was considered muscle. The optical parameters corresponding to the three skin pathologies are listed in Table 1 [19, 20].

Figure 1 displays the optical absorption (J/cm3) in three types of human skin tumors (nodular BCC, infiltrative BCC, and SSC) for the illumination conditions previously showed. As it can be observed optical absorption varies depending on the type of pathology and anatomical position in the tissue sample. In deeper tumor positions optical absorption is lower than in superficial zones and so the photon density that will excite the ground state photosensitizer diminishes with depth. Larger power densities induce larger absorption in all the pathologies.

The tumor optical properties affect final optical distribution and make optical energy vary depending on the biological media. The total energy in the whole tissue model was calculated for the three types of tumor. The results are listed in Table 2, where it can be observed that an infiltrative BCC accumulates less energy than the other two types of skin pathologies for all the irradiances applied to the tumor surface. The SCC tumor presents larger energy accumulation than the others.

The spatial and temporal diffusion of a topical photosensitizer precursor in the tissue sample and the photosensitizer endogenously produced were obtained during an incubation period of three hours. The corneal layer reduces the permeability of the skin, so that its value can be adjusted to characterize different skin conditions. In this case it is damaged or reduced by the skin lesion corneal layer. As a consequence, the value of permeability, K = 10−6 m/s, was chosen greater than that in the case of healthy skin with an intact cornea [18]. The diffusion coefficient through the epidermis and dermis was 0.69·10−10 m2/s, the relaxation time of the prodrug 24 hours, the yield of the conversion process 0.5, the relaxation time of the photosensitizer 84 ms, and the relaxation time of the precursor due to the generation of the photosensitizer 25 hours [18, 21]. The clinical Metvix-PDT protocol only specifies that Metvix is applied in a 1 mm thick layer on the affected area covering an extra 5 mm of healthy skin around the damaged area and has an incubation period of 3 hours before radiation. Therefore in order to estimate the concentration of the PpIX precursor applied on the skin tumor as close as possible to a real case, we have taken into account the area of the pathology and MAL density per gram of Metvix cream. Assuming a circular lesion of radius 1 cm, a MAL density of 160 mg per gram of Metvix cream, and the molecular mass of MAL, we calculated a MAL concentration on the pathology surface of 4.5031·1020 cm−3.

The results are presented in Figure 2. As it can be observed, the barrier blocks the diffusion to deeper layers during the first part of the photosensitizer incubation process. This phenomenon causes an uneven distribution of the photosensitizer in the skin during the period of incubation. The initial photosensitizer concentration when the optical irradiation starts will not be the same at every point of the tissue. Therefore there will be differences in the results of the photodynamic procedure. This is an important aspect to consider when dealing with a topical photosensitizer precursor, since in the case of a systemic photosensitizer a homogeneous distribution is generally considered. The photosensitizer concentration three hours after its administration was used as an initial condition to the photochemical model, which calculates the temporal evolution of the molecular components involved in the photochemical reactions by means of the stiff differential equations system previously shown.

Typical clinical conditions for the photosensitizer Metvix for PDT were considered in this study. Literature related to well-known photosensitizers of the porphyrins family was employed to assign the parameters values when they were not available for this photosensitizer [6]. The absorption cross-section of PpIX molecules (= 0.37·10−15 cm2) at the treatment wavelength was derived from a study of the cellular photosensitizing properties of PpIX carried out in a transformed murine keratinocyte cell line [22]. Due to the fact that MAL is a derivative of 5-aminolevulinic acid, the relaxation time from singlet excited state to ground state to was set to 7.4 ns as reported earlier from fluorescence measurements in cells incubated with 5-aminolevulinic acid induced PpIX [23]. The triplet state lifetime in vivo in skin (τ3) was set to 26 μs and the relaxation time of to to 0.04 μs [24]. Quantum yield transitions between different energetic states were adopted to be similar to those previously considered for the photosensitizer Photofrin = 0.2, = 0.3, = 0.3, and = 0.8 as well as biomolecular photobleaching, cytotoxicity, and scavenging rates that were set to 2·10−10 cm3/s−1, 2·10−9 cm3/s−1, and 1·10−9 cm3/s−1, respectively [9]. At the beginning of the optical irradiation period there are not molecules of photosensitizer in excited state; so the initial concentrations of and are 0 cm−3. In the same way singlet oxygen molecules have not yet been produced and their initial concentration is 0 cm−3. The tissue sample was considered homogeneously oxygenated with an initial concentration of cellular oxygen of 5·1017 cm−3 and diffusion and perfusion rates of 1·1012 cm−3·s−1 [9]. Initial concentration of intracellular molecular singlet oxygen receptors was 5·1017 cm−3 scavenger concentration was 1·103 cm−3, and cell damage repair rate was 2.6·1012 cm−3·s−1 [9].

The temporal evolution of the photosensitizer at different energetic states, oxygen, and singlet oxygen molecular concentrations was obtained by means of the photochemical model for the three different skin tumors and varying optical irradiance. The concentrations of ground state molecules of photosensitizer, oxygen in ground state, and receptors vary slowly in response to the activation light. While applying optical radiation, the concentrations of the triplet state molecules of photosensitizer and singlet oxygen in all the tissue samples begin to increase slowly until they reach a maximum value.

Different maximum singlet oxygen concentrations were reached for different types of skin pathologies and depths; so different treatment outcomes could be expected. Figure 3 shows the maximum singlet oxygen concentration reached at every point of the tissue sample, depending on the target anatomical position. The spatial variation is presented according to the depth from the tissue surface and the distance to the center of the laser beam.

It can be observed for all the pathologies that an irradiance increment induces a larger amount of singlet oxygen in the target tissue. However, this amount is not the same in all of them. Figure 4 shows that the pathology which accumulates a larger singlet oxygen concentration depends on the observation depth from the tumor surface. At the maximum tumor depth of 3 mm singlet oxygen produced is larger in the nodular BCC (2.746·104, 4.204·104, and 5.651·104 cm−3 for 100, 150, and 200 mW/cm2, resp.) than that in the other two pathologies independently of the applied irradiance (see Figures 4(a)–4(c)), whereas the pathology which accumulates the lowest maximum singlet oxygen concentration is infiltrative BCC (1.731·104, 2.654·104, and 3.531·104 cm−3 for 100, 150, and 200 mW/cm2, resp.). However when the observation depth is more superficial, 0.1 mm from the tumor surface (Figures 4(d)–4(f)), infiltrative BCC presents the largest maximum singlet oxygen concentration (7.229·108, 1.091·109, and 1.428·109 cm−3 for 100, 150, and 200 mW/cm2, resp.) and SCC the lowest (5.673·108, 8.459·108 and 1.118·109 cm−3 for 100, 150 and 200 mW/cm2 resp.). The variation in singlet oxygen production at different depths when comparing the lesions is caused essentially by the particular optical radiation distribution in each case. As some pathologies accumulate preferentially optical radiation at lower depths, the singlet oxygen production is lower deeper in the tissue, where energy absorption is smaller.

In Figures 4(d)–4(f) a clear influence of the optical beam radius in the maximum singlet oxygen concentration accumulated in the target tissue can be observed. The most superficial tissue zone shows a significant singlet oxygen concentration decrease as the measurement point gets far from the laser beam impact center. Differences in the maximum amount of singlet oxygen molecules generated in the volume occupied by the tumor can also be observed. They vary depending on both irradiance and type of tumor, as it is shown in Table 3. For a fixed irradiance, the amount of singlet oxygen molecules in the volume occupied by the tumor varies depending on the type of pathology. It is maximum in the infiltrative BCC tumor (1.1995·107, 1.8165·107 and 2.4246·107 molecules of for 100, 150 and 200 mW/cm2 resp.) and minimum in the SCC (9.9338·106, 1.5072·107 and 2.0199·107 molecules of for 100, 150 and 200 mW/cm2 resp.).

When comparing Tables 2 and 3, it can be seen that the total optical energy absorbed by the tumor is not proportionally related with the total singlet oxygen molecules. For instance, SCC presents a higher value of total optical energy, but the lowest number of total singlet oxygen molecules. This effect is related with the nonlinear character of the photochemical process. As optical energy is bigger at a specific point, singlet oxygen production efficiency decreases. Optical energy is preferentially accumulated in the superficial layers of the SCC (see Figure 1), and so deeper layers present lower-energy absorption. This makes the whole singlet oxygen production efficiency decrease.

All the results previously shown indicate that the treatment effect, which is closely related to the amount of singlet oxygen produced, varies depending on the type of tumor we are dealing with. Therefore the type of pathology should be taken into account to adjust the irradiation source parameters, such as the optical laser beam radius and the irradiance applied to the pathology surface to obtain the maximum cancerous tissue necrosis. This approach would serve to predict the damage induced in the types of pathologies studied if a precise singlet oxygen threshold was established as the amount necessary to induce cell toxicity in each case.

4. Conclusions

Rigid protocols employed in PDT, which are independent of the particular lesion and the patient, were shown to compromise the treatment efficiency according to clinical studies. These studies identify different treatment responses and the recurrence of the skin disorders as the main inconvenient of PDT. In order to remove this problem the adjustment of all the parameters involved in the photochemical reactions is proposed, as a function of the type of skin disease. This adjustment was analyzed by means of a complex predictive PDT model.

The PDT model employed takes into account optical radiation distribution, a nonhomogeneous topical photosensitizer spatial-temporal distribution, and the time-dependent photochemical interaction. The results obtained in the present work show a clear influence of the type of skin tumor under treatment in the effect of photochemical reactions induced in PDT, specifically in the amount of reactive oxygen produced. Under the same illumination conditions and assuming the same photosensitizer distribution for the three types of skin tumors modeled (nodular BCC, SSC, and infiltrative BCC), differences in the maximum amount of singlet oxygen molecules accumulated in the tumoral region were observed, depending on the type of pathology. These differences induce different treatment responses, because the amount of reactive oxygen molecules is related to cell destruction processes and necrosis. Therefore it was demonstrated that treatment outcome also depends on the type of pathology. As a consequence actual clinical protocols should be adjusted taking into account factors related to the type and morphology of the tumor we are dealing with.

Acknowledgments

The authors thank María López-Escobar (Dermatology Department, Marqués de Valdecilla University Hospital) for providing the results of the clinical application. This work has been partially supported by the Leonardo Torres Quevedo Foundation and San Cándido Foundation.