1 August 2011 New Monte Carlo model of cylindrical diffusing fibers illustrates axially heterogeneous fluorescence detection: simulation and experimental validation
Author Affiliations +
J. of Biomedical Optics, 16(8), 085003 (2011). doi:10.1117/1.3613920
We present a new Monte Carlo model of cylindrical diffusing fibers that is implemented with a graphics processing unit. Unlike previously published models that approximate the diffuser as a linear array of point sources, this model is based on the construction of these fibers. This allows for accurate determination of fluence distributions and modeling of fluorescence generation and collection. We demonstrate that our model generates fluence profiles similar to a linear array of point sources, but reveals axially heterogeneous fluorescence detection. With axially homogeneous excitation fluence, approximately 90% of detected fluorescence is collected by the proximal third of the diffuser for μs'/μa = 8 in the tissue and 70 to 88% is collected in this region for μs'/μa = 80. Increased fluorescence detection by the distal end of the diffuser relative to the center section is also demonstrated. Validation of these results was performed by creating phantoms consisting of layered fluorescent regions. Diffusers were inserted into these layered phantoms and fluorescence spectra were collected. Fits to these spectra show quantitative agreement between simulated fluorescence collection sensitivities and experimental results. These results will be applicable to the use of diffusers as detectors for dosimetry in interstitial photodynamic therapy.
Baran and Foster: New Monte Carlo model of cylindrical diffusing fibers illustrates axially heterogeneous fluorescence detection: simulation and experimental validation



Cylindrical diffusing fibers (diffusers) are optical fibers in which the distal portion emits light radially along some length. They are designed to provide homogeneous irradiance along the length of the diffusing section. This irradiance is also uniform in all radial directions. These diffusers are used in a number of areas in biomedicine, such as interstitial photocoagulation and interstitial hyperthermia.1

Perhaps the most prominent use of cylindrical diffusing fibers is in interstitial photodynamic therapy (iPDT). Photodynamic therapy (PDT) is an emerging treatment that has been used to treat malignant, premalignant, and benign conditions.2, 3 iPDT refers to the use of PDT for tumors that are not accessible to off-surface or intraluminal illumination. This is generally accomplished with the use of embedded cylindrical diffusers, although linear arrays of light emitting diodes are also used.4 iPDT has been used in humans to treat many types of cancer, including obstructing esophageal5 and bronchial cancer,6 prostate,7, 8 cholangiocarcinoma,9 glioma,10 and large tumors of the head and neck.11, 12 The outcome of iPDT is largely determined by the fluence distribution of the treatment light, the distribution and concentration of a photosensitizer, and the availability of molecular oxygen. Therefore, accurate determination of the fluence distribution from cylindrical diffusers is important for the prediction of the outcome and treatment planning.

A key part of determining the fluence distribution involves knowing the optical properties of the regions being treated. These optical properties depend on tissue properties, such as blood volume, and oxygenation, as well as photosensitizer absorption. The values of optical properties are typically determined in a spectroscopic fashion. Most photosensitizers are naturally fluorescent, so fluorescence measurements are often used.13 In iPDT, this typically involves the insertion of dedicated spectroscopy fibers.14 This increases clinical complexity and can lead to increased bleeding. In some cases, it would therefore be convenient to use the previously inserted treatment diffusers for spectroscopic measurements, especially fluorescence measurements. This is especially true for situations in which only a single fiber can be inserted into the desired volume, such as in the treatment of cholangiocarcinoma.9

In order to examine these problems, we would like to simulate light propagation in and around individual cylindrical diffusers, as well as arrays of diffusers. A common method for the simulation of light propagation in biological tissue is the Monte Carlo (MC) simulation. This technique treats light as a series of photons or photon packets which undergo absorption, scattering, and transmission and reflection at boundaries. Each of these phenomena is handled in a statistical fashion, based on the optical properties of the sample. If enough photons are run, MC simulation provides very accurate results. The accuracy of MC simulations has been demonstrated many times,15, 16 and MC is treated as the “gold standard” for simulation of light propagation in biological tissue. A number of MC codes are freely available, including Monte Carlo modeling of light transport in multilayered tissues (MCML) by Wang 17

In order to use cylindrical diffusers in the MC simulation, a model of them must be created. Previous efforts have treated cylindrical diffusers as linear arrays of ideal point sources.18, 19 While this makes computation more straightforward than in a more complex model, this description may not be fully accurate. Once photon packets are launched from the linear array of point sources, the effect of the diffuser on the propagation of light is neglected. The sources are treated only as launch points for photon packets, and the effect of the index mismatch between the diffuser and the surrounding tissue is ignored. Photon packets that would re-enter the diffuser are treated as still being in the surrounding tissue, which eliminates the possibility of using linear arrays of point sources to model detectors for spectroscopy. In order to model the detection of photon packets, these packets must be followed through the diffuser to the detection point. Therefore, a more complete MC model of cylindrical diffusing fibers is desired.

A disadvantage of conventional, CPU-based MC simulation is the long run-times of the simulations. A MC simulation of a typical tissue volume with an adequate number of photons can sometimes take hours or days.20 This clearly limits the use of MC in clinical practice. This has caused many investigators to use various approximations to represent light propagation in tissue, particularly the diffusion approximation.21 While this approximation can be calculated more rapidly, it is less accurate than MC simulation and places certain limitations on the measurement geometry.

One way to reduce the run-time of a MC simulation is through the use of a graphics processing unit (GPU). GPUs were originally designed to be used in computer graphics, and are therefore very efficient at performing calculations in parallel. Since each photon run is independent of other runs, MC simulation is an ideal candidate for parallelization on a GPU. This allows reduction of run-times by multiple orders of magnitude with no loss in accuracy. A number of GPU-based MC codes exist, capable of representing layered samples,22 three-dimensional (3D) samples with voxels,23 and three-dimensional samples with a triangular mesh.24 However, none of these code packages includes a diffuser model or tracks fluorescence generation and detection.

In this study, we have created a GPU-accelerated MC model of cylindrical diffusing fibers that is based on the physical construction of these fibers. This model has been used to simulate the irradiance delivered by these fibers, as well as the generation and detection of fluorescence by them. The simulation of fluorescence detection was particularly interesting, as heterogeneous detection along the length of the fiber was predicted. This heterogeneous detection sensitivity was experimentally validated using a layered phantom, showing quantitative agreement with simulation.




