## 1.

## Introduction

Fluorescence spectroscopy has been used successfully to discriminate premalignancy and malignancy in a number of organ sites.^{1} However, due to the complex interplay of absorption, scattering, and fluorescence in tissue, it is difficult to separate the intrinsic fluorescence properties from absorption and scattering, thus making these spectra difficult to interpret.

To address this issue, a number of groups have proposed methods to determine the fluorophore concentration or intrinsic fluorescence spectra (i.e., the fluorescence properties independent of absorption and scattering) from a measured fluorescence spectrum.
^{2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14} These approaches range from empirical methods
^{2, 5, 7, 11} to diffusion theory modeling.^{3, 4} These approaches are generally limited in that they are valid for a limited range of absorption and scattering, they require extensive empirical calibration, and/or they are not flexible in their applicability to different optical probe geometries. To address these concerns a Monte-Carlo-based model to extract intrinsic tissue fluorescence is presented here. This method has the advantage that it is based on Monte Carlo modeling, and so, has the theoretical advantage of not having to impose constraints on the range of optical properties being modeled. Furthermore, this method can model the actual fiber optic probe geometry used for the fluorescence measurement and thus, in principle, is flexible in its application to a variety of probe geometries. Moreover, the approach requires only a single phantom measurement to enable adaptation to any optical system configuration (e.g., different probe geometries). This paper demonstrates its application to fluorescence spectra measured from tissue phantoms with one particular probe geometry/instrument. The companion paper^{15} presents the application of this model to extract intrinsic fluorescence from human breast tissue fluorescence spectra.

## 2.

## Methods

## 2.1.

### Theory

## 2.1.1.

#### Foundation of our fluorescence model

The approach for extracting intrinsic fluorescence from a homogeneous, semiinfinite, turbid medium is achieved in two steps. First, a scalable Monte-Carlo-based model of diffuse reflectance previously developed by our group^{16, 17} is used to extract the absorption and scattering coefficients (optical properties) of the medium from a measured diffuse reflectance spectrum. Second, the optical properties are incorporated into a scalable Monte Carlo model of fluorescence to extract the intrinsic fluorescence from the raw turbid medium fluorescence. Scaling is critical to make the Monte Carlo model an efficient inversion tool for modeling tissue fluorescence.

The foundation for our fluorescence model is the work carried out previously by Swartling
^{18} Swartling outlined a technique wherein Monte Carlo simulations of fluorescence can be performed efficiently.^{18} Their approach is to break the Monte Carlo simulation into two separate simulations, one dealing with the excitation light traveling from the light source to the tissue fluorophore and one dealing with the emitted fluorescence traveling from the fluorophore to the detector. Reciprocity can also be taken advantage of to require only a single simulation for both excitation and emission. They demonstrate the use of a “reverse emission” simulation where the excitation simulation is run, and reciprocity is used to determine the emission collection probability. Each of the two baseline simulations (absorption and emission) can then be scaled to any arbitrary set of optical properties using a variety of approaches. One scaling approach that they presented was to run the baseline simulation with zero absorption, and then use Beer’s law to scale to any absorption coefficient by simply recording the path length of each photon. They also mention the possibility of a scaling approach for different scattering properties, but did not demonstrate this in their manuscript.

We have developed a fluorescence model for a semiinfinite medium, with illumination and collection at the surface of the medium, which builds on the work by Swartling
^{18} and Graaff
^{19} Figure 1
shows a flowchart of the model. Briefly, the inputs into the model are the absorption and scattering coefficient pairs of the medium at the excitation wavelengths [
${\mu}_{a}\left({\lambda}^{x}\right)$
,
${\mu}_{s}\left({\lambda}^{x}\right)$
] and emission wavelengths [
${\mu}_{a}\left({\lambda}^{m}\right)$
,
${\mu}_{s}\left({\lambda}^{m}\right)$
], the experimentally measured fluorescence spectrum
$\left({F}_{\mathrm{meas}}\right)$
; and the fiber geometry collection efficiency
$P\left(r\right)$
, where
$r$
is the radial distance from source to exit traveled by each simulated photon. The absorption and scattering coefficients at the excitation wavelength are used to determine the absorbed energy density probability grid
$\mathit{A}(r,z)$
, where
$(r,z)$
is the point of absorption in a cylindrical coordinate system, and
$\mathit{A}(r,z)$
gives the probability per unit volume that a photon will be absorbed at the specified coordinate. The absorption and scattering coefficients at the emission wavelength are used to determine the fluorescence escape density probability grid,
$\mathit{E}(r,z)$
, defined as the probability per unit area that a photon originating at the coordinate
$(0,z)$
,will exit the surface of the medium at coordinate
$(r,z)$
. The measured fluorescence, the absorbed energy grid, the fluorescence escape density grid, and the fiber probe geometry are incorporated into in the model to determine the intrinsic fluorescence properties of a given fluorophore, which are described by the quantum yield
$\left(\varphi \right)$
, concentration
$\left(C\right)$
, extinction coefficient
$\epsilon \left(\lambda \right)$
, and emission line shape
$\eta \left(\lambda \right)$
. The intrinsic fluorescence properties are independent of the effects of absorption and scattering and the instrument configuration.

