1. Introduction

In therapeutic and diagnostic applications of light in medicine, one of the central problems is accurate determination of the optical properties (OPs) of tissues, tumors and exogenous drugs such as photosensitizers. OPs can be used to indicate whether a tissue is in a normal or pathological state [1]. Other physiological parameters such as hemoglobin oxygenation are closely related to the OPs and can be conveniently measured non-invasively [2]. OPs can also provide the basis for contrast in either structural or functional optical imaging [3]. Even if OPs are not the final goal for a particular application, they are still required to describe light propagation whether the measurement geometry is invasive or non-invasive. Calculation of the light fluence distribution and the photosensitizer concentration used in photodynamic therapy (PDT) also requires knowledge of OPs [4]. Hence, accurate knowledge of OPs is essential for optimum use of light in diagnosis and treatment of diseases.

PDT has been demonstrated to be effective for the treatment of cancer at several sites, including skin, retina, lung, and reproductive and urological tissues and organs, and has been clinically approved for some of these sites in Canada, the United States, Europe, and Japan [5]. The treatment of deep-seated tumors such as prostate cancer requires that optical fibers be implanted within the organ. In this interstitial geometry, a spatially resolved steady state technique has been developed to determine the OPs in phantoms and in-vivo, where light was delivered and collected by isotropic fibers [6–8]. As opposed to this steady state approach, in the frequency domain measurements can be undertaken at only two positions to recover OPs (the absorption coefficient, µ_a and the reduced scattering coefficient, µ_s’) since the additional information of phase can be obtained. Figure 1 shows the detection geometry used in this paper. Two isotropic detectors at different radial positions collected light emitted from an isotropic, sinusoidally-modulated source. All fibers were immersed in a phantom sufficiently large to mimic the clinical situation. Optical properties were determined from the heterodyne frequency domain measurements of the phase difference and the ratio of the amplitudes (fluences) at two known positions, using inverse calculation based on either the combined isotropic similarity model (ISM) and standard diffusion approximation (SDA), ISM/SDA for short, [9] or the modified spherical harmonics method [10–11].

 

Fig. 1. The geometry of the two-distance detection scheme. Two detector fibers are placed parallel to each other as well as to the source fiber. For these experiments, a single detector fiber was translated between positions.

Download Full Size | PDF

This article aims to examine the feasibility of this detection scheme and the accuracy in the recovered OPs by experiments in tissue-simulating phantoms. Our previous theoretical studies showed the potential improvement in recovering OPs using the ISM/SDA model, short for the combined isotropic similarity model and standard diffusion approximation [9] instead of the standard diffusion approximation. We will first briefly review another interesting method for solving the radiative transfer equation, the modified spherical harmonics method [10–12], followed by the experimental comparison of the performance of the two models. Source anisotropy and inaccuracy in fiber placement are commonly encountered practical issues and these were investigated through Monte Carlo simulation. Phase amplitude crosstalk in frequency domain systems has received little attention in the literature. In contrast to reports of other researchers, our experimental results showed that this crosstalk is not significant. Misleading results can be obtained by some methods of measuring this crosstalk. Finally, we compared estimates of OPs made under two conditions: 1) two detectors located at 4mm and 7mm, or 7mm and 10mm from the source without a neutral density filter used to attenuate the signal for the detector closer to the source, and 2) two detectors located at 5mm and 10mm from the source with the detector closer to the source equipped with a neutral density filter.

2. Theoretical background

The radiative transport equation (RTE) is widely accepted as an accurate description of light transport in tissue but an exact solution can be found only under restricted conditions. We need a solution for the fluence rate in an infinite medium for which the scattering is highly forward directed and where the intensity of the source is sinusoidally modulated. In a previous paper [9] we developed an approximate solution we will refer to as the combined isotropic similarity model (ISM) and standard diffusion approximation (SDA), ISM/SDA for short. For the inverse problem described above, we found that the SDA recovers the reduced scattering coefficient better than the absorption coefficient, while the opposite is true for the ISM model. Hence in the ISM/SDA, we used ISM to estimate the absorption coefficient and SDA to recover the reduced scattering coefficient. The solution of the diffusion equation for a sinusoidally modulated point source in an infinite turbid medium is fully described in the literature [13] and will not be repeated here. Here we present only the basic idea of the ISM model and the final result for the fluence rate. The first step is to use the principle of similarity to replace the tissue medium with an “equivalent” medium in which light is scattered isotropically. To do this, the absorption coefficient is left unchanged but the scattering coefficient µ_s is replaced with an isotropic scattering coefficient µ′_s =(1-g)µ_s where g is the average cosine of the scattering angle in the tissue, typically greater than 0.9. The second step is to derive an analytic solution for the fluence rate in this equivalent medium due to an isotropic, sinusoidally modulated source. The fluence rate, φ, can be thought of as complex variable where the phase is related to the propagation delay from the source to detector. In reference [9] we developed the following expression for the real and imaginary parts of the complex fluence rate:

where µ′_t =µ_a +µ_s ′; c, the albedo is µ′_s /µ′_t ; the real parameter a ₀ is the root of c=a ₀/tanh-a ₀; k ₀ is imaginary and defined as k ₀=ia ₀ ; r is the distance from source to detector; and the functional form of g(c,ξ) is 1/[(1-cξ tanh^-1 ξ)²+(πcξ/2)²]. The imaginary part, Im(φ), is:

where v is the speed of light in tissue, ω is the angular frequency of the source modulation, and k′ is a dummy integration variable.

An alternative to the ISM/SDA is the modified spherical harmonics method (MSH) [10,11]. It has been found that this method is fast and accurate, thus it is particularly suitable for the inverse calculation because fluence and phase can be expressed in an analytical form [12]. The essence of this method is that the angular part of the Fourier transform of quantities such as the radiance and the phase function in the radiative transfer equation (RTE, also known as the Boltzmann equation), can be expanded in the basis of spherical harmonic functions defined in a reference frame (called the rotated frame) whose z-axis is aligned with the direction of the Fourier wave vector. Because of this specific choice of the basis spherical harmonic functions in the rotated reference frame, the RTE is simplified, compared to the RTE in the laboratory frame. Solving the RTE in the rotated frame reduces to solving the eigenvalues and eigenvectors of a tridiagonal symmetric matrix W. The accuracy of this method depends on the chosen order of W. In the original paper [10], it was argued that for typical biological tissue, the accuracy needed can be satisfied with order less than 50. The computation of the W matrix of this order or less can be done within seconds using Matlab code on a personal computer. In general, MSH is an efficient and accurate numerical method for the RTE in infinite geometry. The formalisms relevant to our work can be found in Appendix A. Since we are interested in the value of the reduced scattering coefficient µs’ instead of the individual scattering coefficient µs and anisotropy factor g, all the calculations in this article used g=0.9 when MSH was employed. Use of the preset value g in our particular two-distance detection scheme is required because this scheme can recover only two independent optical parameters instead of the three parameters (which are µ_a, µ_s’ and g).

3. Experimental methods

3.1. The frequency domain system

Standard heterodyne detection was used in our frequency domain system. Two frequency generators (Marconi Instruments, Model 2022A) were used to provide modulation (typically 100MHz) of the laser source (Oriel Diode Lasers, Model no. 79401) and the high voltage (frequency was typically offset by 200Hz compared to the frequency of the laser modulation) applied to the photomultiplier tube (PMT) (Hamamatsu R928) with the housing (ISS model no. P016). The wavelength of the diode laser is 750nm with maximum power of 4.3mW and 100% modulation depth. Source and detector fibers (0.4mm in diameter, Photoglow) are uniform within ±5% on 88% (symmetrical to the fiber line) of their surfaces. The signal detected from the probe fiber was fed into a lock-in amplifier (Stanford Research Systems model SR850) via a signal preamplifier (Ithaco model 1642). In this analysis we refer to two detector fibers, but for these experiments a single fiber was translated between positions using a stage with 0.02mm resolution.

3.2. Data acquisition

Data collection in the frequency domain measurement was performed with LabView (National Instruments, Austin, TX) through the RS-232 serial port of the lock-in amplifier. For each set of fluence and phase measurements, nine runs were collected and subsequently averaged to reduce the noise and obtain the standard deviation. All measurements were performed with the room lights off. The RF amplifier was positioned about two meters from other instruments to reduce electromagnetic interference.

3.3. Monte Carlo simulation

A detailed description of Monte Carlo simulation can be found in a previous publication [14]. In our work, we allowed for anisotropic light scattering by utilizing the Henyey-Greenstein phase function [15]. Light uniformly emitted from the source (either for the whole or a portion of the 4π solid angle) was simulated using a uniform sampling scheme. The fluence scoring mesh was either spherical shells with 0.1mm radial separation for an isotropic source or spherical shells with additional polar angular resolution of 1 degree for a non-isotropic source (which still preserves azimuthal symmetry, see discussion below). Fluence is just the deposited photon weight per unit volume. The phase is calculated by tracking each photon in the time-domain and Fourier transformation of the total time-dependent fluence rate to the frequency domain.

3.4. Inverse algorithm