Model of Cylindrical Diffusing Fiber

Our diffuser model is based on a patent by E. L. Sinofsky,25 whose construction is shown schematically in Fig. 1. It consists of four components: a fiber core, a diffusive medium, a dielectric reflector, and cladding. The fiber core is the 400 μm core of a 0.22 numerical aperture (NA) jacketed optical fiber. The diffusive medium consists of titanium dioxide (TiO2) particles, with an average radius of 14.7 nm and an index of refraction of 2.488, embedded in a clear epoxy (Mastersil 151, Masterbond Inc., Hackensack, New Jersey) that has a density of 1 g/cm3 and a refractive index of 1.43.25 The number density of TiO2 particles providing a uniform axial fluence distribution is determined through simulation, as described in Sec. 3.1. The reflector is a dielectric stack, which is designed to provide maximum reflection with minimal heating. The cladding is a 0.25 mm layer of Teflon® FEP (DuPont, Wilmington, Delaware) with an index of 1.344 that surrounds the diffusive medium and dielectric reflector.

Fig. 1

Schematic of the cylindrical diffusing fiber, showing the four regions that are modeled. Photon packets are launched from the fiber core, with the diffusive medium being treated as having negligible absorption, and the cladding being treated as having negligible absorption and scattering. The dielectric reflector is treated as a perfect reflector.


Each of these four regions is modeled differently. The fiber core is treated as the launch point for photon packets, with packets being launched within the NA of the fiber. The face of the fiber core also serves as the detection point for fluorescence photon packets. Any fluorescence packet striking the fiber core within the NA of the fiber is recorded as being detected.

Since the diffusive medium consists only of TiO2 spheres and a clear epoxy, it is treated as having negligible absorption. Photon propagation within the diffuser is therefore governed only by scattering. The bulk scattering coefficient (μ s) within the diffusive medium depends on the concentration and size of the TiO2 particles and is assumed to be uniform. This scattering coefficient must be optimized to achieve axially homogeneous irradiance along the diffuser length. In order to determine the optimal scattering coefficient for each diffuser length, μ s was varied over a series of simulations and the homogeneity of the axial fluence profile was examined. The optimal value of μ s was found by minimizing the deviation from an axially homogeneous fluence profile along the length of the diffuser.

The cladding layer consists of clear plastic and does not contain any scatterers. It is therefore treated as having negligible absorption and scattering. Therefore, this layer has the effect of inducing Fresnel reflection and transmission. The dielectric reflector is treated as a perfect reflector.


GPU-Accelerated Monte Carlo Model

The propagation of photons is based upon the variance reduction MCML code by Wang 17 The movement, absorption, and scattering of photon packets are handled in the same way as in MCML. The code has been expanded to include the above diffuser model, a full 3D volume of optical properties, and fluorescence generation and detection.

Diffusers are specified by the three-dimensional coordinates, length, and radius of the diffusive region, and the outside radius of the cladding. Also specified are the NA of the fiber core, the scattering coefficient of the diffusive region, and the refractive indices of the diffusive region, and cladding. The total fluence to be delivered by the diffuser is specified in units of J/cm. An arbitrary number of these diffusers can be included in the simulation, with each diffuser able to launch photons and detect fluorescence. In the current version of the code, these diffusers are constrained to be parallel to the z-axis for simplicity. In clinical practice, diffusers are generally inserted in a parallel fashion,26 so this limitation is not a major one.

For photon propagation, each photon packet is launched at a random location on the specified fiber core at a random angle that is within the NA of the fiber. In the diffusive region, the photon packet propagates by normal MC methods. This includes reflections off of the dielectric reflector and cladding. Since the scatterers present in the diffusive region are much smaller than a wavelength, the scattering can be treated as isotropic (g = 0), assuming random orientations of the scatterers and random polarization of the light. Upon transmission into the cladding, propagation continues without scattering or absorption. If the photon is transmitted into the surrounding tissue, propagation is continued as normal with the specified scattering and absorption in the tissue region.

Boundary interaction with the diffuser is done by monitoring the photon position with respect to the coordinates of the diffuser axis. If the current step size of the photon would take it either into or out of a diffuser, the coordinate system is temporarily converted to cylindrical coordinates. A random number is then compared to the Fresnel reflection coefficient for the boundary. If this random number is less than the reflection coefficient, the photon is reflected and propagation continues. Otherwise, the photon undergoes refraction and then continues propagation.

In order to represent the 3D tissue volume in which the diffuser is embedded, the simulation space is sub-divided into a matrix of cuboid voxels. Each of these voxels stores the absorbed weight at that location, scaled by the user-specified fluence and the size of the voxel, as well as an identifying region number. Each of these region numbers is tied to a set of optical properties that is specified by the user. These optical properties consist of the absorption coefficient (μ a), the scattering coefficient (μ s), the index of refraction, the scattering anisotropy (g), and the quantum yield of fluorescence if that region is fluorescent. These properties are specified at both the excitation and fluorescence emission wavelengths.

The 3D nature of the simulation also means that boundary reflections must be handled in three-dimensions. An excellent overview of this is given in Pfefer 27 In our model, this is done by first determining all possible voxel boundaries that the photon packet will cross. The distance to each boundary is determined based on the current photon trajectory, and the region beyond that boundary is also recorded. The step size is then set to the shortest boundary distance that also leads into a different region. Reflection or transmission is determined based upon Fresnel coefficients, as described above, and photon propagation then continues. If the photon packet would not strike a boundary before executing its full step size, it is allowed to propagate its full step size.

During propagation of photon packets, fluorescence can be generated in the tissue. Whenever photon weight is absorbed, a fluorescence photon packet with a weight equal to the absorbed weight multiplied by the quantum yield of fluorescence is launched at that location. This fluorescence photon packet is immediately propagated by the above MC methods, but uses optical properties specified for the fluorescence emission wavelength. If a fluorescence photon packet strikes the fiber core inside the diffuser within the NA of the fiber, the photon is scored as detected. The position of its creation is recorded, as is the detected weight. This allows for the generation of a map of the origins of detected fluorescence photons, as well as what diffuser they were detected by in the case of an array of diffusers. This provides a description of regions that are sampled by spectroscopic measurements performed using the diffuser as a detector.

