Superficial Radially Resolved Fluorescence and 3d Photochemical Time-dependent Model for Photodynamic Therapy

Photodynamic therapy (PDT) dosimetric tools are crucial for treatment planning and noninvasive monitoring by means of fluorescence. Present approaches consider usually a 1D problem, a simple photochemical process, or a spatially homogeneous photosensitizer. In this work, a radially resolved superficial photosensitizer fluorescence and 3D photochemical time-dependent PDT model are presented. The model provides a time-dependent estimation of tissue fluorescence and the photosensitizer and singlet oxygen 3D concentrations. The model is applied to a basal cell carcinoma treated by Metvix topical photosensitizer protocol. The analysis shows the potentiality in treatment planning and monitoring. The fluorescence results are in agreement with previous measurements. Nowadays photodynamic therapy (PDT) dosimetry remains a challenging problem [1]. Although different dosi-metric factors have been studied, the photodynamic process complexity and the dynamic variation of the parameters involved make this goal difficult to achieve [2]. Accurate predictive models would enable an optimal treatment outcome, based on a personalized treatment planning and monitoring, including a decrease in the tumor recurrence rate after PDT treatment. A thorough PDT model must take into account several underlying phenomena, occurring at different spatial and temporal scales. The most relevant ones are the optical propagation in tissue, the photosensitizer spatial-temporal distribution, and the photochemical interaction among the excitation light, the photosensitizer and the tissue oxygen [2]. Most of the recent models are focused on the use of the aminolevulinic acid (ALA) in dermatology, and include approaches for these basic phenomena, as well as the protoporphyrin IX (PpIX) total superficial fluorescence [3]. A recently reported model studied the relationship between a photosensitizer fluorescence-based dosimetry and singlet oxygen distribution during treatment [4]. This model is 1D and was applied to only healthy skin, so it does not consider the photodynamic process radial dependence, nor the influence of the pathological tissue optical properties. A very recent work proposed the use of photosensitizer fluorescence to estimate the local oxygen concentration [5]. This work proves that the steady-state fluorescence signal of AlPcS4 depends on oxygen concentration and optical fluence rate. The analysis is constrained to a simple step function fluence rate, from which molecular components involved in photochemical interaction are evaluated. As a consequence, it does not consider the geometry or the optical properties of the biological medium. Both issues have been formerly included in a Monte Carlo (MC) method to calculate the PpIX fluorescence and singlet oxygen production during treatment in superficial basal cell carcinomas (BCC) [6]. However, …

