2.1.

Clinical Topical Photodynamic Therapy Treatment

PDT treatments were carried out at Ninewells Hospital and Medical School in Dundee, Scotland, United Kingdom, using the photosensitizer prodrug, 5-ALA. Formal consent was obtained from the patients before the study was undertaken. T-PDT was preformed using either of two light-emitting diode (LED) light sources, Aktilite CL16 and CL128 (PhotoCure ASA, Oslo, Norway). These commercial devices have identical spectra, where the peak irradiance occurs at a wavelength of 632 nm with a full width at half maximum of 19 nm. Both light sources were calibrated and traceable to national measurement standards.¹⁷ Typical treatment times and skin surface irradiances are in the region of 15 min and 80 mW/cm², respectively. At our centre, a total treatment light dose (LD) of 75 J/cm² is administered during T-PDT. 5-ALA was topically applied to the lesional area. Following a 6-h incubation period; one skin lesion per patient was treated.

2.2.

Clinical Fluorescence Measurements from Superficial Basal Cell Carcinomas

PpIX is naturally occurring and present in small quantities in cells. Uptake and conversion of 5-ALA to PpIX in diseased tissue leads to increased levels of PpIX in lesions. When PpIX is excited with light at a wavelength of ∼405 nm, characteristic fluorescence is induced with a readily observable peak at ∼635 nm. A high fluorescence ratio of tumor to surrounding tissue allows for demarcation of the tumor. In our study, fluorescence spectra were acquired using a custom-built optical biopsy system (University of Glasgow, Glasgow, United Kingdom) formed by a fiber-coupled fluorescence spectroscopy system.¹⁸ The output from an LED emitting at 405 nm was coupled to a 400 μm core optical fiber probe used for both excitation and collection of fluorescent emission. PpIX fluorescence was induced in the skin, and spectra were recorded by placing the probe perpendicularly in contact with the tissue. From each individual lesion, eight fluorescence spectra were recorded before 5-ALA application, immediately before treatment, approximately halfway through treatment, and immediately after treatment. The mean spectrum was then calculated from each set of eight spectra for each lesion. Treatment was stopped briefly for ∼1 min to allow for the halfway measurements to be carried out. All measurements were performed in a darkened room to minimize artifacts from ambient room light.

2.3.

Description and Validation of the Monte Carlo Radiation Transfer Code

A fixed weight MCRT code was written and developed in FORTRAN 77 to simulate the scattering, absorption, fluorescence, and photobleaching processes present in T-PDT. It was possible to account for the generation and propagation of excitation photons, and of fluorescence photons emanating from a tumor embedded in normal skin tissue. We have also incorporated an algorithm in the code to determine the production of singlet oxygen in the tumor. Although other MC simulations of photon transport and fluorescence have been performed in multilayered tissues,^{19, 20} our 3-D MC code can calculate light distributions for a clinical T-PDT geometry where a tumor is embedded in normal skin tissue. Earlier studies have used this type of geometry to model the distribution of light within a tumor.^{21, 22} Our code, however, also incorporates different optical properties for normal and tumor tissue and the addition of a photosensitizer, namely, PpIX into the tumor. Furthermore, photobleaching has also been included. This enables the efficacy of T-PDT to be investigated.

Within our code, μ_a and μ_s are the tissue-absorption and -scattering coefficients, which can be varied depending on the depth within our phantom, and g is the scattering anisotropy factor. The Henyey–Greenstein phase function²³ is adopted to approximate scattering in tissue and is mathematically expressed in the form

Fig. 1

Three-dimensional (3-D) Monte Carlo model geometry, where x, y are mutually orthogonal axes in the skin surface and z represents depth within the skin tissue.

The Monte Carlo code was validated by comparing a range of simulations to the results generated by Keijzer,¹⁶ where light distributions in artery tissue were examined. For example, using a collimated incident-beam diameter of 1 mm, the fluence rates recorded from our simulations at the surface and at a depth of 0.5 and 1 mm into the tissue were 1.8, 0.55, and 0.11 W/cm², respectively. These results are comparable to those obtained by Keijzer ¹⁶