All of the above functions are performed on a GPU (GeForce GTX 570, NVIDIA Corporation, Santa Clara, California). It was programmed using NVIDIA's Compute Unified Device Architecture (CUDA) extensions to C. This allows for the code to run much faster than on a CPU. The GPU version of our model runs approximately 30 times faster than a CPU version of the same code. The parallel random number generator used in the code is based on that by Alerstam 22 This is required to ensure that photon packets being run simultaneously on separate threads use independent random numbers. The simulation parameters are stored in shared memory for faster access, while the absorption and detected fluorescence matrices are stored in the global device memory due to their large size. Atomic operations are required to update the absorption or detected fluorescence matrices in order to eliminate race conditions in memory access.


Phantom Preparation

In order to validate the MC results, multilayered phantoms were created. These phantoms consisted of either two or three layers, as shown in Fig. 2. In both cases, the bottom layer was a solid layer comprised of agar. 500 mg of agar (SELECT Agar®, Invitrogen Corporation, Carlsbad, California) was added to 50 ml of de-ionized water and heated to 95 °C while stirring. Heat was then removed and the solution was stirred until the temperature fell to 80 °C. At this point, 5 mL of 10% Liposyn (Liposyn® II 10%, Abbott Laboratories, Abbott Park, Illinois), 14.6 μl of India ink (Higgins No. 4418, Chartpak Incorporated, Leeds, Massachusetts), and 36.6 μl of 1 mg/ml fluorescein (Sigma-Aldrich, St. Louis, Missouri) were added. Stirring continued until the solution reached 60 °C, at which point it was refrigerated at 5 °C overnight to solidify. The addition of this concentration of Liposyn provides a scattering coefficient of ∼90 cm−1 at 488 nm.28 Assuming a scattering anisotropy of 0.82, this produces a reduced scattering coefficient of 16.2 cm−1. To ensure consistency, all phantoms were prepared from the same bottle of Liposyn. The addition of India ink brings the absorption coefficient to 2 cm−1 at 488 nm. The amount of fluorescein added gives a fluorophore concentration of 2 μM. This will be referred to as the fluorescein layer.

Fig. 2

Experimental set-up for measurement of fluorescence in (a) two-layer and (b) three-layer phantoms. The diffuser was inserted such that each layer bordered on an equal length of the diffuser. Only one diffuser is shown for clarity, but experiments used two diffusers inserted in parallel with a separation of 1 cm. One diffuser was used for delivery of an axially homogeneous fluorescence excitation profile, while the other was used for detection of fluorescence.


For the two layer phantom shown in Fig. 2a, the top layer was a liquid layer. This layer consisted of 10 ml of octanol (1-octanol, Sigma-Aldrich), 1 ml of 10% Liposyn, 2.8 μl of India ink, and 20.5 μl of 1 mg/ml 1,1-Dioctadecyl-3,3,3,3-tetramethylindocarbocyanine perchlorate (DiI) (Sigma-Aldrich). This was again designed to have a reduced scattering coefficient of 16.2 cm−1, an absorption coefficient of 2 cm−1, and a fluorophore concentration of 2 μM. This will be referred to as the DiI layer.

For the three layer phantom shown in Fig. 2b, the top layer was the DiI layer and the bottom layer was the fluorescein layer, both of which are described above. The middle layer consisted of 10 ml of 0.9% saline, 1 ml of 10% Liposyn, 3.06 μl of India ink, and 110.4 μl of 201.3 μM 2-[1-hexyloxyethyl]-2-devinyl pyropheophorbide-a (HPPH) dissolved in 1% Tween/phosphate buffered saline (PBS). HPPH is a second generation photosensitizer29 that is also fluorescent. It was prepared for us by Ravindra Pandey at the Roswell Park Cancer Institute in Buffalo, New York, using the procedure described in Pandey 30 This will be referred to as the HPPH layer.

Before use, each of the two liquid layers was stirred vigorously and poured on top of the fluorescein layer. Due to the immiscibility of octanol and saline solution, the DiI and HPPH layers remained distinct. HPPH is hydrophobic, but we dissolved it in Tween/PBS before diluting this solution in 0.9% saline. This keeps the HPPH confined to the saline solution, and we have seen no evidence of HPPH crossing into the DiI layer. Both liquid layers also remained distinct from the solid fluorescein layer.


Experimental Validation of Simulated Fluorescence Detection by Diffusers

For fluorescence measurements, two diffusers were inserted parallel into the above phantoms, with a separation of 1 cm. For two layer experiments, two 1-cm diffusers (Pioneer Optics Company, Bloomfield, Connecticut) were inserted 0.5 cm into the fluorescein layer shown in Fig. 2a. The DiI layer was designed to have a depth of 0.5 cm when poured into a 150 ml flask, which meant that the proximal half of the diffuser was fully in the DiI layer and the distal half was fully in the fluorescein layer. Therefore, any detected DiI fluorescence is known to have been generated along the proximal half of the diffuser, and any detected fluorescein fluorescence is known to have been generated along the distal half of the diffuser.

For three layer experiments, two 1.5-cm diffusers (PhotoGlow Inc., South Yarmouth, Massachusetts) were inserted 0.5 cm into the fluorescein layer shown in Fig. 2b. Both the DiI and HPPH layers were designed to have a depth of 0.5 cm when poured into a 150 ml flask, which meant that the diffuser length was divided into three equal sections.

In both cases, excitation light at 488 nm was provided by an argon-ion laser (Innova 70, Coherent, Santa Clara, California), filtered by a band-pass filter (Z488Trans-pc-xr, Chroma Technology Corp, Bellows Falls, Vermont). Excitation light was delivered at 40 mW/cm for two-layer experiments and at 15 mW/cm for three-layer experiments. The optical power was lower for three-layer experiments due to larger coupling losses in the 1.5 cm diffusers. Optical power was measured using an integrating-sphere-based laser power measurement system (LPMS-060-SF-Si, Labsphere, North Sutton, New Hampshire). Excitation light was delivered through one of the diffusers in the experiment, and fluorescence was detected by the other. This was done to ensure homogeneous excitation along the diffuser length, and to minimize the detection of autofluorescence generated within the fiber. Light captured by the second diffuser was filtered through a long-pass filter (HQ500LP, Chroma) before being detected by a TE-cooled, 16 bit spectrometer (Compass X, B&W Tek, Newark, Delaware) using an integration time of 20 s.

For calibration purposes, fluorescence spectra were also collected from phantoms consisting of only one of the fluorophore layers. In order to do this, larger volumes of the layers described above were created while maintaining the concentrations of each component. This allowed the diffusers to be fully submerged in each layer.


