1.

Introduction

Photodynamic therapy (PDT) is a treatment modality for malignant and nonmalignant diseases that relies on the combination of photosensitizer, light, and molecular oxygen to generate reactive oxygen species and cause local cell death. PDT is approved by the Food and Drug Administration (FDA) for the treatment of esophageal and endobronchial cancers, actinic keratosis, and Barrett’s esophagus with high-grade dysplasia and has been used to treat a number of malignancies, often at the surface of accessible organs.¹ For treatment sites that are bulky or not accessible to surface illumination, optical fibers are inserted directly into the tumor. These treatment fibers typically take the form of cylindrical diffusers, which provide roughly homogeneous emission along their distal length, and whose geometry is well suited for the treatment of large volumes. This is known as interstitial PDT (iPDT) and has been performed to treat tongue base carcinoma,² prostate cancer,³ and cholangiocarcinoma,⁴ with some applications in the treatment of large skin,⁵ esophageal, duodenal, and colorectal tumors.⁶ For the FDA-approved treatment of endobronchial cancer, where either intraluminal or interstitial illumination is possible, interstitial placement of fibers is the preferred method of light delivery due to improved efficacy and reduced exposure of healthy tissue.⁷

Since iPDT does not allow for direct visualization of the treatment light and is often performed in proximity to healthy tissue, careful treatment planning is required to ensure that a therapeutic dose is delivered to the tumor while minimizing damage to surrounding healthy tissue.⁸ This dose is a combination of photosensitizer concentration (drug dose) and fluence (light dose). For a given set of tissue optical properties, absorption and scattering, this fluence distribution can be calculated or modeled for a given source geometry and strength. Optical properties can vary significantly between and within patients, so it is necessary to determine these for each individual to ensure accurate treatment planning. Without such a determination, the portion of the tumor receiving a prescribed therapeutic dose can vary widely, even for modest errors in assumed tissue absorption and scattering.

In previous iPDT studies, optical properties have been determined by a number of methods. Researchers at the University of Pennsylvania have used isotropic fibers as sources and detectors for optical property determination in iPDT of prostate cancer.⁹ Cylindrical diffusers also have been used as sources for optical property determination, but a separate, translatable isotropic fiber was required for detection.¹⁰ Techniques involving separate spectroscopy instruments require the insertion of additional needles into the target volume, increasing the risk and complexity of the procedure. Investigators at Lund University have replaced the treatment diffusers with flat-cleaved fibers that are also used for optical property determination,¹¹ reducing the complexity. However, when flat-cleaved fibers are used for treatment, we have previously shown that treatment time and the number of fibers required to adequately cover the tumor volume are both increased.¹² It is therefore attractive to utilize treatment diffusers for spectroscopy, to take advantage of their illumination advantages for iPDT while eliminating unnecessary needle insertions.

In this study, we demonstrate recovery of optical properties at a single wavelength from turbid samples using arrays of interstitial cylindrical diffusing fibers as sources and detectors. Verification of the methods described here was performed in liquid tissue-simulating phantoms containing a scattering emulsion and ink.

2.

Materials and Methods

3.

Results

Normalized, detected fluence as a function of radial source–detector separation at the proximal diffuser ends and the resulting fit using Eqs. (1)–(3) are shown in Fig. 1(a) for a single phantom with ${\mu }_{\mathrm{eff}}=4.4{\mathrm{cm}}^{-1}$. Note that four of the detection diffusers were 2 cm in length and three were 1 cm in length, and diffuser axes were not necessarily parallel. So the detected signal does not fall off smoothly with radial source–detector separation at the proximal diffuser ends due to the increased surface area of detection for the 2-cm diffusers and varying source–detector separation at distal diffuser ends due to diffuser axes not being parallel. For the curve shown in Fig. 1(a), the extracted value of ${\mu }_{\mathrm{eff}}$ was $4.3{\mathrm{cm}}^{-1}$. The extracted values of ${\mu }_{\mathrm{eff}}$ for all phantom measurements are shown in Fig. 1(b), with the solid red line representing perfect agreement. The mean error in recovery of ${\mu }_{\mathrm{eff}}$ was 4.4%, with a maximum error of 8.9%.

Fig. 1