The inversion scheme to estimate OPs from the frequency domain data is the same as that previously described [9]. Fluence and phase data and trial OPs were used to calculate χ ² using the selected forward model. The trial OPs were adjusted until 15 $\frac{{\chi }_{\mathit{prev}}^2-{\chi }_{\mathit{cur}}^2}{{\chi }_{\mathit{cur}}^2}\mid <{10}^{-5}$, where ${\chi }_{\mathit{\text{cur}}}^{\mathit{2}}$ and ${\chi }_{\mathit{\text{prev}}}^{\mathit{2}}$ refer to the values for the current and previous iterations. The value of 10–15 was selected because further reduction yielded insignificant change (less than 0.1%) in estimated OPs. χ² is defined as: ${\chi }^2={\left(\frac{R_m-R_c}{\delta R_m}\right)}^2+{\left(\frac{{\varphi }_m-{\varphi }_c}{\delta {\varphi }_m}\right)}^2$, where R (ϕ) is the ratio of the fluences (phase difference) at two known distances; δR_m (δϕ_m ) is the standard deviation of the fluence (phase); m and c denote measured and calculated respectively. The Levenberg-Marquardt optimization algorithm [16] was used to minimize χ² .

3.5. Accurate measurement of OPs of Intralipid and India ink

The tissue-simulating phantom used in this work was made using various concentrations of Intralipid (IL 20% Pharmacia Corporation, Peapack, NJ) and India ink (Demco, St-Lambert, Quebec, Canada). Determination of OPs of Intralipid and India ink is not a trivial task. Some of the previous work was summarized in the paper by Dimofte et al [8]. They showed that there exists a significant variation in the results reported by different investigators. The difference could be as large as a factor of two (e.g. µ_s’ values from van Staveren [17] are twice those from Flock [18]). The discrepancy can be partly ascribed to the variation in batches. It is also pointed out by Flock [18] that appropriate measurement techniques must be used for characterization of OPs of phantoms. In our work, narrow beam attenuation was employed to determine the absorption coefficient (µ_a) of India ink using a fiber optic spectrophotometer (Ocean Optics, SD2000) in transmission mode. In order for the ink concentration to fall in the linear region of Beer’s law, the absorbance is usually chosen to be 0.05 to 0.7. The measurements for the absorption for the India ink at 750nm can be described by:

The intercept indicates the absorption due to water (milliQ water), which is close to the published data: 0.00261mm^-1 [19]. It has been shown that India ink may have non-negligible scattering which can result in underestimation of scattering or overestimation of absorption of ink-based phantoms [20]. This scattering is due primarily to large particle aggregates (~1µm), which can scatter light significantly. We ignored the scattering component due to ink and show below that this assumption was justified.

The added-absorber method was employed to determine the reduced scattering coefficient (µ_s’) of Intralipid and check the values obtained above for µ_a of ink. For broad beam irradiation of a turbid medium the fluence rate varies with depth z as:

where ${\mu }_{\mathit{eff}}=\sqrt{3{\mu }_a\left({\mu }_a+{{\mu }_s}^{\text{'}}\right)}$ under conditions where standard diffusion theory applies. In order to obtain a satisfactory uniform beam, the tip of an isotropic fiber (which can be regarded as an isotropic point source) was located at the focal point of a lens with a diameter of 10cm. A second isotropic fiber was placed within the phantom and used to measure the depth dependence of the fluence rate. The phantom container used was 15×15×11cm³. Care must be exercised because the phantom volume cannot be regarded as (semi) infinite if 1/µ_eff is less than about 12 times its diameter, in which case light losses at the boundaries cause an overestimation of µ_eff [21]. At the same time, µ_a should be kept less than 5% of µ_s’ in order for diffusion theory to be valid. Assuming the absorption of Intralipid itself is negligible, the following equation results:

where ${\mathit{\mu }}_a^{\mathit{\text{ink}}}$ =ξL(ξ is a constant in mm^-1, L is the concentration of the India ink in percent) and ${\mathit{\mu }}_a^w$ =0.00261 mm ^-1. By fitting the measured (µ_eff)² to a range of concentrations of added ink (0.02% to 0.04%) for a phantom with fixed µ_s’, the intercept gives the reduced scattering of IL and the slope gives the absorption of ink. The results are: µ_s’=1.034mm^-1 for 1% IL and ξ=3.09mm^-1 at 750nm, which agrees with the ink formula obtained in the spectrometer measurements (Eq. 3). By varying the IL concentration (0.3% to 2.0%) with fixed ink absorption (0.01mm^-1) and repeating the added absorber experiment, we determined that µ_s’ due to IL is:

at 750nm.

3.6. Phase amplitude crosstalk

There is limited discussion of phase amplitude crosstalk in the literature [22–26]. Ramanujam et al [22] showed that this phenomenon occurs in both homodyne and heterodyne detection techniques using either PMT (R928) or avalanche photodiode (APD) detectors. According to the paper by Ramanujam et al [22], the crosstalk was 12.23o/O.D. (i.e. factor of 10 change in amplitude) using an R928 PMT (high voltage was 740V) and heterodyne technique (modulation frequency was 200MHz). The accepted explanation is that detector rise time changes with input light intensity [25]. The methods proposed to rectify this problem were 1) maintenance of a fairly constant amplitude at the output of the detector using a dynode feedback loop, or 2) use of a correction curve as the crosstalk is an intrinsic and reproducible property of the detector [26]. However, we found that the measured crosstalk effect is dependent on the geometry of the light path and the method by which light is attenuated.

 

Fig. 2. Phase shift vs. amplitude at high voltage of 700V on PMT. The light was attenuated using a variable optical attenuator. Error bars show the standard deviations.

Download Full Size | PDF

Figure 2 shows the relationship between the phase shift and amplitude when the amplitude was reduced at the output of the laser using a fiberoptic coupled variable optical attenuator (OZ Optics). The high voltage of the PMT was set to 700V. The phase shift can be as large as 2.5o if the PMT anode current is reduced from 0.35 to 0.05µA. This curve qualitatively agrees with that obtained by other groups such as Ramanujam et al [22]. We next used neutral density filters between two collimating lenses to attenuate the laser light. The distance between the two lenses was adjustable and light attenuation was achieved by inserting an appropriate neutral density filter (NDF) between them. It should be noted that because the index of refraction of the NDF is different from that of the air, a phase change is introduced by the filter. This effect has been taken into account by a correction factor determined by the thickness and refractive index of the NDF, modulation frequency, and the light speed. Figure 3 shows the phase shift due to the attenuation through four neutral density filters (O.D.=0.38, 0.43, 0.92, and 0.99) when the separation between the two lenses was 3.0cm and the PMT anode current with no filters was about 0.35µA. It can be seen that the phase shift steadily increases as attenuation increases. One order of magnitude reduction in amplitude results in about 0.8° phase shift if the curve is projected back to optical density equal to 0. We also examined how the phase shift behaved if the separation between the two lenses was varied. Figure 4 shows the phase shift due to four different separations (0.0, 3.1, 5.5, and 8.0cm) between the two lenses. The O.D. of the neutral density filter used was 0.92. Figure 4 shows that larger separation gave a larger phase shift, and the phase shift was less than 0.20° if the lenses were not separated. One possible reason that the optical attenuator gave more phase shift than did the NDF is that the attenuator blocks part of the laser beam and the NDF attenuates the beam uniformly. This could result in preferential attenuation of some modes of the optical fibers and a change in the observed phase. For our setup, we concluded that phase amplitude crosstalk was comparable to phase noise, and we have ignored it in the analysis.

 

Fig. 3. Phase shift due to the attenuation through four neutral density filters (O.D.=0.38, 0.43, 0.92, and 0.99). Separation between the two collimating lenses was 3.0cm. The output current at the PMT was maintained at 0.35uA. Error bars show the standard deviations.

Download Full Size | PDF

 

Fig. 4. Phase shift caused by the insertion of 0.92 OD neutral density filter at four different separations (0.0, 3.1, 5.5, and 8.0cm) between the two collimating lenses. Error bars show the standard deviations.

Download Full Size | PDF

4. Results

4.1. MC simulation of anisotropic source

As stated earlier, light intensity on a spherical surface centered at the tip of the fibers used in our work is only uniform (within 5%) for the polar angle between 40° and 180° (see Fig. 5). The detailed light intensity profile inside the cone (between 0° and 40°) of the fiber is not provided by the manufacturer. This anisotropic source will cause the fluence rate to depend not only on the distance to the measurement point but also its angular coordinates. In highly scattering media this dependence will be less evident as the source-detector distance increases. This was demonstrated by Monte Carlo simulations for three sets of OPs, (µ_a, µ_s’)=(0.01, 1), (0.08, 0.4), and (0.2, 0.3) in mm^-1, which cover µ_s’/µ_a from 1.5 to 100. In the Monte Carlo simulation, the source was a sphere of diameter 0.4mm. Photons were launched from random locations on the sphere surface, outside the 40° cone. No photon was launched inside that cone (see the insert in Fig. 5). Detector fibers were assumed to have isotropic response to simplify the MC code. The scattering anisotropy factor was set to 0.9. For the medium with (µ_a, µ_s’)=(0.01, 1.00)mm^-1, Fig. 5 shows the typical variation in fluence rate versus polar angle at a distance of 3mm from a unit strength source. The insert in Fig. 5 shows the source orientation and definition of angle. The fluence rate was 0.042mm^-2s^-1 at zero degrees (i.e. a point “behind” the fiber) and steadily increased to 0.049mm^-1 at 180°. In order to examine how this variation in the fluence rate affects the result of inverse calculation, three pairs of distances (2, 3), (2, 8), (4, 8)mm were chosen. The detectors were located as shown in Fig. 1. The choice of this geometry is based on the fact that in most clinical PDT, fibers are inserted with the guide of a brachytherapy template in which fibers are parallel to each other. The expected recovered OPs, (µ_a, µ_s’)=(0.08, 0.4) mm^-1, are summarized in Table. 1. Similar results were observed for the other sets of OPs (data not shown). The pair (4, 8) using MSH model shows OPs can be recovered within 3%. For the same positions, ISM/SDA gives about 8% error in estimating µ_s’. For the other two position pairs, the errors in recovered OPs can be as large as 30% (using ISM/SDA) and 12% (using MSH). Hence, in order to minimize the effect of anisotropy of the source fiber, the detection positions should be at least 3mm away from the source.

 

Fig. 5. Fluence rate (mm^-2s^-1) versus polar angle defined in the insert. Fluence rate was calculated at 3mm from the source. In the insert, the bold line represents the fiber, and the angle denoted by θ varies from 0 (positive z direction) to 180 (negative z direction).

Download Full Size | PDF

Table. 1. Recovered OPs using ISM/SDA and MSH for three pairs of positions. The expected OPs: µs’=0.4mm-1, µa=0.08mm-1

View Table | View all tables in this article

4.2. Experimental comparison of MSH and ISM/SDA models under “extreme” conditions

Previous studies showed that ISM/SDA can recover OPs from MC generated data within 5–10% for µ_s’/µ_a>5. Monte Carlo simulation performed above confirmed that point and has also shown that MSH can recover OPs very accurately for µ_s’/µ_a>1, and the computation can be done rapidly (within seconds). In order to investigate the effect of factors that MC simulation did not incorporate such as non-uniform response of detector fiber, noise in instruments, perturbation to fluence induced by the presence of the fibers, we performed experiments with tissue-simulating phantoms. In this part of the study the reduced scattering coefficient was chosen to be 0.2mm^-1. Higher values could be used, but high absorption then has to be used in order to obtain low albedo and the signal is not measurable at larger distances. The absorption coefficients were set at 0.02, 0.04, 0.1, 0.2, 0.3mm^-1 by adding appropriate amount of ink to the phantom. The detection position pair was chosen to be (5, 8)mm. The results for the recovered OPs are shown in Table. 2. It can be seen that the performance of MSH is comparable to ISM/SDA for µ_s’/µ_a>5. For the phantom with µ_s’/µ_a=9.6, the errors caused by both methods are greater than (or close to) 10%. Since MC simulation has shown that for this µ_s’/µ_a ratio, the recovered OPs using both models should be in error by less than 10% if the noise in amplitude and phase is about 10% [9] comparable to our experimental noise, we suspect that there was some other uncertainty introduced in the data collection process. The most likely uncertainty is in measuring the distances of the detectors from the source, as will be discussed later in this section. MSH works well for 1<µ_s’/µ_a<5, while ISM/SDA can introduce errors as large as 23.6% in recovered OPs. For the phantom with µ_s’/µ_a<1, generally MSH introduces errors greater than 10% (1.0% error is probably fortuitous), but still works better than ISM/SDA for which the error is generally greater than 20%. For media with µ_s’/µ_a<1, MSH should work satisfactorily in principal. The reason that fairly large errors (>15%) occur is probably due to the inaccurate measurement of small phase. The experimental phase was as small as 0.6°, which is not substantially greater than its standard deviation 0.18° (based on our data). To increase the phase entails an increase in detector separation, which in turn will dramatically increase the difference between the amplitudes measured at the two detection positions. Fortunately, the ratio µ_s’/µ_a encountered in biological tissues, even those with high blood content (implying high absorption), is greater than 5. Use of either MSH or ISM/SDA should be satisfactory in most cases.

Table. 2. Experimental comparison of MSH, and ISM/SDA. Detection position pair was (5, 8). The expected µ_s’=0.2mm^-1

View Table | View all tables in this article

4.3. MC study of inaccurate fiber placement

Needles inserted into tissue for optical fiber placement may bend, so that even if a template is used, the real positions of the needles and optical fibers may not be those expected. Hence an analysis of the resulting errors in estimated OPs due to the misplacement of fibers is desirable. First, we found that the recovered OPs can be off by at least 10% if one of the fibers is misplaced by 0.1mm and the other is accurately positioned. The Monte Carlo generated data for (µ_a, µ_s’) were OP1: (0.005, 1.00), OP2: (0.015, 0.9), and OP3: (0.2, 0.5) in mm^-1. The recovered OPs were obtained by offsetting the probe at 7mm (the other probe is at 10mm) by 0.1, 0.2, 0.3, 0.4, 0.5mm away from the source and using those data in the inverse algorithm. The results are shown in Table 3 and can be understood as follows. It is known that the fluence decreases sharply with distance close to the source and more slowly beyond around 10mm (the actually transition distance depends on the OPs). Hence small changes in distance close to the source result in relatively large changes in fluence.

Table 3. Monte Carlo study of relative difference of recovered OPs using ISM/SDA model in µ_a (left) and µ_s’ (right).

View Table | View all tables in this article

The probe at 7mm was offset by 0.1, 0.2, 0.3, 0.4, and 0.5mm away from the source. The other probe is fixed at 10.0mm. OP. 1, 2, and 3 are (1.00, 0.005), (0.90, 0.015), and (0.50, 0.20) in mm-1, where 1^st (2^nd) number refers to absorption (reduced scattering).

Second, we investigated the error in recovered OPs in the situation where the separation between the two fibers is fixed at 3mm but both fibers are displaced. The expected OPs were µ_a=0.01mm^-1, µ_s’=1.00mm^-1. The nominal positions are (3, 6), (5, 8) and (7, 10) mm. Table. 4. shows the recovered OPs for the nominal positions of (3, 6) when both detectors were shifted by 0.1, 0.3, 0.5 1.0 1.5mm. The upper row in µ_a and µ_s’ corresponds to the two detectors shifted away from the source. The lower row in µ_a and µ_s’ shows the two detectors shifted toward the source. Two points can be observed: 1) the same shift will give larger error closer to the source. The shift starts to give significant error (10%) if it is more than 0.3mm when close to the source, and more than 1.0mm far from the source; 2) the error in recovered OPs due to the simultaneous shift of the two detectors is far less than that due to the offset of one fiber (with the other in the correct position). Hence the two-distance technique has the advantage that precise measurement of the absolute position of the two fibers is not required, but the relative distance between the two fibers is important. It also shows that the fibers should be placed relatively far from the source (point (1) above-mentioned) to minimize the errors.