Data Correction and Spectral Fitting

All fluorescence spectra collected were corrected for background and system response. This was done by subtracting a dark measurement from the raw fluorescence measurement and dividing by a wavelength-dependent system response. Dark measurements were acquired by integrating dark signals for the durations described above. System responses were acquired by placing a single diffuser in the center of an integrating sphere (3P-LPM-060-SF, Labsphere) and shining a NIST-traceable lamp (Model No. LS-1-CAL, Ocean Optics, Dunedin, Florida) through one of the detection ports on the sphere. The measured spectrum was then background subtracted and divided by the known lamp spectrum to obtain the system response.

After fluorescence spectra were corrected, they were fit using a singular value decomposition (SVD) algorithm based on the work of Press 31 Basis spectra for each of the fluorophores were used in the fits. These basis spectra were acquired by measuring the fluorescence emission of each layer in a commercial fluorometer (Varian Eclipse, Palo Alto, California). In addition to the basis spectra, a series of 61 Fourier terms was also used in fitting to account for unknown possible contributions to the measured fluorescence. These Fourier terms were given a smaller weight in the fitting in order to favor the basis spectra. Fitting followed the SVD scheme found in MATLAB (Mathworks, Natick, Massachusetts).

After SVD fitting, the fit magnitudes for each fluorophore were further corrected for absorption and fluorescence yield effects. Each of the three fluorophores used has some absorption at 488 nm, and emits fluorescence beyond 500 nm. However, each fluorophore has different absorption at 488 nm and a different quantum yield of fluorescence. As mentioned previously, fluorescence spectra were acquired for each layer individually. These “pure” measurements were made using the same excitation wavelength and power as in the layered measurements, and with the same concentrations of scatterer, absorber, and fluorophore. So, any differences in absorption at 488 nm or fluorescence quantum yield were also present in the pure measurements. Therefore, any difference in magnitude between the pure and layered measurements was due to the detection capability of the diffuser at the height of the appropriate layer. To account for this, we divide the fit magnitude of the layered measurement by that of the pure measurement. This gives a relative measurement of the detected fluorescence that is weighted for absorption and fluorescence quantum yield.




Determination of the Scattering Coefficient Within the Diffusive Region

As mentioned previously, the scattering coefficient within the diffuser determines the homogeneity of the axial fluence profile at the surface of the diffuser. If the value of μ s is too large, most of the light launched from the fiber face will scatter out of the diffuser near the proximal end. If the value is too small, then the scattering out of the diffuser will be weak, which decreases the efficiency of the device. Figure 3 shows the simulated degradation in axial fluence profile for a 1 cm diffuser that is induced by use of an improper scattering coefficient. Figure 3a shows the homogeneous fluence profile along the surface of a 1 cm diffuser for the determined optimal value of μ s = 0.2009 cm−1 inside the diffuser. Figure 3b shows the effect of raising the value of μ s to 10 cm−1 within the diffuser, which results in too much light scattering out of the diffuser near the proximal end. The value of μ s within the diffuser is therefore crucial to the generation of a homogeneous axial fluence profile.

Fig. 3

Simulated irradiance profiles along the surface of a 1 cm diffuser, with its proximal end at z = 1 cm, illustrating the effect of changing μs. Simulation parameters are identical in both cases except for the value of μ s inside the diffusive region, which was (a) 0.2009 cm−1 and (b) 10 cm−1.


The optimal μ s for homogeneous axial fluence changes with diffuser length. The optimal μ s values computed using our MC model for diffuser lengths between 1 and 5 cm are shown in Table 1. As can be seen, the required μ s for homogeneous fluence decreases with increasing diffuser length. This corresponds to a reduced concentration of scatterers through μ s = σN, where σ is the scattering cross-section and N is the number density of scatterers.32 The scattering cross-section can be calculated using