While our model builds on the work by Swartling, ^{18} it has a few notable differences. The fluorescence escape density function grid is simulated, and reciprocity is used to determine the absorption energy density grid which is the opposite of that carried out by Swartling
^{18} This approach was found to yield better agreement with the original Monte Carlo simulations. Next, scaling relations described by Graaff were incorporated into this model.^{19} These scaling relations are used to compute the absorption and fluorescence escape density grids separately for any set of optical properties. This is performed by storing sufficient information about the path of each photon to enable its exit weight and location to be scaled to any set of optical properties. Specifically, the exit weight of each photon is given by
${W}_{\mathrm{new}}={W}_{\mathrm{sim}}({a}_{\mathrm{new}}\u2215{a}_{\mathrm{sim}}{)}^{N}$
, and the exit location is given by
${r}_{\mathrm{new}}={r}_{\mathrm{sim}}({\mu}_{t,\mathrm{sim}}\u2215{\mu}_{t,\mathrm{new}})$
, where
$a$
is the albedo,
$W$
is the photon weight,
$N$
is the number of interactions the photon underwent within the medium,
${\mu}_{t}$
is the transport coefficient
$({\mu}_{s}+{\mu}_{a})$
, and
$r$
is the radial distance from source to exit traveled by each photon, with the parameters specified for both the original simulation (sim) and the desired optical properties (new). These scaling relationships describe the travel of a photon between two locations in tissue. To handle differences in the anisotropy factor
$g$
the scattering coefficient
${\mu}_{s}$
was scaled such that the reduced scattering coefficient
${\mu}_{s}(1-g)$
was conserved, as follows,
${\mu}_{s,\mathrm{new},\mathrm{scaled}}={\mu}_{s,\mathrm{new}}(1-{g}_{\mathrm{new}})\u2215(1-{g}_{\mathrm{sim}})$
. These scaling relationships have the advantage of being able to more efficiently describe the effects of absorption outside of the diffusion regime, by taking into account the change in mean step size, which is proportional to
$1\u2215{\mu}_{t}$
of surviving photons associated with higher absorption coefficients, as compared to a simple Beer’s law scaling, which implicitly assumes the step sizes are fixed and proportional to
$1\u2215{\mu}_{s}$
.

The second notable difference is that the quasidiscrete Hankel transform^{20} was used to convolve the absorption energy density and fluorescence escape function grids, enabling greater speed, which is important for clinical applications. The Hankel transform is equivalent to the 2-D Fourier transform of a circularly symmetric function. This enables convolution of the absorption energy density and fluorescence escape function simulations to be performed efficiently. We followed the approach of Li, ^{20} who outline a computationally efficient means by which this can be performed. A third new feature of this work is that convolution was also used to model the optical probe geometry. This enables this model to be applied to any arbitrary probe configuration.

Finally, a framework was developed by which these forward simulation techniques can be applied to tissue spectra, by incorporating a scalable Monte Carlo model to calculate optical properties at both the excitation and emission wavelengths, which can then be used to solve for the intrinsic fluorescence properties of the tissue or other media independent of absorption and scattering. The fast, scalable Monte Carlo model extracts optical properties from the measured reflectance spectrum of the medium.^{16} The model makes it a must to know *a priori* what absorbers are present in the tissue and assumes that the scattering can be approximated by Mie theory scattering properties monodisperse, spherical particles.^{21} The diffuse reflectance spectrum measured from the medium is calibrated to the diffuse reflectance spectrum from a single phantom reference measurement of “known” optical properties to calibrate the system throughput and wavelength dependence. This ratio enables a direct comparison of the measured spectra to that simulated by the scalable Monte Carlo model. A nonlinear optimization algorithm is used to minimize the sum of squares error between the measured and modeled diffuse reflectance spectra by varying the absorber concentration and scatterer size/density to extract the optical properties of the medium.

## 2.1.2.

#### Implementation of the fluorescence model