(a) Measured, normalized detection versus source–detector radial separation at the proximal diffuser ends for measurements made in a single phantom with ${\mu }_{\mathrm{eff}}=4.4{\mathrm{cm}}^{-1}$ (filled circles). The solid line represents the calculated detection pattern for the fitted value of ${\mu }_{\mathrm{eff}}=4.3{\mathrm{cm}}^{-1}$. Note that four of the detectors were 2 cm in length while three were 1 cm, and diffuser axes were not necessarily parallel. So the detected signal does not fall off smoothly with radial source–detector separation at the proximal diffuser ends. (b) Fitted ${\mu }_{\mathrm{eff}}$ values (open circles) determined for measurements made in 15 liquid phantoms show good agreement with known phantom optical properties. Data points shown are mean values for three repeated measurements in each phantom, with errors bars representing standard deviation being smaller than the data points in all cases. The solid line represents perfect agreement.

Extracted ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ values are shown in Fig. 2. For the recovered values of ${\mu }_{\mathrm{a}}$ shown in Fig. 2(a), each data point corresponds to the mean value determined from phantom measurements, in which the absorption coefficient was fixed at one of five values, but the amount of scatter varied over three concentrations. Similarly, each data point in Fig. 2(b) represents the mean ${\mu }_{\mathrm{s}}^{\prime }$ value found for phantoms with equal reduced scattering coefficient but varying absorption. Mean error in the recovery of absorption coefficient was 12.6%, with a maximum error of 23.9%. Mean error in the determination of reduced scattering coefficient was 7.6%, with a maximum error of 20.7%.

Fig. 2

Recovery of (a) ${\mu }_{\mathrm{a}}$ and (b) ${\mu }_{\mathrm{s}}^{\prime }$ from measurements made in liquid phantoms. Each data point in (a) corresponds to the mean value determined from phantom measurements, in which ${\mu }_{\mathrm{a}}$ was fixed at one of five values, but ${\mu }_{\mathrm{s}}^{\prime }$ varied over the three values shown in (b). Each data point in (b) represents the mean ${\mu }_{\mathrm{s}}^{\prime }$ value found for phantoms with equal reduced scattering coefficient, but with the varying values of ${\mu }_{\mathrm{a}}$ shown in (a). Error bars in both cases are standard deviation, with solid lines showing perfect agreement.

For recovered optical properties, the error tended to increase with decreasing transport albedo [$a^{\prime }={\mu }_{\mathrm{s}}^{\prime }/({\mu }_{\mathrm{a}}+{\mu }_{\mathrm{s}}^{\prime })$], as shown in Fig. 3. Error in recovery was negatively correlated with $a^{\prime }$ for both ${\mu }_{\mathrm{a}}$ (Pearson $r=-0.59$) and ${\mu }_{\mathrm{s}}^{\prime }$ (Pearson $r=-0.33$).

Fig. 3

The relationship between error in recovery of ${\mu }_{\mathrm{a}}$ (filled circles) and ${\mu }_{\mathrm{s}}^{\prime }$ (open squares) and the transport albedo ($a^{\prime }$) shows a trend of increasing error with decreasing transport albedo. Error in recovery was negatively correlated with $a^{\prime }$ for both ${\mu }_{\mathrm{a}}$ ($r=-0.59$) and ${\mu }_{\mathrm{s}}^{\prime }$ ($r=-0.33$).

Error in the recovery of ${\mu }_{\mathrm{eff}}$ increased with a decreasing number of detected signals used in the fitting algorithm, as shown in Fig. 4. When at least three detectors were included in fitting, the mean error did not exceed 6.5%, compared to a mean error of 4.4% when all seven detectors were used. For the case of only two detectors, the mean error increased markedly to 54.6%, perhaps representing a fundamental limit for this technique.

Fig. 4

The error in recovery of ${\mu }_{\mathrm{eff}}$ increases with a decreasing number of detector fibers utilized. Data points are means of results from all experiments, with the algorithm run for a given number of detector fibers, while error bars are standard deviation.

4.

Discussion

A number of research groups have examined the recovery of optical properties from interstitial measurements based on steady-state reflectance spectroscopy. Chin et al.²¹ demonstrated accuracy to within $\sim 20\%$ for measurements made in phantoms containing Intralipid and Naphthol Green using relative steady-state reflectance measurements. Dimofte et al.¹⁰ showed recovery of optical properties using diffusers as sources, but with isotropic point detectors. When the value of ${\mu }_{\mathrm{s}}^{\prime }$ was not fixed a priori, the standard deviation for determining ${\mu }_{\mathrm{a}}$ in a homogeneous phantom was 18%. In our previous work using a custom, interstitial optical probe, we were able to achieve a mean error of 9% in the recovery of absorption and 19% in the recovery of scattering.²² The mean errors of 12.6% and 7.6% in determination of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$, respectively, that we demonstrate here compare favorably to these results, without requiring separate fibers for spectroscopy.