[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \sigma = \frac{{8\pi }}{3}k^4 a^6 \left| {\frac{{n_s ^2 - n_m ^2 }}{{n_s ^2 + 2n_m ^2 }}} \right|, \end{equation}\end{document} σ=8π3k4a6ns2nm2ns2+2nm2,
which is valid for the scatterer size and wavelength we are considering.32 Here, k is the wave vector, a is the scatterer radius, n s is the index of the scatterer, and n m is the index of the surrounding medium. The scatterer concentrations in Table 1 were calculated using a = 14.7 nm, n s = 2.488, and n m = 1.43. The values calculated for 1 and 5 cm diffusers match those used in commercial diffusers manufactured by Pioneer Optics Company.33 Number densities of scatterers were converted to parts per million (ppm) by weight using the known density and volume of the diffusive region and the molecular weight of TiO2.

Fig. 9

Comparison between simulated and experimental fluorescence detection using (a) 1 cm and (b) 1.5 cm diffusers. Heights of experimental bars (□) indicate the mean value of (a) n = 6 and (b) n = 4 experiments, with error bars representing standard deviation. Values used were corrected for background, system response, fluorescence quantum yield, and absorption at 488 nm. Heights of simulated bars (■) indicate the mean value of 3 simulations, with error bars representing the standard deviation. The value of μ a was set to 2 cm−1 for both simulation and experiment.


Table 1

Optimal scattering coefficient (μ s) for homogeneous irradiance for multiple diffuser lengths, and the corresponding calculated number density of scatterers (N). N is given in units of cm−3 and ppm by weight. Number densities in ppm by weight used in commercial diffusers manufactured by Pioneer Optics Company are shown for comparison.

DiffuserμsNCalculatedPioneer diffuser
length (cm)(cm−1)(cm−3)N (ppm)N (ppm)


Fluence Distributions Generated by Monte Carlo Model of Diffuser

A typical fluence map generated by a 1 cm diffuser is shown in Fig. 4. The air-tissue boundary is located at z = 0, and the proximal end of the diffuser is at z = 0.5 cm. For this particular simulation, the tissue sample is homogeneous and semi-infinite with μ a = 2 cm−1, μ s = 100 cm−1, g = 0.9, and n = 1.395 (Ref. 34).

Fig. 4

3D rendering of the fluence distribution around a 1 cm cylindrical diffusing fiber with its proximal end at z = 0.5 cm and an air-tissue boundary at z = 0. Tissue optical properties were set to μ a = 2 cm−1, μ s = 100 cm−1, and g = 0.9. Voxel size was 0.02×0.02×0.02 cm.


In order to compare to a model based on a finite linear array of point sources, the radial degradation of fluence from the diffuser was examined for both models. In these simulations, we assumed a refractive index of 1.43 for the linear array to match that of the epoxy in the diffuser. Cuts through the fluence distribution were made radially from the axial mid-point of the diffuser for both models. The results of this are shown in Fig. 5. As can be seen, the fluence cuts are virtually identical for the two models. Similar results were obtained for cuts at other axial positions on the diffuser, and at other angles (data not shown). This indicates that a model based on a linear array of point sources is sufficient for determination of fluence distribution.

Fig. 5

Comparison between linear array of point sources model and our MC model in terms of radial degradation of fluence for μ a = 2 cm−1, showing substantial overlap between the two methods. Shown is a cut through the fluence at the axial center of a 1 cm diffuser. Simulation parameters were identical, except for the source model. The arrow indicates the position of the outer radius of the diffuser.



Modeling of Fluorescence Generation and Detection

As mentioned previously, fluorescence photon packets are generated at the location of each absorption event within the tissue medium surrounding a diffuser. Generated fluorescence photon packet weights consist of the absorbed weight, scaled by the quantum yield of fluorescence. Therefore, the distribution of generated fluorescence is a scaled version of the absorption distribution. This is shown in Fig. 6a for a 1 cm diffuser. As expected, this generated fluorescence is homogeneous along the length of the diffuser. Figure 6b shows the origins and weights of generated fluorescence photons that cross into the diffuser. This is simply a scaled version of the generated fluorescence, and is again homogeneous along the length of the diffuser. The distribution of the origins of detected fluorescence photon packets is shown in Fig. 6c. Each pixel magnitude is the sum of the detected weight of fluorescence photon packets generated at that location that migrated through the diffuser, reached the fiber core, and were scored as detected. This distribution is highly heterogeneous, with a majority of the detected fluorescence being generated near the proximal end of the diffuser.

Fig. 6

(a) Fluorescence generated in tissue by a 1-cm diffuser with its proximal end at z = 2 cm and an outer radius of 0.05 cm. (b) Origins of fluorescence photons that crossed into the diffuser after being generated in the surrounding tissue. (c) Origins of fluorescence photons that were detected by the diffuser. μ a was set to 2 cm−1 in tissue. Planar cuts are shown through the simulated volume at the center of the diffuser. Only the right half of this plane is shown for clarity. The left half is identical.


Axial profiles of the simulated detected fluorescence are shown in Fig. 7. These profiles were created by summing the detected fluorescence at each z position within the simulated tissue volume. Optical properties of the tissue sample were set to μ s = 90 cm−1, g = 0.82, and n = 1.395, with either μ a = 2 cm−1 to match experimental conditions or μ a = 0.2 cm−1 to represent a typical value in tissue.35 As can be seen, the majority of detected fluorescence comes from the proximal portion of the diffuser, with a slightly increased contribution from the distal portion of the diffuser relative to the center. Detected fluorescence, detailed by the diffuser segment, is shown in Table 2. In all cases, the majority of detected fluorescence comes from the proximal third of the diffuser, with contributions ranging from 86 to 92% for μ a = 2 cm−1 and from 70 to 88% for μ a = 0.2 cm−1. In all but the 1 cm case for both values of μ a and the 1.5 cm case for μ a = 2 cm−1, the distal third of the diffuser is responsible for the next largest contribution to detected fluorescence, with the center segment having the smallest contribution. Relative fluorescence collection by the distal third ranges from 5.38 to 6.66% for μ a = 2 cm−1 and 9.11 to 11.98% for μ a = 0.2 cm−1, with the center section collecting 1.63 to 6.11% for μ a = 2 cm−1 and 2.03 to 18.15% for μ a = 0.2 cm−1.

Fig. 7

Simulated detected fluorescence by axial position along (a) 1 to 2 cm diffusers with μ a = 2 cm−1, (b) 3 to 5 cm diffusers with μ a = 2 cm−1, (c) 1 to 2 cm diffusers with μ a = 0.2 cm−1, and (d) 3 to 5 cm diffusers with μ a = 0.2 cm−1. All simulations used μ s = 90 cm−1 and g = 0.82, and placed the diffuser's proximal end at z = 0.5 cm. The arrows indicate the location of the proximal end of the diffusers.


Table 2

Percentage of generated fluorescence that is collected by various segments of different diffuser lengths, as determined by our MC model. Thirds correspond to the three-layer experiment, and halves correspond to the two-layer experiment. Uncertainties given are the standard deviation of three successive runs of the simulation.

DiffuserμaDiffuser segment (thirds)Diffuser segment (halves)
10.269.88 ± 0.0218.15 ± 0.0311.98 ± 0.0181.17 ± 0.0118.83 ± 0.01
2.086.22 ± 0.177.76 ± 0.086.02 ± 0.0991.54 ± 0.158.46 ± 0.15 ± 0.0812.76 ± 0.019.64 ± 0.0986.06 ± 0.0413.94 ± 0.04
2.090.48 ± 0.123.92 ± 0.075.60 ± 0.1293.01 ± 0.136.99 ± 0.13
20.282.54 ± 0.118.26 ± 0.089.21 ± 0.0288.10 ± 0.0611.90 ± 0.06
2.092.03 ± 0.042.58 ± 0.075.39 ± 0.0993.52 ± 0.076.48 ± 0.07
30.286.77 ± 0.044.08 ± 0.019.16 ± 0.0389.46 ± 0.0210.54 ± 0.02
2.092.71 ± 0.121.91 ± 0.085.38 ± 0.1393.75 ± 0.106.25 ± 0.10
40.288.49 ± 0.042.40 ± 0.079.11 ± 0.0389.89 ± 0.0110.11 ± 0.01
2.092.92 ± 0.041.57 ± 0.055.50 ± 0.0493.73 ± 0.036.27 ± 0.03
50.288.66 ± 0.922.03 ± 0.019.32 ± 0.1189.73 ± 0.1310.27 ± 0.13
2.093.34 ± 0.071.56 ± 0.065.10 ± 0.1394.10 ± 0.105.90 ± 0.10

Simulations were run with 1,000,000 photon packets on a 100×100×100 grid of voxels. Due to the statistical nature of MC modeling, simulations were repeated three times in order to evaluate reproducibility. The variation in detected fluorescence varies by at most 4% between runs, with most values varying by less than 1%. The run-times were approximately 1 min for simulations of fluence distribution and approximately 30 min for simulations of fluorescence detection.


Experimental Validation of Fluorescence Detection

Experiments were performed in two- and three-layer phantoms with 1 and 1.5 cm diffusers, respectively. Typical SVD fits to corrected fluorescence spectra are shown in Fig. 8. The two-layer fit in Fig. 8a shows contributions from DiI and fluorescein, while the three-layer fit in Fig. 8b shows contributions from DiI, fluorescein, and HPPH. In both cases, the magnitude of the Fourier terms is relatively small, indicating a good fit to the data.

Fig. 8

Results of SVD fitting to representative fluorescence spectra collected from (a) two-layer and (b) three-layer phantoms, with spectra corrected for background and system response. Fit magnitudes shown are not corrected for the effects of fluorescence quantum yield and absorption at 488 nm.


Fit coefficients for the fluorophores were further corrected for absorption and fluorescence quantum yield as described above. The corrected values were then used to determine the relative contributions of each segment of the diffuser to the overall detected fluorescence. The results of this are shown in Fig. 9, with the simulated values shown for comparison. For both two- and three-layer experiments, there is agreement between the simulation and experiment to better than 6%. The results shown are for a set of 6 two-layer experiments, and a set of 4 three-layer experiments. Error bars are standard deviations. Error bars for simulated data are the standard deviations from three successive runs of the simulation.



Simulations of both the linear array of point sources and our diffuser model show very similar results for radial fluence profiles, demonstrating that it is appropriate to use the linear array of point sources model when only the fluence distribution for a

single diffuser is desired. This less complex model will run slightly faster, since the propagation of photons within the diffuser and boundary collisions with the diffuser are not computed.

To the best of our knowledge, this is the first time that heterogeneous collection of fluorescence by a cylindrical diffuser has been demonstrated, either experimentally or by simulation. As shown here, most of the fluorescence collected by a diffuser is generated near its proximal end with a relatively small spike near the distal end. Given the long scattering mean free path within the diffuser, the large contribution from the proximal end may be explained simply in terms of the diffuser geometry. The NA of the fiber core and the index mismatch between the tissue, cladding, and diffusive medium define the range of angles over which fluorescence entering the diffuser is detected. The axial-coordinate of a photon incident on the diffuser and the index mismatch combine to constrain the range of angles over which the photon can enter the diffuser and be detected. As the axial-coordinate of the incoming photon moves further from the fiber core, this range of acceptable angles decreases. Because fluorescence in the tissue is emitted isotropically, photons are incident on the diffuser at a random angle. Therefore, if the range of acceptable angles decreases with increasing distance from the fiber core and the incident angles are random, the detection of photons also decreases with increasing distance from the fiber core. The spike in collection at the distal end can be explained by the presence of the dielectric reflector. Photons that enter the diffuser close to the distal end are more likely to strike the reflector, which directs the trajectory of a fraction of these toward the proximal end of the diffuser. Coupled with reflections off of the cladding, this can result in increased probability of detection. In geometric terms, this has the effect of increasing the range of incoming photon angles that can reach the fiber core.

It is important to note that the detection distribution of a cylindrical diffuser is sensitive to the optical properties of the tissue. In particular, the absorption of the tissue will change the axial detection profile. As we have shown, the proximal portion of the diffuser is most sensitive to the detection of fluorescence. With lower tissue absorption, the likelihood of a fluorescence photon reaching the proximal portion of the diffuser before being absorbed is increased. This allows fluorescence that is generated a larger axial distance from the fiber core to be detected. For example, a value of 2 cm−1 for μ a results in approximately 86% of detected fluorescence being generated along the proximal third of the diffuser. For μ a = 0.2 cm−1, this value is reduced to 70%, with a greater proportion of the detected fluorescence being generated further from the proximal end of the diffuser. While the fluorescence detection profile is still heterogeneous, the sharp degradations shown in Figs. 7a and 7b are smoothed out for lower values of μ a, as seen in Figs. 7c and 7d.

Knowledge of the detection distribution of a cylindrical diffuser is important if these fibers are to be used for dosimetry in iPDT. Given a detected fluorescence spectrum, spectroscopic methods can be used to estimate optical properties and photosensitizer concentrations. Having a MC model of the origins of detected fluorescence photons allows for these quantities to be mapped to a specific tissue volume. However, a disadvantage of using cylindrical diffusing fibers as spectroscopy probes is that an incomplete map of optical properties would be generated. As we have shown, the detection sensitivity of a cylindrical diffuser is highly heterogeneous, with certain regions near the diffuser not being sampled by spectroscopic measurements. This means that spectroscopic data would not be available for all locations within the tissue, requiring assumptions of optical properties in un-sampled regions. In schemes that use separate spectroscopy fibers, a full map of optical properties can be generated by translating the spectroscopy fibers.36

The use of separate spectroscopy fibers is not always an option, as in some cases, only a single fiber can physically be inserted into the desired volume. This is typically the situation in iPDT for cholangiocarcinoma,9 as well as in the treatment of other hollow organs. In these cases, the only way to obtain spectroscopic information about the tissue would be to use the single treatment diffuser as both a source and detector.

Of critical importance to the clinical deployment of this technique is the speed of the MC simulation. As previously noted, traditional MC simulations can take several hours to run. With our GPU-accelerated model, we have reduced this run-time to the order of 1 min. While this is a large improvement, it would be desirable to further reduce this to the order of 1s. Alerstam have shown a 300× speed-up by moving to a GPU-accelerated version of MCML.22 This would be excellent for our model, but it is unclear whether this speed-up is achievable for our scenario. The model presented by Alerstam only allows for layered samples and does not incorporate a diffuser model or fluorescence. This eliminates the need for boundary checks at each voxel and at the diffuser boundary, and removes the need for atomic operations in the detection of fluorescence. These differences make a large difference in run-time due to the nature of the GPU. Each of the boundary checks currently requires a substantial number of conditional statements. When executed on the GPU, all branches of conditionals are serially executed before converging after the conditional.37 This results in an increase in run-time. The use of atomic operations for fluorescence detection also results in a significantly increased run-time. Other studies have shown that the use of atomic operations in GPU-accelerated MC can result in a slow-down of 75%.23 These factors do not mean that further speed-up is impossible for our code, just that it will require significant removal of branching and intelligent usage of atomic operations.

Our MC model may also be useful in the design of cylindrical diffusing fibers. As mentioned previously, the output irradiance of the diffuser is largely determined by the scattering properties of the diffusive region. These scattering properties are determined by the size and index of the scatterer used, as well as the density of the scatterer. Once a scatterer is chosen, its size and index are known. By simulating a range of μ s values and optimizing homogeneity, the density can be computed. This can then be used to manufacture the desired diffuser.

Alternatively, some other nonuniform irradiance profile along the diffuser may be desired. Our model can again be used to predict the scatterer concentration required to achieve this profile. This could be done by adjusting a uniform scatterer concentration, or by allowing a heterogeneous distribution of scatterers. Having a nonuniform distribution of scatterers is often used in the design of two-dimensional diffusers to generate a desired illumination pattern,38 and could therefore be applied to the design of cylindrical diffusers.

We have demonstrated a new MC model of cylindrical diffusing fibers. This model predicts a fluence distribution similar to that of a linear array of point sources and a heterogeneous fluorescence detection distribution. The heterogeneous detection of fluorescence has been experimentally demonstrated, and indicates agreement with the model.


This work was supported by NIH Grant Nos. CA68409 and CA122093 awarded by the National Cancer Institute. The authors wish to thank Ron Hille from Pioneer Optics Company for providing parameters used in commercial diffusers.



L. M. Vesselov, W. Whittington, and L. Lilge, “Performance evaluation of cylindrical fiber optic light diffusers for biomedical applications,” Lasers Surg. Med. 34, 348–351 (2004).10.1002/lsm.20031Google Scholar


P. Agostinis, K. Berg, K. A. Cengel, T. H. Foster, A. W. Girotti, S. O. Gollnick, S. M. Hahn, M. R. Hamblin, A. Juzeniene, D. Kessel, M. Korbelik, J. Moan, P. Mroz, D. Nowis, J. Piette, B. C. Wilson, and J. Golab, “Photodynamic therapy of cancer: an update,” CA Cancer J. Clin. 61, 250–281 (2011).10.3322/caac.20114Google Scholar


A. E. O’Connor, W. M. Gallagher, and A. T. Byrne, “Porphyrin and nonporphyrin photosensitizers in oncology: preclinical and clinical advances in photodynamic therapy,” Photochem. Photobiol. 85, 1053–1074 (2009).10.1111/j.1751-1097.2009.00585.xGoogle Scholar


R. A. Lustig, T. J. Vogl, D. Fromm, R. Cuenca, R. A. Hsi, A. K. D’Cruz, Z. Krajina, M. Turic, A. Singhal, and J. C. Chen, “A multicenter phase I safety study of intratumoral photoactivation of talaporfin sodium in patients with refractory solid tumors,” Cancer 98(8), 1767–1771 (2003).10.1002/cncr.11708Google Scholar


C. J. Lightdale, S. K. Heier, N. E. Marcon, J. James, S. McCaughan, H. Gerdes, B. F. Overholt, J. Michael, V. Sivak, G. V. Stiegmann, and H. R. Nava, “Photodynamic therapy with porfimer sodium versus thermal ablation therapy with Nd:YAG laser for palliation of esophageal cancer: a multicenter randomized trial,” Gastrointest. Endosc. 42(6), 507–512 (1995).10.1016/S0016-5107(95)70002-1Google Scholar


K. Moghissi, K. Dixon, M. Stringer, T. Freeman, A. Thorpe, and S. Brown, “The place of bronchoscopic photodynamic therapy in advanced unresectable lung cancer: experience of 100 cases,” Eur. J. Cardiothoracic Surg. 15, 1–6 (1999).10.1016/S1010-7940(98)00295-4Google Scholar


H. Patel, R. Mick, J. Finlay, T. C. Zhu, E. Rickter, K. A. Cengel, S. B. Malkowicz, S. M. Hahn, and T. M. Busch, “Motexafin lutetium-photodynamic therapy of prostate cancer: short- and long-term effects on prostate-specific antigen,” Clin. Cancer Res. 14(15), 4869–4876 (2008).10.1158/1078-0432.CCR-08-0317Google Scholar


C. M. Moore, T. R. Nathan, W. R. Lees, C. A. Mosse, A. Freeman, M. Emberton, and S. G. Brown, “Photodynamic therapy using meso tetra hydroxy phenyl chlorin (mTHPC) in early prostate cancer,” Lasers Surg. Med. 38, 356–363 (2006).10.1002/lsm.20275Google Scholar


M. E. J. Ortner, K. Caca, F. Berr, J. Liebetruth, U. Mansmann, D. Huster, W. Voderholzer, G. Schachschal, J. Mossner, and H. Lochs, “Successful photodynamic therapy for nonresectable cholangiocarcinoma: a randomized prospective study,” Gastroenterology 125(5), 1355–1363 (2003).10.1016/j.gastro.2003.07.015Google Scholar


T. J. Beck, F. W. Kreth, W. Beyer, J. H. Mehrkens, A. Obermeier, H. Stepp, W. Stummer, and R. Baumgartner, “Interstitial photodynamic therapy of nonresectable malignant glioma recurrences using 5-aminolevulinic acid induced protoporphyrin IX,” Lasers Surg. Med. 39, 386–393 (2007).10.1002/lsm.20507Google Scholar


P.-J. Lou, H. Jager, L. Jones, T. Theodossy, S. Brown, and C. Hopper, “Interstitial photodynamic therapy as salvage treatment for recurrent head and neck cancer,” Br. J. Cancer 91, 441–446 (2004).10.1038/sj.bjc.6601993Google Scholar


M. A. Biel, “Photodynamic therapy treatment of early oral and laryngeal cancers,” Photochem. Photobiol. 83, 1–6 (2007).10.1111/j.1751-1097.2007.00153.xGoogle Scholar


J. C. Finlay, T. C. Zhu, A. Dimofte, D. Stripp, S. B. Malkowicz, T. M. Busch, and S. M. Hahn, “Interstitial fluorescence spectroscopy in the human prostate duing motexafin lutetium-mediated photodynamic therapy,” Photochem. Photobiol. 82, 1270–1278 (2006).10.1562/2005-10-04-RA-711Google Scholar


S. R. H. Davidson, R. A. Weersink, M. A. Haider, M. R. Gertner, A. Bogaards, D. Giewercer, A. Scherz, M. D. Sherar, M. Elhilali, J. L. Chin, J. Trachtenberg, and B. C. Wilson, “Treatment planning and dose analysis for interstitial photodynamic therapy of prostate cancer,” Phys. Med. Biol. 54, 2293–2313 (2009).10.1088/0031-9155/54/8/003Google Scholar


L. H. P. Murrer, J. P. A. Marijnissen, and W. M. Star, “Ex vivo light dosimetry and Monte Carlo simulations for endobronchial photodynamic therapy,” Phys. Med. Biol. 40, 1807–1817 (1995).10.1088/0031-9155/40/11/003Google Scholar


Q. Liu, C. Zhu, and N. Ramanujam, “Experimental validation of Monte Carlo modeling of fluorescence in tissues in the UV-visible spectrum,” J. Biomed. Opt. 8(2), 223–236 (2003).10.1117/1.1559057Google Scholar


L. Wang, S. L. Jacques, and L. Zheng, “MCML - Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47, 131–146 (1995).10.1016/0169-2607(95)01640-FGoogle Scholar


L. H. Murrer, H. P. Marijnissen, and W. M. Star, “Monte Carlo simulations for endobronchial photodynamic therapy: the influence of variations in optical and geometrical properties and of realistic and eccentric light sources,” Lasers Surg. Med. 22, 193–206 (1998).10.1002/(SICI)1096-9101(1998)22:4<193::AID-LSM3>3.0.CO;2-KGoogle Scholar


B. Farina, S. Saponaro, E. Pignoli, S. Tomatis, and R. Marchesini, “Monte Carlo simulation of light fluence in tissue in a cylindrical diffusing fibre geometry,” Phys. Med. Biol 44, 1–11 (1999).10.1088/0031-9155/44/1/002Google Scholar


D. A. Boas, J. P. Culver, J. J. Stott, and A. K. Dunn, “Three-dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002).Google Scholar