First, Monte Carlo simulations were run to generate the emitted fluorescence traveling from the emitter to the detector, i.e., an escaping fluorescence energy density map $\mathit{E}(r,z)$ , as already defined. Photons were launched over a series of depth increments, ranging from 0 to $1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ , in 0.005-cm increments, with 100,000 photons being launched at each increment. This provides an increment approximately on the order of the mean free path in tissue. A depth of $1\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ was chosen, since depths greater than this were found to have negligible contribution to the measured fluorescence for the probe geometry used in this study. The following parameters were used: absorption coefficient ${\mu}_{a}=0.01\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ; scattering coefficient ${\mu}_{s}=150\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ ; anisotropy factor $g=0.8$ ; refractive index of medium; 1.335 (corresponds to that of the phantoms); and refractive index of fused silica fiber; 1.47. Photons that escape to the surface of the medium and which are within the numerical aperture (NA) of the detector fiber (in our case, the NA is 0.12) were stored in an array, containing the exit weight, launch depth $\left(z\right)$ , net radial distance traveled $\left(r\right)$ , and number of interactions $\left(N\right)$ within the medium prior to escape.

Next, using the principle of reciprocity, it is possible to determine the course of the excitation light traveling from the light source to the tissue fluorophore, i.e., how much light energy would be absorbed at a given location within the medium were the photon launched at the surface of the tissue using a fiber with the same NA. However, note that one must account for the difference in solid angle between the acceptance/emission of light by the absorber/fluorophore within the tissue
$\left(4\pi \right)$
and the limited acceptance angle of the source/detector fiber (defined by the NA). The acceptance/emission solid angle for a single molecule would not be
$4\pi $
, but for a suspension of randomly oriented molecules, the collective acceptance/emission is uniform. Swartling derived the effect of a difference in solid angle between the source and detector and the acceptance/emission of light by the absorbers/fluorophores,^{18} and using this factor
$(4\pi \u2215\mathrm{\Delta}\Omega )$
, the relationship between
$\mathit{E}(r,z)$
and
$\mathit{A}(r,z)$
is given by

Not all of the absorbed energy is converted to fluorescence. The probability of an absorbed photon generating a fluorescent photon at a given emission wavelength is given by Swartling
^{18} as

## Eq. 2

$${\varphi}_{\mathrm{eff}}({\lambda}^{x},{\lambda}^{m})=\varphi \frac{{\mu}_{a}^{f}\left({\lambda}^{x}\right)}{{\mu}_{a}\left({\lambda}^{x}\right)}\frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}{\lambda}^{m}}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\mathrm{d}{\lambda}^{m}},$$Given $\mathit{A}(r,z)$ and Eq. 2, it is then possible to determine the generated fluorescence energy density at all points within the medium for any given set of optical properties at the excitation wavelength and for a given set of fluorescence properties [ $\varphi $ , ${\mu}_{a}^{f}\left({\lambda}^{x}\right)$ , $\eta \left({\lambda}^{m}\right)$ ] as follows:

where ${\mathit{F}}_{\mathrm{gen}}(r,z)$ is the fluorescence energy density generated within the medium as a function of depth $z$ and radial distance $r$ . The fluorescence generated at a given depth is then convolved with the fluorescence escape energy density to obtain the exit probability of fluorescent photons escaping the surface of the medium to be collected. Summing over the entire range of depths within the medium gives the total emitted fluorescence as a function of radial distance. Thus, the emitted fluorescence is given by## Eq. 4

$${F}_{\mathrm{exit}}\left(r\right)=\sum _{j}\mathrm{\Delta}{z}_{j}{F}_{\mathrm{gen}}(r,{z}_{j})*E(r,{z}_{j}),$$The inverse problem is solved by noting that since ${\varphi}_{\mathrm{eff}}$ is a scalar, owing to the associative property of convolution, substituting Eq. 3 into Eq. 4, it can be taken out of the summation in Eq. 4, as

## Eq. 5

$${F}_{\mathrm{exit}}\left(r\right)={\varphi}_{\mathrm{eff}}({\lambda}^{x},{\lambda}^{m})\sum _{j}\mathrm{\Delta}{z}_{j}A(r,{z}_{j})*E(r,{z}_{j}).$$The terms within the summation depend on only the optical properties of the medium at the excitation and emission wavelengths, and not on the fluorescence properties of the medium. In particular, $\mathit{A}(r,z)$ is a function of the optical properties at the excitation wavelength, and $\mathit{E}(r,z)$ is a function of the optical properties at the emission wavelength. Furthermore, note that in Eq. 2, for a known fluorophore, the wavelength-dependent extinction coefficients and emission spectra are known. Thus, ${\varphi}_{\mathrm{eff}}({\lambda}^{x},{\lambda}^{m})$ can be expressed as the product of the fluorophore concentration, the quantum yield, and a set of wavelength dependent constants, as