2.4.

Monte Carlo Model Assumptions

With clinical T-PDT in mind, it was necessary to make some assumptions in our model. We have assumed mismatched boundary conditions, and in a similar method to Farrell,¹⁰ that the tumor had a homogeneous distribution of PpIX. The fluence rate distribution decays with tissue depth. Photobleaching alters the spatial distribution of PpIX in the tissue. The top layers experience the highest fluence rate and therefore photobleach faster.¹⁰ We have assumed that the PpIX concentration has a negligible effect on the absorption coefficient of the tissue, therefore, resulting in a time-independent fluence-rate distribution.²⁷ PpIX photobleaching was assumed to follow a first-order exponential decay,²⁷ where the PpIX concentration, C(x, y, z, t) decreased at a rate proportional to the local fluence rate, Ψ(x, y, z), as is shown in

The dosimetric parameter, namely, the PD, is defined as the number of photons absorbed by the photosensitizer per unit volume of tumor tissue and is assumed to be proportional to the local yield of singlet oxygen production.¹⁰ A threshold photodynamic dose (PD_T), leads to tissue necrosis if the number of photons absorbed by the photosensitizer per unit volume of tumor tissue exceeds a certain value. A threshold value of approximately 8.6 × 10¹⁷ photons/cm³ absorbed by a porphyrin photosensitizer in tumors was used.²⁸ Again, following the work of Farrell,¹⁰ we assumed that the local yield of singlet oxygen production following absorption by the photosensitizer defined the PD and tissue necrosis was assumed to occur when PD > PD_T. The PD is directly related to the photosensitizer concentration and the light fluence rate, as is shown in

2.5.

Application of the Monte Carlo Model to Clinical Photodynamic Therapy

On the basis of our assumptions, we modeled a clinical situation, i.e., the treatment itself and the fluorescence measurements undertaken during treatment requiring two different sets of simulations to be undertaken. The former consisted of simulations where the singlet oxygen produced in the tumor was modeled as a function of depth and time, whereas the latter included modeling the PpIX fluorescence detected at particular time points in the treatment. In this way, our model could provide information about both the fluorescence information used as a diagnostic and the possibility of tumor cell kill. We sampled the fluorescence emission in the code from clinically measured 5-ALA-induced PpIX fluorescence spectral data. This approach is different from other models, where only excitation and fluorescence emission wavelengths have been used rather than full spectral data.^{27, 29} Figure 2 depicts the optical properties used in the Monte Carlo simulations, which were taken from data published in the literature.^{30, 31, 32} Figure 3 illustrates our absorption coefficient for PpIX as a function of wavelength. A stock solution of PpIX in dimethylsulfoxide was produced with a concentration of 45 μg/ml. The solution was then placed in a cuvette, where a spectrophotometer (Hitachi Spectrophotometer U-3010) was used to measure the optical absorbance and the absorption coefficient for PpIX was determined to be 0.06 cm⁻¹ at 630 nm, using an extinction coefficient for PpIX, ɛ_{630 nm} = 0.0014 μg/ml⁻¹ cm⁻¹. This is within the range of values published in the literature.¹⁰

Fig. 2

Spectral characteristics and model inputs of (a) normal skin and basal cell carcinoma scattering coefficients and (b) normal skin and basal cell carcinoma absorption coefficients. Data from Salomatina (Ref. 30).

Fig. 3

Absorption coefficient for PpIX as a function of wavelength.

For the PDT treatment simulations, we based the inputs for our model on the well-established treatment parameters used in the Scottish Photodynamic Therapy Centre, Ninewells Hospital & Medical School, Dundee. The entire top surface of the cube was irradiated uniformly, at the treatment wavelength of 632 nm. A surface irradiance, Ψ₀, of 82 mW/cm² was delivered over a simulated treatment time of ∼30 min, thus administrating a simulated total treatment LD of 150 J/cm². The photons propagated through the normal skin and tumor using the scattering and absorption coefficient values at this wavelength. After photon absorption in the tissue, the amount of energy deposited was calculated. The energy absorbed in the tumor was determined at each time step, corresponding to increments of 40 s. On the basis of our assumptions of the PD, it was possible to infer the 3-D distribution of photons and, hence, 3-D distribution of singlet oxygen production in the tumor as a function of LD.