As shown in Fig. 3, the error in the recovery of optical properties increases with decreasing transport albedo. This is expected, as the diffusion approximation utilized here is valid only for values of $a^{\prime }$ that are close to 1, with a common cutoff being $a^{\prime }>0.98$.²³ If phantoms for which the transport albedo was $<0.98$ are excluded, the mean error in recovery improves to 11.1% for ${\mu }_{\mathrm{a}}$ and 7.0% for ${\mu }_{\mathrm{s}}^{\prime }$. While this is not a large improvement in the accuracy of the algorithm, it does motivate a possible examination of higher order approximations to the Boltzmann transport equation for representation of diffusers. One such model is the ${\mathrm{P}}_3$ approximation, which expresses the terms in the Boltzmann equation using a third-order expansion in spherical harmonics. We have previously demonstrated the utility of this approximation in the context of reflectance²⁴ and interstitial fluorescence spectroscopy.¹⁶ Applying this technique to the current problem would require representing the diffuser emission and detection as a superposition of ${\mathrm{P}}_3$ point sources. While this could improve accuracy, the computational complexity would also be increased. Even with the simpler diffusion model described here, the optical property recovery algorithm requires $\sim 8\min$ to determine ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ for a single set of measurements on a modern desktop computer. Using the ${\mathrm{P}}_3$ approximation would further slow this process, possibly precluding the technique from use in a clinical setting. As an example, for the simple case of steady-state, spatially resolved, and surface-contact diffuse reflectance, the use of the ${\mathrm{P}}_3$ approximation results in $\sim 50\%$ longer runtime compared to the diffusion approximation. Additionally, the straightforward determination of $D$ from the ratio of calculated and measured data would not be possible due to the multiple exponential terms in the ${\mathrm{P}}_3$ point-source representation.²³ Development of an algorithm utilizing the ${\mathrm{P}}_3$ approximation would require further investigation to determine its feasibility and potential improvement in accuracy.

The impact of attenuation by the other diffusers was ignored in Eqs. (1)–(3). This would be the case in which light reaching a remote detection fiber has passed through another diffuser. We have neglected this for two reasons, the first being that absorption is negligible in the diffusive regions of these fibers. So the attenuation due to diffuser absorption is expected to be low. Scattering through intervening diffusers is also a potential issue. For the diffusers used in this study, we have previously shown that the scattering coefficient within the diffusive region is small for all but the most distal regions.²⁰ Using the Beer–Lambert law, this results in 95% of light remaining unscattered during transmission through the diffusive region for all but the most distal 5 mm. Since scattering within the diffuser is assumed to be isotropic due to the scatterer size, some amount of this scattered light will continue propagating in a similar direction to the incident light. Upon re-entering the surrounding turbid medium, these photons will again experience multiple scattering, resulting in the change in propagation direction being effectively negated over short length scales in the sample. Therefore, we expect the effects of diffuser scattering to be minimal.

The accuracy of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ separation described here was based on calibration of the absolute detected signal in phantoms with the same index of refraction as those used for optical property extraction experiments. This technique would therefore not be as accurate in the determination of optical properties for a sample of unknown refractive index. Further characterization would be required to determine the effects of index on the emission and detection behavior of the diffusers, as has been performed for isotropic dosimetry probes.²⁵ It is anticipated that the recovery of ${\mu }_{\mathrm{eff}}$ would be more robust to changes in refractive index, as this determination is made based on relative measures of detected fluence, rather than an absolute signal.

The goal of the optical property determination described here is to allow for the creation of accurate, individualized treatment plans for iPDT. It is hypothesized that this treatment planning could improve tumor response and reduce damage to surrounding healthy tissue. This has been undertaken by a number of research groups, largely in the context of prostate cancer.^8,13 The methods described here could be readily integrated with treatment planning software to allow for customized therapy. We have previously demonstrated treatment planning software for iPDT that can optimize diffuser length, placement, and total energy delivered based on patient CT images.¹² The impact of measured optical properties can be evaluated using the dose volume histogram (DVH) generated by this treatment planning software, which shows what portion of the tumor volume received a particular dose. This was performed for two scenarios based on the same CT scans of a head and neck cancer patient, with a difference of $0.2{\mathrm{cm}}^{-1}$ between the assumed and measured ${\mu }_{\mathrm{a}}$. For the erroneous assumed absorption, this resulted in a plan utilizing eight cylindrical diffusing fibers each delivering 400 to $1000\mathrm{J}/\mathrm{cm}$. The simulated DVH resulting from this is shown as the solid line in Fig. 5. In a clinical environment, this difference in tissue absorption could not be accounted for without a pretreatment measurement of optical properties. When the treatment plan was updated using the correct value of ${\mu }_{\mathrm{a}}$, representing the case in which a determination of optical properties was performed pretreatment, the portion of the tumor receiving the desired light dose was improved. In this case, an incorrect assumption of absorption coefficient in treatment planning led to the portion of the tumor receiving the desired light dose being reduced from 90% to 76%.