## Eq. 6

$${\varphi}_{\mathrm{eff}}({\lambda}^{x},{\lambda}^{m})=\varphi \frac{{\mu}_{a}^{f}\left({\lambda}^{x}\right)}{{\mu}_{a}\left({\lambda}^{x}\right)}\frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}\lambda}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}\lambda}=\varphi \frac{2.303\epsilon \left({\lambda}^{x}\right)C}{{\mu}_{a}\left({\lambda}^{x}\right)}\frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}{\lambda}^{m}}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\lambda}^{m}},$$Because
$E(r,z)$
, and hence
$A(r,z)$
, were simulated using a point source, the convolution of the two terms results in the radial distribution of fluorescence occurring using a point source illumination. To account for the probe, one must then multiply this with the probe specific collection probability as a function of radial distance and sum over all distances. The collection probability for a probe geometry consisting of separate illumination and collection fibers is given in Refs. 16, 22 (note errata^{22} to original manuscript^{16}) as

## Eq. 7

$$P\left(r\right)=\frac{2}{{\pi}^{2}{R}_{i}^{2}}{\int}_{\mathrm{max}(-{R}_{i},s-r-{R}_{c})}^{\mathrm{min}({R}_{i},s-r+{R}_{c})}(s-x)\phantom{\rule{0.2em}{0ex}}{\mathrm{cos}}^{-1}\phantom{\rule{0.2em}{0ex}}\left[\frac{{s}^{2}+{(s-x)}^{2}-{R}_{i}^{2}}{2(s-x)s}\right]{\mathrm{cos}}^{-1}\phantom{\rule{0.2em}{0ex}}\left[\frac{{r}^{2}+{(s-x)}^{2}-{R}_{c}^{2}}{2(s-x)r}\right]\phantom{\rule{0.2em}{0ex}}\mathrm{d}x,$$## Eq. 8

$${F}_{\mathrm{meas}}({\lambda}^{x},{\lambda}^{m})=S\varphi \frac{2.303C\epsilon \left({\lambda}^{x}\right)}{{\mu}_{a}\left({\lambda}^{x}\right)}\frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}{\lambda}^{m}}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\lambda}^{m}}\underset{i}{\sum}\left\{{a}_{i}P\left({r}_{i}\right)\sum _{j}[\mathrm{\Delta}{z}_{j}A({r}_{i},{z}_{j})*E({r}_{i},{z}_{j})]\right\},$$## Eq. 9

$$2.303C\epsilon \left({\lambda}^{x}\right)\varphi \frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}{\lambda}^{m}}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\lambda}^{m}}=\frac{{\mu}_{a}\left({\lambda}^{x}\right){F}_{\mathrm{meas}}({\lambda}^{x},{\lambda}^{m})}{S\sum _{i}\left\{{a}_{i}P\left({r}_{i}\right)\sum _{j}[\mathrm{\Delta}{z}_{j}A({r}_{i},{z}_{j})*E({r}_{i},{z}_{j})]\right\}}.$$
${a}_{i}$
is the surface area of each radial grid element. Just to clarify, all of the parentheses indicate functional relationships or summations, and are not products. Also
${F}_{\mathrm{meas}}({\lambda}^{x},{\lambda}^{m})$
is the fluorescence measured using a particular probe geometry, given the optical
$\left[{\mu}_{a}\right({\lambda}^{x},{\lambda}^{m}),{\mu}_{s}^{\prime}({\lambda}^{x},{\lambda}^{m}\left)\right]$
and fluorescence properties
$[\varphi ,\epsilon ({\lambda}^{x}),\eta ({\lambda}^{m}),C]$
of the medium, and
$S$
is introduced here as a scaling factor necessary to account for the difference in magnitude between the Monte Carlo simulations, which are on an absolute scale, and the measurement, which is typically relative to a fluorescence standard. This factor must be determined using a phantom measurement to pool data collected using different instruments or probe geometries. For a single instrument setup, it can be set to unity to leave the intrinsic fluorescence spectra on a relative scale. The quantum yields, extinction coefficients, and emission efficiencies could be independently determined, allowing an absolute determination of fluorophore concentration in tissue. The result of this model thus allows for the extraction of intrinsic fluorescence line shape, and intensity, independent of the effects of absorption and scattering. This quantity is proportional to the product of the concentration and fluorescence quantum yield. This can in turn be used to determine the product of the fluorophore concentration provided that some *a priori* information is known about the fluorophore.