Table 4. Monte Carlo study of the effect of the shift of two detectors on recovered OPs using ISM/SDA model in µa and µs’.

View Table | View all tables in this article

The upper row in µ_a (µ_s’) corresponds to the two detectors shifted away from the source. The lower row in µ_a and (µ_s’) shows the two detectors shifted toward the source.

4.4. Choosing OPs of prostate tissue-simulating phantoms

The in vivo optical properties of prostate demonstrate fairly large inter- and intra- animal variability compared to other tissues such as brain. The absorption and reduced scattering coefficients vary between 0.003–0.058 and 0.1–2.0 mm^-1 (human prostate) at 732nm according to Zhu et al [6], and between 0.002–0.053 and 0.4–2.0mm^-1 (canine prostate) at 660nm as reported by Lilge et al [27]. During PDT, the presence of photosensitizer contributes to the absorption as well. Some drugs such as Tookad have high molar extinction coefficient (~10⁵ M^-1cm^-1) at the peak absorption wavelength [28]. All these form the basis for our choice of OPs used for the phantoms.

Table 5. OPs for the phantoms used in Figs. 4, 5 (top two rows), and 6 (bottom two rows).

View Table | View all tables in this article

Table 6. The average of root-mean-square percent errors and standard deviations (in bracket) in recovered OPs for nine phantoms (for OPs, see text) using ISM/SDA and MSH.