M. S. Thompson, A. Johansson, T. Johansson, S. Andersson-Engels, S. Svanberg, N. Bendsoe, and K. Svanberg, “Clinical system for interstitial photodynamic therapy with combined on-line dosimetry measurements,” Appl. Opt. 44(19), 4023–4031 (2005).10.1364/AO.44.004023Google Scholar


E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing graphics processing units for high-speed Monte Carlo simulation of photon migration,” J.Biomed. Opt. 13(6), 060504 (2008).10.1117/1.3041496Google Scholar


Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178 (2009).10.1364/OE.17.020178Google Scholar


N. Ren, J. Liang, X. Qu, J. Li, B. Lu, and J. Tian, “GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues,” Opt. Express 18(7), 6811 (2010).10.1364/OE.18.006811Google Scholar


E. L. Sinofsky, “Phototherapeutic apparatus with diffusive tip assembly,” U.S. Patent No. 5,947,959 (1999).Google Scholar


M. D. Altschuler, T. C. Zhu, J. Li, and S. M. Hahn, “Optimized interstitial PDT prostate treatment planning with the Cimmino feasibility algorithm,” Med. Phys. 32(12), 3524–3536 (2005).10.1118/1.2107047Google Scholar


T. J. Pfefer, J. K. Barton, E. K. Chan, M. G. Ducros, B. S. Sorg, T. E. Milner, J. S. Nelson, and A. J. Welch, “A three-dimensional modular adaptable grid numerical model for light propagation during laser irradiation of skin tissue,” IEEE J. Sel. Top. Quantum Electron. 2(4), 934–942 (1996).10.1109/2944.577318Google Scholar