## 2.2.

### Phantom Validation

The model just described was tested using synthetic liquid tissue phantoms. Liquid phantoms were made using hemoglobin (H0267, Sigma-Aldrich Corp.) as the absorber, furan-2 as the fluorophore, and polystyrene spheres (07310, Polysciences, Inc.) as the scatterer. Three sets of phantoms were prepared, having low, medium, and high scattering properties. Into each phantom,
$12\phantom{\rule{0.3em}{0ex}}\mathrm{\mu}\mathrm{g}\u2215\mathrm{ml}$
of furan-2 was added (which contributes negligible absorption), and varying concentrations of hemoglobin were added in to change the absorption properties of the medium, yielding a total of 11 phantoms. One phantom was an outlier, having significantly lower (25 to 50%) fluorescence and reflectance intensity relative to the other phantoms (likely due to incomplete probe contact) and was excluded. The absorption coefficient of hemoglobin was determined using a spectrophotometer (Cary 300, Varian, Inc.), and the scattering coefficient of the polystyrene spheres was calculated using Mie theory,^{21, 23} taking into account the wavelength dependent real refractive index of polystyrene spheres^{24} and water.^{25} Both the spheres and water were assumed to be nonabsorbing in the analysis, although the effects of absorption by spheres were evaluated separately. Table 1
summarizes the optical properties of the phantoms over the wavelength range of 330 to
$600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. Phantoms 1 to 4 correspond to low scattering, 5 to 8 correspond to medium scattering, and 9 to 11 correspond to high scattering. Within each scattering level, as the phantom number increases, more hemoglobin has been added, so the absorption coefficients increase, while the scatterer is diluted slightly, so the scattering coefficients decrease. These phantoms are designed to mimic tissue optical properties in the UV to visible spectral range, and specifically correspond to the optical properties we have observed in breast tissues.^{17} Phantom 6 was used as the reference phantom for the purposes of extracting the optical properties of all phantoms because it is in the middle of the range for both absorption and scattering.

## Table 1

Summary of the optical properties of the phantoms used to validate the fluorescence model over the wavelength range of 330 to 600nm .

Absorption Coefficient μa | Reduced Scattering Coefficient μs′ | |||||||
---|---|---|---|---|---|---|---|---|

Phantom Number | Min (cm−1) | Max (cm1) | Mean (cm−1) | Mean Relative Error of Extraction (%) | Min (cm1) | Max (cm1) | Mean (cm1) | Mean Relative Error of Extraction (%) |

1 | 0.0 | 0.9 | 0.2 | 26.1 | 3.9 | 5.9 | 4.6 | 4.9 |

2 | 0.1 | 6.8 | 1.4 | 24.9 | 3.8 | 5.8 | 4.5 | 8.0 |

3 | 0.1 | 12.4 | 2.6 | 27.1 | 3.7 | 5.7 | 4.4 | 12.3 |

4 | 0.2 | 22.9 | 4.8 | 26.5 | 3.6 | 5.4 | 4.2 | 16.1 |

5 | 0.0 | 0.9 | 0.2 | 11.9 | 11.7 | 17.7 | 13.9 | 0.8 |

6 | 0.0 | 6.8 | 1.4 | 0.0 | 11.4 | 17.3 | 13.5 | 0.0 |

7 | 0.0 | 12.4 | 2.6 | 4.4 | 11.2 | 16.9 | 13.2 | 1.9 |

8 | 0.0 | 22.9 | 4.8 | 5.8 | 10.7 | 16.1 | 12.6 | 3.2 |

9 | 0.0 | 6.7 | 1.4 | 7.0 | 19.0 | 28.7 | 22.4 | 18.8 |

10 | 0.0 | 12.3 | 2.6 | 1.2 | 18.5 | 28.1 | 21.9 | 18.8 |

11 | 0.0 | 22.7 | 4.8 | 2.4 | 17.7 | 26.8 | 21.0 | 16.7 |