View Table | View all tables in this article

4.4.1. Case where a neutral density filter was not used

We have chosen three pairs: location 1: (4, 7) mm, location 2: (7, 10) mm and location 3: (5, 10) mm. The first and the second number in the brackets are distances of the first and second probe locations respectively, measured from the source fiber. The OPs are chosen as follows (for reasons, see above): µ_a=0.008, 0.023, 0.17mm^-1, and µ_s’=0.3, 1.05, 1.8mm^-1, out of which nine possible combinations can be formed and Table 5. shows how they are numbered. The measurements were not conducted for phantom #9 at location 3 due to the weak signal at the farther detector. Figures 6 and 7 show the percent errors in the OPs recovered using ISM/SDA and MSH respectively at the three locations for nine phantoms. To better understand which location pair is best, we calculated all the root-mean-square (RMS) percent errors (for repeated runs on the same phantom) and their standard deviations for all nine phantoms at each location pair, (see Table 6). The standard deviations are measures of how the errors spread among nine phantoms for each pair location. The errors at location 2 are the smallest: 13% in µ_s’ and 9% in µ_a on average. Thus two detectors should be placed “far” from the source to minimize errors in the estimated OPs.

 

Fig. 6. Percent errors in estimated OPs using ISM/SDA model for nine phantoms at locations 1, 2, and 3 that are shown by sub-Fig (a), (b) and (c) respectively. Phantom identification numbers can be found in Table 5.