H. J. v. Staveren, C. J. M. Moes, J. v. Marle, S. A. Prahl, and M. J. C. v. Gemert, “Light scattering in Intralipid-10% in the wavelength range of 400-1100 nm,” Appl. Opt. 30(31), 4507–4514 (1991).10.1364/AO.30.004507Google Scholar


D. A. Bellnier, W. R. Greco, H. Nava, G. M. Loewen, A. R. Oseroff, and T. J. Dougherty, “Mild skin photosensitivity in cancer patients following injection of Photochlor (2-[1-hexyloxyethyl]-2-devinyl pyropheophorbide-a; HPPH) for photodynamic therapy,” Cancer Chemother. Pharmacol. 57, 40–45 (2005).10.1007/s00280-005-0015-6Google Scholar


R. K. Pandey, F.-Y. Shiau, A. B. Sumlin, T. J. Dougherty, and K. M. Smith, “Structure/activity relationships among photosensitizers related to pheophorbides and bacteriopheophorbides,” Bioorg. Med. Chem. Lett. 2(5), 491–496 (1992).10.1016/S0960-894X(00)80176-6Google Scholar


W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, New York (1992).Google Scholar


A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE Press, New York (1997).Google Scholar


R. Hille, personal communication (2010).Google Scholar


F. P. Bolin, L. E. Preuss, R. C. Taylor, and R. J. Ference, “Refractive index of some mammalian tissues using a fiber optic cladding method,” Appl. Opt. 28(12), 2297–2303 (1989).10.1364/AO.28.002297Google Scholar