The SkinSkan fluorometer (J.Y. Horiba) and fiber optic probe were used for making all experimental diffuse reflectance and fluorescence measurements from the synthetic liquid tissue phantoms. Briefly, the instrument consists of a 150-W xenon lamp, dual excitation and emission grating monochromators, and a photomultiplier tube (PMT) detector. Diffuse reflectance measurements were collected using a “synchronous scan,” wherein the excitation and emission monochromators are scanned in tandem across the wavelength range of interest (330 to $600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ). Fluorescence emission spectra were acquired by fixing the excitation monochromator to $330\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ and scanning the emission monochromator from 350 to $500\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . All scans were run with a 5-nm increment. The fiber optic probe consisted of a central collection core with a diameter of $1.52\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , surrounded by an illumination ring with an outer diameter of $2.18\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Both the illumination core and surrounding ring consisted of 31 fibers, each of which had a core/cladding diameter of $200\u2215245\phantom{\rule{0.3em}{0ex}}\mathrm{\mu}\mathrm{m}$ .

The diffuse reflectance (330 to $600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ ) and fluorescence (330-nm excitation) were measured from each sample by placing the probe submerged just below the surface of the phantom. The tip of the probe was sealed with epoxy between the fibers, preventing any capillary action. Each phantom was placed in a container of approximately 1.5-cm diameter, with a total volume of $8\phantom{\rule{0.3em}{0ex}}\mathrm{ml}$ , which satisfied the criterion of a semiinfinite medium for the model. This was tested by measuring with the probe at the surface of the medium, and submerged halfway down, with no observable difference in the signal intensities. Then the scalable Monte Carlo model was used to extract the optical properties of each sample from the measured diffuse reflectance using procedures described in Ref. 16. The extracted optical properties were input to the fluorescence model to retrieve the intrinsic fluorescence from the measured phantom fluorescence. In addition, the “known” optical properties, taken from Mie theory and the spectrophotometer measurement, were also used as inputs into the fluorescence model to retrieve the intrinsic fluorescence. This allows for an assessment of the propagation of errors from the reflectance model to the fluorescence model. To illustrate the effectiveness of the model in determining the intrinsic fluorescence properties within the medium, Eq. 9 was solved for the product:

## Eq. 10

$$\varphi \frac{\eta \left({\lambda}^{m}\right)\mathrm{\Delta}{\lambda}^{m}}{{\int}_{0}^{\infty}\eta \left({\lambda}^{m}\right)\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\lambda}^{m}}.$$## 3.

## Results

Table 1 shows the extraction errors for each of the phantom’s optical properties, as compared to the expected values. Phantom 6, a middle-level scatterer, was used as the reference phantom for the extraction, hence it has low errors, as do the phantoms with similar optical properties. The low- and high-scattering phantoms each show larger errors, relative to the middle-scattering level. There are relatively smaller increases in the extraction errors as the absorption coefficient is varied for the same scatterer level. This indicates that varying the density of polystyrene spheres causes larger deviations in the extracted optical properties than does varying the absorber concentration. This may be attributed to uncertainty in the optical properties of polystyrene spheres, as discussed further in the following.

Figure 2 shows the product in Eq. 10 for the nonnormalized [Figs. 2a, 2b, 2c] and normalized [Figs. 2d, 2e, 2f] fluorescence spectra, acquired at 330-nm excitation. The nonnormalized spectra plots, shows the raw fluorescence spectra [Fig. 2a], intrinsic fluorescence using the extracted optical properties [Fig. 2b], and intrinsic fluorescence using the “known” optical properties [Fig. 2c]. The normalized plots show the same data sets in Figs. 2d, 2e, 2f], respectively, normalized to their peak intensity. The fluorescence spectral line shape acquired with no absorber or scatterer present is also shown in the normalized plots (black line) and can be seen to be very similar to the extracted line shape for the cases where the known, and extracted optical properties were used. We can see that the model corrects for the original large differences in magnitude and line shape present due to the differences in absorption and scattering in each of the phantoms.

Table 2 quantifies the errors in line shape and intensity extracted for the 11 phantoms. Line shape was evaluated by calculating the fraction of variance accounted for by the line shape of the fluorophore that is measured in a dilute solution (this provides the known intrinsic line shape), using a linear fit. This is given by the relationship 1-SSE/SST, where SSE is the sum of squares of the residuals in this fit, and SST is the sum of squares of the original spectra. This represents the ${R}^{2}$ values, all of which are above 0.98. The percent deviation from the mean peak intrinsic fluorescence intensity of furan for all phantoms [determined from Eq. 10] was $9.3\pm 9.8\%$ for the case where the optical properties were extracted from the diffuse reflectance and $9.1\pm 4.7\%$ where the known optical properties were used. The phantoms with the lowest scattering properties (1 to 4) appear to have generally higher errors than the more scattering phantoms. This is significant for the use of known optical properties ( $p<0.001$ using a $t$ test), but not significant for the use of extracted optical properties.