For the PDT fluorescence simulations, PpIX fluorescence photons were induced at specific stages during the PDT treatment simulations. The entire top surface of the cube and cylinder was irradiated uniformly at the fluorescence excitation wavelength of 400 nm. Fluorescence photons that exited the top surface were tracked from their place of origin, and it was then possible to know exactly where they came from inside the tumor. After each PDT treatment simulation was carried out, the PDT fluorescence simulation followed immediately afterward. In other words, the two simulations were run simultaneously for given time steps (i.e., every 40 to ∼1800 s). In this manner, PpIX fluorescence was detected in the tumor as a function of depth and LD in the presence of photobleaching.

3.1.

In vivo 5-Aminolaevulinic Acid–Induced Protoporphyrin IX Fluorescence Spectra

Figure 4 illustrates the mean fluorescence of the eight spectra from each of the six patients before the application of the 5-ALA cream. Figure 5 presents PpIX fluorescence spectra noninvasively measured in six different patients as a function of time during treatment. Note that in both sets of figures, the overall mean of all six patients is highlighted by a thick solid line. As depicted by Fig. 4a, collected peak autofluorescence signals of ∼500 nm are highly variable between patients. Collected PpIX fluorescence signals of ∼635 nm can also vary between patients, as illustrated in Fig. 4b. Figure 5 illustrates that, after topical application of 5-ALA, there is a wide variability seen in the collected PpIX fluorescence signals among the six patients examined at each time point of treatment. Each baseline fluorescence spectrum measurement was subtracted from their corresponding fluorescence spectrum measurement at the three different time points of treatment (i.e., before, halfway through and after treatment).

Fig. 4

Baseline fluorescence spectra recorded before 5-ALA application for six patients; the mean is depicted in bold. (a) Spectral region showing tissue autofluorescence and (b) spectral region where PpIX fluorescence peak will feature.

Fig. 5

PpIX fluorescence spectra recorded for six patientsl the mean is depicted in bold. (a) immediately prior to treatment, (b) halfway through treatment, and (c) immediately after PDT treatment.

3.2.

Photobleaching Observed during Clinical Photodynamic Therapy

In vivo photobleaching was monitored through the use of the clinical patient sBCC fluorescence data collected and discussed above. Figure 6a depicts the mean PpIX fluorescence intensity from the six patients at the time points of treatment. As expected, this decreases from the outset to the end of treatment and illustrates how the peak PpIX fluorescence reduces as treatment progresses and photobleaching occurs. For each patient, the PpIX fluorescence intensity values recorded at ∼635 nm were normalized to each of their individual maximum PpIX fluorescence intensity value. Each point, therefore, in Fig. 6b represents normalized PpIX fluorescence intensity values recorded at ∼635 nm from patients at particular time points during their respective treatment times. Assuming photobleaching follows a single exponential decay curve, a time constant, τ, of 172 s represents the best fit to the data presented in Fig. 6b. The photobleaching dose constant, β, was then calculated according to Eq. 5. At a mean surface irradiance, (Ψ₀), of 82 mW/cm², a mean value of β = 14 J/cm² with an associated standard deviation of 1 J/cm² was obtained.

In the literature, a range of values of β has been provided (1.8–33 J/cm²) in various studies.^{9, 10, 27, 33} Such a wide variation in the value of β may arise from a range of factors, including the time between administration of the drug to treatment beginning, the fluence rate used for illumination, the type of animal system investigated, and the type of tumor being treated. In our case, this value was derived from clinically obtained data for human patients presenting with sBCC.

3.3.

Monte Carlo Predictions of Protoporphyrin IX Fluorescence Originating in a Tumor