Download Full Size | PDF

4.4.2. Case where a NDF was used for the detector closer to the source

In a relatively high absorption turbid medium, the phase shift between two detectors is small if their separation is not sufficient. This dependence of phase shift, θ, on µ_a, r₁, r₂ can be expressed analytically using the diffusion approximation:

On the other hand, large separation, which gives larger phase shift, entails much weaker signal at the second detector, especially in high absorption media. This in turn will introduce relatively large error in the amplitude ratio. One simple way to obtain large phase shift and reduce the dynamic range of the amplitude measurement is to use a calibrated NDF to attenuate light reaching the PMT when it is detecting light from the near fiber. The correction to phase shift due to the difference in the speed of light between the NDF and the air was taken into account. The OPs chosen for this experiment are listed in Table. 5: µ_a and µ_s’ range from 0.008mm^-1 to 0.10 mm^-1, and 0.3mm^-1 to 1.6mm^-1 respectively. The O.D. of the NDF was 0.90±0.02. The two detectors were located at 5 and 10mm from the source. Figure 8 shows the error in estimated OPs using ISM/SDA [sub-Fig. (a)] and MSH [sub-Fig. (b)] models for the nine phantoms. Again, it can be seen that errors in µ_a were less than those in µ_s’ in most cases. Using the ISM/SDA model, the errors in recovered OPs were greater than 10% for phantom 3, 5, 6, and 9. The error in µ_a for phantom 5 is slightly greater than 10%. The other three phantoms show larger errors because: 1) ISM/SDA breaks down for phantom 3 because µ_s’/µ_a<5; 2) typically, larger fluctuations in measured phase and amplitude were found in these experiments, and this resulted in larger errors. Using MSH, errors in µ_a and µ_s’ for phantoms 6 and 9 were larger than 10%, and this can be explained by reason #2 mentioned above. However, we still observe that errors in µ_s’ for phantom 3 and 4 were slightly larger than 10%.

 

Fig. 7. Percent errors in estimated OPs using MSH model for nine phantoms at locations 1, 2, and 3 that are shown by sub-Figs. (a), (b) and (c) respectively. Phantom identification numbers can be found in Table 5.

Download Full Size | PDF

 

Fig. 8. Percent errors in estimated OPs using ISM/SDA [sub-Fig. (a)] and MSH [sub-Fig. (b)] models for nine phantoms at location (5, 10)mm. Phantom identification numbers can be found in Table 5.

Download Full Size | PDF

Two points are worthy of note: 1) the relative error in µ_s’ is always greater than that in µ_a, and 2) their errors are higher than those found using MC generated data where input noise of about 10% was comparable to experimental noise [9]. It is helpful to understand which can be thought of as random errors and which are not. The OPs can be expressed in terms of modulation frequency, light speed in tissue, probe locations relative to the source and measured amplitudes and phase difference. In our two-distance detection scheme, it is generally assumed that the errors in modulation frequency and light speed are negligible. Hence, errors originate mainly from the measurements of distances, phase and amplitudes. Once two locations are chosen, they are fixed during the measurements of phase and amplitudes. Thus, errors in phase and amplitudes can be regarded as random, and the standard deviation can be calculated using several repeated measurements. However, errors due to distances should be regarded as systematic error in our scheme. This is different from the work done by most investigators, where multi-distance schemes were used and the effect on estimated OPs due to the inaccuracy in measured distances may be “averaged out” [29,30]. We can use the diffusion approximation (simple yet qualitatively correct) to calculate the random error in recovered µ_a and µ_s’ due to measured phase difference and fluences (error due to distance measurements will be treated separately and shown below). At the low frequency limit (ω/2π<150MHz), µ_a can be rewritten as follows:

where θ, A, and r are the phase difference, fluence, and distance from the source respectively; subscript denotes the position numbers; ω is the angular modulation frequency and v is light speed in the phantom. In accordance with the law of random error propagation, the random error in µ_a is:

Similarly, µ_s’ is given by:

where d is the distance between the two detector positions. The random error in µ_s’ is:

The fact that Δµ′_s /µ′_s is always larger than Δµ_a/µ_a is obvious from Eq. (11). Although a translation stage with 0.02mm precision was used in measuring the displacement, there are still factors which could potentially contribute to the errors, such as the difficulty of determining the absolute values of the distance between the probe and the source. Error due to these distance measurements cannot be calculated using the above formula because this error is not random. Tables 3 and 4 demonstrate that estimates of µ_s’ are affected more by systematic errors than estimates of µ_a.

It can be seen that both random and non-random errors contribute more to µ_s’ than µ_a. This is found to be generally true for determining OPs in infinite media using fluence based detection schemes in the steady state or in the frequency domain. Other sources of systematic errors include the distortion of photon propagation in the presence of the detector tip and the fiber which is connected to it, and non-uniform acceptance of light by the supposedly isotropic probe. In high absorption media, the benchmark values of OPs for India ink, (Eq. 3) and Intralipid, (Eq. 6) may become less accurate because large aggregates of ink particles may form as depicted in the paper by Madsen et al [20]. This may potentially contribute to systematic errors.

5. Discussion and conclusion

The work presented here is intended as a tool to predict the biological effect induced by PDT treatment. By measuring the OPs of the target tissue we are able to calculate light fluence and drug concentration. Knowledge of light fluence and photosensitizer concentration (called explicit dosimetry) is one possible method to predict the biological damage during PDT. Other alternatives exist, such as measuring phosphorescence of the photosensitizer or directly measuring the amount of singlet oxygen using the oxygen luminescence [31], which we will not discuss here.

The drug, Tookad is currently under investigation for treating prostate cancer and shows promising results [32, 33]. It appears that there is very limited information about the in vivo values of absorption coefficient of Tookad. Weersink et al [34, 35] used a steady state diffuse reflectance method applied on skins of rabbits and pigs, and determined the absorption at its peak value could be as high as 0.2mm^-1 for 5mg/kg injection dose of WST11 (a hydrophilic form of Tookad). Considering the reduced scattering of prostate can be as low as 0.3 mm^-1 [6] in canine prostate (although the wavelengths in two cases are slightly different), diffusion theory cannot give satisfactory results. Our work has shown that ISM/SDA or MSH can successfully address this situation. It is well known that PDT effects are reduced in hypoxia, and that hypoxia can be caused by high light fluence rate that leads to rapid consumption of oxygen [36]. If severe oxygen-depletion occurs during the treatment, a dose metric that uses only light fluence and photosensitizer concentration will not be a satisfactory indicator of true biological damage. Hence monitoring pO₂ during the treatment is desired. Many methods, such as Eppendorf probes or phosphorescence or fluorescence of exogenous probe molecules can be used to monitor partial oxygen tension pO₂. If we intend to add no more instrumentation to the current system, it may be possible to use the absorption spectrum to determine hemoglobin saturation and pO₂ [37]. Such an instrument would require a radio-frequency modulated source at several wavelengths, or a combination of frequency domain and continuous wave measurements as suggested by Bevilacqua et al [38] for reflectance.

In conclusion, our phantom experiments showed that interstitial measurement at two positions in the frequency domain can recover absorption within 10% and reduced scattering to slightly over 10% using either the data collected from two detectors 3mm apart or using those from two detectors 5mm apart with light collected by one detector attenuated by a neutral density filter. This was true for a wide range of OPs that can possibly be encountered in clinical PDT of the prostate, µ_a=0.008 to 0.17mm^-1 and µ_s’=0.3 to 1.8mm^-1. The experimental results confirmed the theoretical prediction that ISM/SDA can recover OPs within 5–10% for µ_s’/µ_a>5. The MSH method can recover OPs within 10% for µ_s’/µ_a>1, but generally greater than 10% if µ_s’/µ_a<1. We believe this to be caused by the large relative uncertainty in measured phase in high absorption media. Monte Carlo simulation showed that detectors should be situated at least 3mm away from the source in order to avoid the effects of source anisotropy. The simultaneous shift of two detectors either toward or away from the source by an amount S is less important than the situation where one detector is in the right position and the other is off by the same amount S. This suggests that accurate measurement for the distance between two detectors is required. The optimum strategy is probably to design a probe using 2 detector fibers in a fixed geometry, or even to fix all 3 optical fibers (one source, two detectors) in a single probe. Phase amplitude crosstalk, according to our findings, was practically insignificant. It is anticipated that the method described in this article will be used in future PDT dosimetry.