## Table 2

The variance of the intrinsic line shape is well described by the fluorescence spectrum of furan-2 in a dilute solution, as indicated by high R2 values; and the percent deviations from the mean intensity for each of the phantoms.

Phantom | Line Shape—Extracted Optical Properties (R2) | Line Shape—Known Optical Properties (R2) | Deviation in Intensity—Extracted Optical Properties (%) | Deviation in Intensity—Known Optical Properties (%) |
---|---|---|---|---|

1 | 0.99 | 0.99 | $-23.8$ | $-13.5$ |

2 | 0.98 | 0.99 | $-12.8$ | 12.0 |

3 | 0.98 | 0.99 | $-18.4$ | 16.0 |

4 | 0.98 | 0.98 | $-27.0$ | 12.3 |

5 | 0.99 | 0.99 | $-8.1$ | $-12.3$ |

6 | 0.99 | 0.99 | 0.2 | 0.2 |

7 | 0.99 | 0.99 | $-1.2$ | 3.4 |

8 | 0.99 | 0.99 | $-2.3$ | 6.1 |

9 | 1.00 | 1.00 | 5.8 | $-8.7$ |

10 | 1.00 | 1.00 | 0.3 | $-7.5$ |

11 | 0.99 | 1.00 | $-2.7$ | $-8.0$ |

Following up on this finding, we evaluated the effects of the assumption that the polystyrene spheres contribute negligible absorption to the phantoms. There are no data to our knowledge characterizing the complex refractive index of polystyrene spheres at UV wavelengths down to
$330\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. However, using the findings of Ma
^{24} as a guide, we experimented with two different imaginary components of the refractive index,
$-0.0004$
, which is approximately the value Ma reported throughout the visible,^{24} and
$-0.0001$
, which corresponds to lower absorption, which was arbitrarily chosen. Table 3
shows the results of using these imaginary components in calculation of the absorption and scattering properties of these phantoms. We can see that the incorporation of absorption by the spheres substantially improves the accuracy with which the intrinsic fluorescence intensity can be extracted for the “known” optical properties, but not in the case of the extracted optical properties. This may be due to the inverse problem becoming more poorly conditioned when the concentration of polystyrene spheres affects both the absorption and scattering properties, resulting in greater errors in the extracted optical properties. Interestingly, using the imaginary component at only the excitation wavelength results in similar error reduction as the case where it is used at both the excitation and emission wavelengths.

## Table 3

Mean relative error of the extracted fluorescence intensity using the “known” and extracted optical properties.

Imaginary Refractive Index of Polystyrene Spheres | Mean Relative Error Using “Known” Optical Properties | Mean Relative Error Using Extracted Optical Properties |
---|---|---|

0 (nonabsorbing) | $9.1\pm 4.7\%$ | $9.3\pm 9.8\%$ |

$-0.0001$ at all $\lambda $ | $7.0\pm 4.0\%$ | $10.5\pm 10.8\%$ |

$-0.0004$ at all $\lambda $ | $4.5\pm 4.7\%$ | $15.7\pm 12.6\%$ |

$-0.0001$ at excitation only | $7.7\pm 4.2\%$ | $10.0\pm 10.2\%$ |

$-0.0004$ at excitation only | $4.4\pm 4.1\%$ | $12.8\pm 11.4\%$ |

$-0.0001$ at emission only | $8.5\pm 4.5\%$ | $9.7\pm 10.3\%$ |

$-0.0004$ at emission only | $6.5\pm 4.1\%$ | $11.4\pm 11.3\%$ |

## 4.

## Discussion

This paper outlined a methodology by which intrinsic fluorescence spectra can be extracted from combined fluorescence and diffuse reflectance spectra, which are commonly measured from tissues in a variety of clinical and preclinical studies. This enables more biologically relevant interpretation of these data, since they can be directly related to underlying tissue properties (fluorophore concentrations or microenvironment). The significant advantages of this model are (1) it does not require any assumptions regarding the underlying fluorescence properties in order to calculate the intrinsic fluorescence spectra; (2) it is able to account for any arbitrary probe geometry; (3) it does not require extensive empirical calibration; (4) it is appropriate to use for cases of high absorption and small source-detector separation; and (5) since the intrinsic fluorescence properties are instrument independent, it allows for data acquired from multiple instruments and/or probe geometries to be pooled. This pooling would require the use of a phantom measurement to determine the instrument dependent factor $S$ , as already described, or in lieu of that, the data could be normalized and pooled on a relative scale to eliminate this requirement. This model can be used to extract endogenous sources of fluorescence as well as characterize the prevalence of exogenous sources of fluorescence in turbid media.