E. Hull, D. Conover, and T. Foster, “Carbogen-induced changes in rat mammary tumor oxygenation reported by near infrared spectroscopy,” Br. J. Cancer 79, 1709–1716 (1999).10.1038/sj.bjc.6690272Google Scholar


T. C. Zhu, J. C. Finlay, and S. M. Hahn, “Determination of the distribution of light, optical properties, drug concentration, and tissue oxygenation in-vivo in human prostate during motexafin lutetium-mediated photodynamic therapy,” J. Photochem. Photobiol. B. 79, 231–241 (2005).10.1016/j.jphotobiol.2004.09.013Google Scholar


“NVIDIA CUDA C Programming Guide,” NVIDIA Corporation (2010).Google Scholar


A. Tagaya, M. Nagai, Y. Koike, and K. Yokoyama, “Thin liquid-crystal display backlight system with highly scattering optical transmission polymers,” Appl. Opt. 40(34), 6274–6280 (2001).10.1364/AO.40.006274Google Scholar

Timothy M. Baran, Thomas H. Foster, "New Monte Carlo model of cylindrical diffusing fibers illustrates axially heterogeneous fluorescence detection: simulation and experimental validation," Journal of Biomedical Optics 16(8), 085003 (1 August 2011). http://dx.doi.org/10.1117/1.3613920



Optical fibers

Monte Carlo methods

Fluorescence spectroscopy



Back to Top