As illustrated in Fig. 7a, PpIX fluorescence photons generated at deeper tumor depths have a lower probability of escape from the tumor than PpIX fluorescence photons generated at shallower tumor depths. When the treatment simulation progressed and photobleaching increased, the PpIX fluorescence emitted from the surface of the tumor decreased and more PpIX fluorescence that originated from deeper within the tumor was observed at the surface. Initially, fluorescence detected at the surface decreased rapidly with depth. As treatment progressed, the fluorescence detected from photons near the surface decreased. From halfway through treatment (75 J/cm²), almost all surface fluorescence was derived from photons deep within the tumor. The output power of PpIX fluorescence detected at the tumor surface at LD = 0 J/cm², was denoted as ^PF_{t = 0} = 6.8 ×10⁻⁵ W. The total fraction of PpIX fluorescence detected at the tumor surface at LD = 37.5, 75, 112.5, and 150 J/cm² compared to ^PF_{t = 0} was 1.95 × 10⁻³, 1.2 × 10⁻⁴, 2.77 × 10⁻⁵, and 3 × 10⁻⁶, respectively. Figure 7b illustrates the effect of increasing LD and photobleaching in the treatment simulation demonstrating a reduction in the total PpIX fluorescence detected at the tumor surface. There is a difference of almost six orders of magnitude in the PpIX fluorescence detected from the surface of the tumor between LD = 0 J/cm² (t = 0 s) and LD = 150 J/cm² (∼t = 1800 s). A decrease of approximately three orders of magnitude is evident between 0 and 37.5 J/cm², indicating a large decrease in PpIX fluorescence pertaining to rapid photobleaching followed by a slower photobleaching thereafter.

In order to compare our model, which used clinical data reported in the present paper, to clinical data reported by Cottrell,³⁴ results were reanalyzed using excitation and emission wavelengths of 632 and 705 nm, respectively [Fig. 8a]. Our simulated normalized PpIX fluorescence curve—which corresponds to 82 mW/cm²—was compared to two other fluence rates, namely, 150 and 60 mW/cm², taken from the clinical data for sBCC reported by Cottrell ³⁴ In the MC model, we used our clinically determined photobleaching dose constant, β, of 14 J/cm² and found that our simulated PpIX fluorescence data were in close agreement with the clinical data of Cottrell ³⁴ up to a fluence of ∼8 J/cm². Our simulated normalized PpIX fluorescence signal then decreased below the clinical normalized PpIX fluorescence signal corresponding to 60 mW/cm².

Figure 8b is an extension of Fig. 8a and illustrates how the normalized PpIX fluorescence continues to decrease up to 150 J/cm². At this point we found that a small quantity of PpIX fluorescence was still evident in the tumor. The data in both these graphs use the logarithmic y-axis to display the normalized PpIX fluorescence more clearly. Interestingly, Cottrell ³⁴ illustrated that PpIX fluorescence was still present up to a fluence of 200 J/cm², when treated with 150 mW/cm², suggesting that photobleaching may still be occurring.

3.4.

Monte Carlo Simulations of Singlet Oxygen Produced in a Tumor

In addition to the PpIX fluorescence detected, singlet oxygen produced at depths in the tumor was also modeled at each time step. At LD = 0 J/cm² (t = 0 s)—in other words, before photobleaching—a large quantity of singlet oxygen was evident, particularly at the surface of the tumor. As treatment progressed, the superficial layers of the tumor became photobleached and less singlet oxygen was produced in these layers. When the treatment continued further, a larger quantity of singlet oxygen could be seen to be produced deeper within the tumor than at the surface. This became more apparent during the latter stages of the treatment simulation. The majority of the PD occurs early on in PDT treatments and reduces as the treatment progresses. However, there is still a small but potentially useful quantity of singlet oxygen produced with an increasing LD in deeper layers of the tumor as shown in Table 1. This singlet oxygen produced is expressed as a fraction of the maximum singlet oxygen produced at the surface of the tumor, z = 0, and at the treatment light dose, LD = 0 J/cm², denoted as, ^SO_{LD = 0,z = 0} and equal to 5.79 × 10¹⁷ photons/cm³.