The application of this model to phantoms was relatively straightforward, since all of the constituents were known and well characterized, and reasonably good accuracy was obtained in extracting the intrinsic fluorescence line shape and intensity. This model performs well for the entire range of optical properties tested. A significant source of error appears to relate to the accuracy with which the optical properties of polystyrene spheres are quantified. This is reflected by the fact the intrinsic fluorescence intensity extraction accuracy was highly sensitive to the imaginary refractive index of polystyrene used in the model (see Table 3). Unfortunately, there is limited information on the optical properties of polystyrene spheres, and these optical properties may vary depending on the source, since absorption is likely highly influenced by contaminants (see Ref. 26 for a detailed discussion). There is no information to our knowledge reporting the optical properties of polystyrene spheres at the excitation wavelength of
$330\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
, so it was necessary to extrapolate the findings of Ma
^{24} It was found that using an imaginary refractive index of
$1e\text{-}4$
, substantially reduced the errors in the extracted intrinsic fluorescence intensity from
$9.1\pm 4.7$
down to
$4.5\pm 4.7$
(see Table 3). It appears that the errors found in the model are most highly influenced by the optical properties at the excitation wavelength. Using the complex refractive index at the excitation wavelength alone, produced the same reduction in error as using it over the entire wavelength range. It is also possible that the real component of the refractive index is inaccurate at the excitation wavelength, as we had to employ extrapolation here as well. Any such analysis is rather speculative, but serves to point out potential sources of error.

Another potential source of error is the fact that the Henyey-Greenstein phase function was used in the Monte Carlo simulations, rather than the Mie phase function. The robustness of this model is well demonstrated by the fact that, despite these uncertainties, in the simplest case of nonabsorbing spheres, the model was found to be able to extract the intrinsic fluorescence intensity with mean percent errors of $9.1\pm 4.7\%$ using the “known” optical properties, and $9.3\pm 9.8\%$ using the extracted optical properties, thus demonstrating the usefulness of this model for a wide range of optical properties. Errors are reported in terms of the variability of the intensity of the intrinsic fluorescence. However, note that if the quantum yield, extinction coefficient, and emission spectrum could be accurately determined, or could be assumed to be constant, equivalent relative errors would be obtained in extracting the concentration of the fluorophore.

Limitations of this model include the fact that it assumes a homogeneous distribution of fluorophores with homogeneous optical properties. However, it is possible to adapt this approach to accommodate a nonhomogeneous distribution, particularly for layered media. This model as implemented also requires knowledge of the absorbers and scatterers present in the medium to retrieve the absorption and scattering properties of the medium, although other techniques for optical property retrieval could be used for this purpose that do not impose such a limitation.^{27} This model could also be made computationally more efficient by running only a single simulation similar to the “reverse emission” approach outlined by Swartling
^{18} (described in the methods section of this manuscript) that could also record information required for the diffuse reflectance baseline simulation, simply by also recording photons exiting the medium.

In the companion manuscript we will explore the application of this model to fluorescence spectra measured from human breast tissues.^{15} Application to tissue is more challenging due to the presence of multiple absorbers, scatterers, and fluorophores of uncertain properties. However, we have previously applied the scalable Monte Carlo reflectance model to extract optical properties from breast tissue diffuse reflectance spectra,^{17} and the addition of the fluorescence model is relatively straightforward since it does not require any assumptions regarding the properties of the fluorophores in order to retrieve the intrinsic fluorescence. The intrinsic fluorescence defined by Eq. 10 was calculated, taking
$S$
to be unity. This enables comparison of relative fluorescence intensities in tissue, where absolute calibration would be extremely difficult due to differences in microenvironment etc. that may affect fluorescence properties. In addition, the presence of multiple fluorophores requires that one separate the influence of each fluorophore spectrally. In the companion paper, multivariate curve resolution^{28} was utilized to achieve this. It may be desirable to obtain the actual concentration of the fluorophore in future applications. However, this does require some additional knowledge of the fluorophore properties that steady state fluorescence measurements alone cannot provide, particularly the quantum yield.

## Acknowledgments

This work was funded by University of Wisconsin Radiological Sciences Training Grant No. 5T32CA009206-27, sponsored by the National Institutes of Health, the Department of Health and Human Services, and the Public Health Service. Additional funding was provided by the National Institutes of Health through Grant No. 1R01CA100559-01A1 and the Department of Defense Breast Cancer Research Program Grant No. DAMD17-02-1-0628.