Nowadays photodynamic therapy (PDT) dosimetry remains a challenging problem [1].Although different dosimetric factors have been studied, the photodynamic process complexity and the dynamic variation of the parameters involved make this goal difficult to achieve [2].Accurate predictive models would enable an optimal treatment outcome, based on a personalized treatment planning and monitoring, including a decrease in the tumor recurrence rate after PDT treatment.
A thorough PDT model must take into account several underlying phenomena, occurring at different spatial and temporal scales.The most relevant ones are the optical propagation in tissue, the photosensitizer spatial-temporal distribution, and the photochemical interaction among the excitation light, the photosensitizer and the tissue oxygen [2].Most of the recent models are focused on the use of the aminolevulinic acid (ALA) in dermatology, and include approaches for these basic phenomena, as well as the protoporphyrin IX (PpIX) total superficial fluorescence [3].A recently reported model studied the relationship between a photosensitizer fluorescence-based dosimetry and singlet oxygen distribution during treatment [4].This model is 1D and was applied to only healthy skin, so it does not consider the photodynamic process radial dependence, nor the influence of the pathological tissue optical properties.A very recent work proposed the use of photosensitizer fluorescence to estimate the local oxygen concentration [5].This work proves that the steady-state fluorescence signal of AlPcS4 depends on oxygen concentration and optical fluence rate.The analysis is constrained to a simple step function fluence rate, from which molecular components involved in photochemical interaction are evaluated.As a consequence, it does not consider the geometry or the optical properties of the biological medium.Both issues have been formerly included in a Monte Carlo (MC) method to calculate the PpIX fluorescence and singlet oxygen production during treatment in superficial basal cell carcinomas (BCC) [6].
However, this work did not consider the inhomogeneous distribution of the photosensitizer, the PDT induced oxygen reduction, or the radial dependence of singlet oxygen.
The constraints of previous works make it difficult to get an accurate quantification of the whole photodynamic process in a real clinical setting, including a thorough evaluation of the fluorescence measured at the tumor surface.This work presents a 3D model for PDT, able to predict the temporal evolution of the superficial radially resolved photosensitizer fluorescence of a skin tumor.It also provides the radially and depthresolved photosensitizer and singlet oxygen concentrations.These components are of great interest to estimate the treatment progress and 1 O 2 -mediated oxidative damage.The model has been applied to a BCC, taking into account the current Metvix protocol specifications [7], to show the potential use in the clinical practice.The results are in agreement with previous measurements.
Light distribution in a 3D tissue can be obtained by means of the radiation transport theory (RTT) [Eq.(1), where Ir; ŝ is the specific intensity and pŝ • ŝ0 is the scattering phase function].The MC method by Wang et al. [8] was employed to solve the equation.A superficial BCC of 3 mm depth, lying on healthy tissue, was considered.The optical properties of a BCC tumor at the excitation wavelength (635 nm) are μ a 1.5 cm −1 , μ s 104.76 cm −1 , g 0.79, and n 1.5 [9].A top hat laser beam with a 0.3 cm radius, perpendicular to the tissue sample, was used to deliver an irradiance of 100 mW∕cm 2 : The inhomogeneous spatial-temporal distribution of the topical photosensitizer precursor methyl aminolevulinate (MAL) in the tissue and the consequent endogenously produced photosensitizer (PpIX) during an incubation period of 3 h were calculated from Fick's law [10].The temporal evolution of the photosensitizer precursor concentration is modeled according to Eq. ( 2).MAL initial superficial concentration, M 0 , was calculated taking into account the Metvix protocol specifications, mainly the precursor density and its molecular mass.For a 1 cm radius circular lesion, M 0 was set to 4.5 • 10 20 cm −3 .D is the diffusion coefficient, whose common value through the epidermis and dermis is 0.69 • 10 −10 m 2 ∕s.K is the permeability of the stratum corneum, which acts as the main diffusion barrier, and was set to 10 −6 m∕s, greater than that for healthy skin with an intact stratum corneum: If the photosensitizer relaxation time is short compared to its precursor diffusion time (τ p ≪ t), then the concentration of photoactive compound is proportional to the instantaneous value of the precursor concentration and can be obtained by Eq. ( 3).In this equation, ε p is the conversion process yield and τ a→p the relaxation time of the photosensitizer precursor due to the generation of the photoactive compound.Here the conversion process yield was set to 0.5, the relaxation time of the photosensitizer to 84 ms, and the relaxation time of the precursor to 25 h [11]: The photochemical pathways during a type II PDT process were modeled by the differential equations system in Eqs. ( 4)- (9).It was first proposed by Foster et al. [12], and later employed by others [13,14], as the basis to describe the temporal evolution of the molecular components involved in the photochemical interaction: In these equations, S 0 , S 1 , and T are the photosensitizer concentration in the ground state, singlet excited state, and triplet excited state, respectively. 3O 2 is the oxygen concentration in ground state; 1 O 2 is the singlet oxygen concentration; R is the concentration of 1 O 2 receptors; C is the scavengers concentration; τ1 is the relaxation time from state S 1 to S 0 ; τ3 is the relaxation time from state T to S 0 ; τ0 is the relaxation time from state 1 O 2 to 3 O 2 ; η 10 is the quantum yield of the S 1 to S 0 transition; η 13 is the quantum yield of the S 1 to T transition; η 30 is the quantum yield of T transition to S 0 ; η 0 is the quantum yield of 1 O 2 transition to 3 O 2 ; αs is the efficiency factor for energy transfer from T to 3 O 2 ; kpb stands for the biomolecular photobleaching rate; kcx is the biomolecular cytotoxicity rate; ksc is the rate of reaction of 1 O 2 with various oxygen scavengers; ν is light speed in tissue; ρ is the photon density; σ psa is the absorption crosssection of S 0 molecules; P is the rate of oxygen diffusion and perfusion; and U is the cell damage repair rate.The values to the previous parameters can be found in [11].The stiff differential equations system was solved by means of a differential equation solver within the MATLAB platform.The time-resolved 3D photosensitizer fluorescence emission at 705 nm was obtained, taking into account the excitation photon density absorbed by the photosensitizer molecules, and the fluorescence quantum yield.The fluorescence power density generated can be expressed as [15] where E photon λem h • c∕λ em is the photon energy at the fluorescence emission wavelength λ em , h is Planck's constant, and c is the speed of light in vacuum.
The superficial fluorescence due to the photosensitizer molecules was calculated as the treatment progresses.For this purpose, another MC approach was used, based on [16].The fluorescence escape was due to a timedependent spatially nonuniform fluorophore distribution.The local fluorescence sources were dynamically updated over time and space by Eq. (10).From the spatial 3D distribution of the fluorescence sources within the tissue, it is possible to obtain their contribution to the total escaping fluorescence at tissue surface.For a particular treatment time, the total escaping fluorescence, J f , at a certain superficial radial position, r, can be obtained by the accumulation of the escaping flux produced by each of the fluorescent power sources.The power density of fluorescence escaping at the tissue surface is calculated according to Eq. (11).In this equation, r s and z s are the coordinates of the fluorescence source, ΔV is the incremental volume associated to the source position, and T is the fluorescence transfer function, from source to surface: J f r; t X r s X z s P f r s ; z x ; tΔV r s ; z s Tr s ; z s ; r: (11) Figure 1 shows the total escaping fluorescence irradiance at the BCC surface as the irradiation time increases.As expected, it decays with time as a consequence of the PpIX photobleaching, which is higher in the superficial zones.This fact results in higher fluorescence at the beginning of treatment, coming mostly from the superficial zones.As the treatment evolves, the potential fluorescence sources are located deeper, causing a pathlength increase of photons to reach the surface and, as a consequence, a decrease in the PpIX fluorescence observed at the tissue surface.A slight increase in the fluorescence out of the optical beam radius can be observed, which becomes evident at approximately the halfway treatment (t 300 s).The fluorescence increase is caused by the scattered excitation and emission photons beyond the directly irradiated area.As it can be seen by comparing t 60 s with t 300 s or t 600 s, this effect intensifies as the treatment progresses, due to the increment of scattering events of fluorescence photons coming from deeper positions.These results highlight the relevance of considering the spatial optical beam profile and the optical properties of the target medium.
The model provides also the time-and spatially dependent photosensitizer and singlet oxygen concentrations.The photosensitizer concentration at different depths within the tumor is shown in Fig. 2 for different temporal instants during the irradiation period.As it can be seen by comparing t 6 s and t 60 s, the higher photobleaching in superficial zones provokes the relocation of the fluorescence photosensitizer sources at deeper positions.This is in agreement with the previous analysis of Fig. 1.Generated singlet oxygen concentration is depicted in Fig. 3.As it can be observed by comparing t 6 s and t 60 s, the maximum concentration is reached in deeper areas as time progresses.This concentration is related with tumor oxidative damage.As a consequence, the adequate and accurate estimation of this concentration is crucial for treatment planning.
Figures 4 and 5 show the relevance of the photodynamic process radial dependence.For a particular depth within the tumor, the time-dependent photosensitizer degradation (Fig. 4) and singlet oxygen concentration (Fig. 5) are evidently correlated with the radius of the optical beam (0.3 cm).The areas with increased concentration of singlet oxygen will be more affected by the oxidative damage.As it is shown in Figs. 4 and 5, this cytotoxic effect extends beyond the optical beam radius, particularly for longer times (t 600 s), due to the scattering of the biological tissue.Again it is shown that a complete PDT analysis should include the 3D optical distribution, considering the optical source profile and the tissue optical properties.This would allow the avoidance of cytotoxic effects out of the area of the tumor.
The total normalized PpIX fluorescence signal obtained by means of the proposed PDT model, applied to a BCC tumor, was represented as a function of treatment time in Fig. 6 (red line).As shown, the model is in good agreement with previously reported measurements in patients affected by the same type of nonmelanoma skin cancer [17], Fig. 6 (blue line).The difference between the blue and red lines in Fig. 6 is due to the distinct irradiance of the clinical Metvix protocol in our model, 100 mW∕cm 2 and in the measurements in [17], 150 mW∕cm 2 .Nevertheless, this difference indicates that a greater amount of irradiation (150 mW∕cm 2 ) induces a fluorescence irradiance increase (Fig. 6, blue line), as expected.This fact demonstrates the accuracy of the model to direct noninvasive fluorescence measurements.The presented model allows the prediction and interpretation of the superficial radially resolved photosensitizer fluorescence.The model also permits us to relate this information to the spatial-temporal evolution of the key molecular components involved in the development and the outcome of treatment such as the photosensitizer or singlet oxygen concentrations.Several dosimetric parameters can be varied, including the optical source radius or irradiance, the photosensitizer type, the photochemical interaction parameters or the biological media geometry and optical properties.It could be employed as a personalized planning and noninvasive monitoring tool in the dermatooncologic clinical practice in order to optimize PDT treatment and to avoid tumor recurrence.