Fig. 5

DVHs showing the percentage of a tumor volume derived from CT scans of a head and neck cancer patient who received a particular simulated fluence (light dose). The solid line represents a scenario in which the treatment plan was made assuming an erroneous ${\mu }_{\mathrm{a}}=0.2{\mathrm{cm}}^{-1}$, while the dashed line represents a scenario in which the treatment plan was made assuming the correct ${\mu }_{\mathrm{a}}=0.4{\mathrm{cm}}^{-1}$.

The methods described here are vulnerable to a number of clinical considerations that could reduce accuracy. As seen in Eqs. (2) and (3), this algorithm relies upon knowledge of the positions of the proximal and distal ends of each diffuser and assumes that the diffusers are straight between these bounds. For the experiments described here, direct measurement of these coordinates was possible. In a clinical scenario, this would not be possible. Fortunately, many of the diffusers utilized, such as the Model RD (Medlight SA, Ecublens, Switzerland), are offered with radiopaque bands at both ends of the diffusive region. This would allow for accurate determination of diffuser positioning through the use of a CT scanner in the procedure room, which is present in our hospital and in many other centers that perform iPDT. These diffusers are often contained in a clear, rigid catheter during treatment, which would ensure that linearity of the diffuser axis is maintained.

Another clinical concern is the presence of heterogeneities within the tumor volume, particularly those induced by bleeding. When the catheters containing the diffusers are inserted into the tumor volume, this could result in blood pooling around the source and detector fibers. Since the model assumes homogeneous optical properties, this localized increase in absorption would be seen only as decreased detection by other fibers. This would result in the extracted ${\mu }_{\mathrm{a}}$ being larger than expected for all regions, as the effect of this bleeding is averaged over the entire path between source and detector. Additionally, tissue optical properties can be heterogeneous, meaning that measurements made over a large array of diffusers will provide results that are a superposition of the sampled absorption and scattering values. Performing measurements with a smaller number of diffusers would shrink the volume interrogated by spectroscopy and improve sensitivity to spatially heterogeneous optical properties. As we have shown here, three detection fibers are required to achieve a mean error of $<6.5\%$ in the determination of ${\mu }_{\mathrm{eff}}$. So for a clinical scenario in which more than three fibers are utilized, subsets of these fibers could be addressed in groups of four (one source and three detectors) to coarsely sample these heterogeneities.

For iPDT applications, the number of diffusers utilized can vary widely based on the tumor volume and patient anatomy. In the treatment of large tumors of the head and neck, the number of simultaneous treatment fibers is often four.^26,27 Similarly for iPDT of the prostate, generally between 4 and 12 fibers have been employed for therapeutic illumination,^13,28 with 4 or more separate detector fibers inserted. Although similar accuracy was shown for the technique described here with fewer than seven detectors used, there was an increase in error with decreasing number of detectors. So for larger tumor volumes in which a greater number of diffusers are used, more of these could act as detectors for spectroscopy, potentially increasing the accuracy of optical property extraction. As discussed above, this would result in a larger area being sampled, reducing the sensitivity to heterogeneities. A trade-off between accuracy and localization would therefore be required in clinical applications.

The techniques described here will not allow for alteration of fiber placement after spectroscopy measurements, but will allow for modification of the illumination duration at each treatment diffuser. For the case described above and illustrated in Fig. 5, the illumination duration would need to be increased by approximately a factor of 4 to account for the measured increase in absorption, resulting in each fiber delivering 1600 to $4400\mathrm{J}/\mathrm{cm}$. In this way, we expect that the methods described here could lead to improved tumor coverage during iPDT and therefore improved outcome. Since this technique does not require the insertion of additional fibers, spectroscopic measurements could also be performed during treatment illumination by temporarily decoupling specific fibers from the treatment laser and routing them to a spectrometer via optical switching. This would allow for real-time determination of tissue optical properties in the volume being treated, which could be used to detect changes in hemoglobin oxygen saturation and photosensitizer concentration for photosensitizers with sufficient extinction at the treatment wavelength. This real-time measurement could be coupled to dynamic treatment planning software to provide safer and more accurate therapy to patients.