Table 1

Singlet oxygen produced at varying depths in the tumor at specific time points in the treatment simulation as a fraction of the maximum singlet oxygen produced at the tumor surface. (SOLD = 0,z = 0 = 5.79 × 1017 photons/cm3).

The singlet oxygen produced at all time points of the treatment simulation was then added together to obtain the singlet oxygen produced in a tumor over the entire treatment simulation as a function of tumor depth and LD. This is shown in Fig. 9a. The horizontal line indicates the necrosis threshold photodynamic dose, PD_T, assumed to be 8.6 × 10¹⁷ photons/cm³ absorbed by the photosensitizer per unit volume of tissue.³⁰ Therefore, we found that the maximum depth achieved of singlet oxygen produced in the tumor after administrating a total treatment LD of 37.5, 75, 112.5, and 150 J/cm² was 2.0, 2.7, 3.0, and 3.3 mm, respectively. This is illustrated in Fig. 9b. PDT with red light is an effective treatment for nonmelanoma skin cancers with a thickness of 1–3 mm. The results here are in good agreement with those quoted in the literature.^{4, 35}

3.5.

Effects of Optical Properties

The results presented in Figs. 7, 8, 9 are based on optical properties described by Salomatina,³⁰ derived from nonmelanoma human skin cancer in vitro. These optical properties appear to be well suited for the present study and thus were incorporated in the MC model. Further simulations were performed to investigate how the collected PpIX fluorescence signal and the generated singlet oxygen changed with different optical properties. These optical properties—used to represent tumor optical properties—were derived from human skin tissue in vitro, sourced from two independent studies and cover a wide spectral range, which is of interest to us.^{31, 32}

Fig. 7

(a) PpIX fluorescence detected at the surface that has originated from varying depths in the tumor with increasing LD in the treatment and fluorescence simulations and (b) total PpIX fluorescence detected at the surface of the tumor with increasing LD in the treatment and fluorescence simulations.

Fig. 8

(a) Comparison of normalized PpIX fluorescence (705 nm) detected at the surface of the tumor as a function of treatment light dose from our model using 82 mW/cm² to clinical data reported by Cottrell ³⁴ at 60 and 150 mW/cm². (b) Total normalized PpIX fluorescence (705 nm) detected at the surface of the tumor with increasing LD in the treatment and fluorescence simulations.

Fig. 9

Singlet oxygen production in units of photons/cm³ absorbed by the photosensitizer in the tumor (a) as a function of depth and LD, and (b) maximum depth achieved after a given time. The threshold photodynamic dose, PD_T, assumed to be 8.6 × 10¹⁷ photons/cm³ absorbed by the photosensitizer is signified by a horizontal line.

Figures 2a and 2b depict the curves of the optical properties used to produce Fig. 9. The scattering anisotropy factor, g, was taken as 0.8 when using the optical properties derived by Salomatina ³⁰ and 0.9 when using the optical properties derived by both Bashkatov ³¹ and Chan ³² Bashkatov ³¹ had the lowest absorption coefficients and highest scattering coefficients as a function of wavelength for the three optical property sets. Furthermore, Chan ³² had the highest absorption coefficients as a function of wavelength. Figure 10a illustrates how different sets of optical properties affected the absolute value of the collected PpIX fluorescence signal, while Figs. 10b and 10c represent there effect on the generated singlet oxygen.^{31, 32} The photobleaching dose constant, β, was fixed at 14 J/cm² for all the simulations when considering different sets of optical properties as inputs for the model. In light of this, optical property variations could explain the differences in the empirical PpIX fluorescence spectra recorded from the six patients presented in Fig. 5.

Fig. 10

(a) illustrates how different sets of optical properties affected the absolute value of the collected PpIX fluorescence signal, (b) represents the effects of optical properties obtained from Bashkatov (Ref. 31) on the singlet oxygen production, and (c) represents the effects of optical properties obtained from Chan (Ref. 32) on the singlet oxygen production.