Recovering intrinsic fluorescence by Monte Carlo modeling

Fluorescence, the re-emission of light by some molecules (called fluorophores) that have absorbed light of a shorter wavelength, is widely used to investigate biological tissues.^1,2 It can be induced by native fluorophores (autofluorescence) or by administered fluorophores. The presence of fluorophores, their concentration, and their changes over time are indicative of the tissue state. Fluorescence spectroscopy can therefore be used to identify different tissue types and pathological states. It complements other optical spectroscopic techniques, such as diffuse reflectance spectroscopy (DRS) or Raman spectroscopy.

Each fluorophore has a characteristic emission spectrum that depends on the excitation wavelength. By comparing or fitting the emitted spectra of individual fluorophores to a measured fluorescence spectrum of tissue the concentration of the fluorophores present in the tissue can, in principle, be deduced. This is not straightforward, however, since the measured fluorescence spectrum will be strongly distorted by scattering and absorption, both at the excitation and at the emission wavelengths. As a consequence, quantitative information about fluorophore concentrations can usually only be gained if correction techniques are used to compensate for the effects of scattering and absorption. This is called recovering the intrinsic fluorescent spectra (IFS) of the tissue. The intrinsic fluorescence is the fluorescence that is only due to fluorophores independent of the interference of absorption and scattering.

Several techniques to recover the intrinsic fluorescence have been developed. Bradley et al. reviews over 50 different publications that addressed the recovery of intrinsic fluorescence properties.³ Although promising results have been obtained with various methods, thus far all available methods have their limitation. For example, probably one of the best available methods for recovering the intrinsic fluorescence spectrum was developed by Zhang, Müller et al.^4,5 based on photon-migration theory. The method requires a DRS spectrum measurement taken at the same location with the same measurement geometry as used for the fluorescence measurement. The requirement that the DRS spectrum has to be taken with the same measurement geometry is, however, limiting its applicability in practice. Autofluorescence signals are normally faint compared to DRS signals and are weakened further by absorption in the tissue. It is therefore important for the source-detector fiber distance to be as small as possible to achieve a high signal-to-noise ratio. Conversely, DRS measurements typically achieve better contrast with longer optical path lengths. Also, the widely used diffusion approximation is not applicable for short distances.⁶ Therefore, in practice the requirement to measure fluorescence and DRS spectra under the same geometry requires either a compromise in the spectra quality or an additional DRS measurement. Other very promising techniques were developed, e.g., by Sinaasappel and Sterenborg for multiple excitation wavelengths⁷ and by Palmer and Ramanujam based on fast Monte Carlo calculations,⁸ which is more flexible but too computationally expensive for real-time intrinsic fluorescence recovery.

In this paper, we propose an alternative way to recover the intrinsic fluorescence spectra, based on a look-up table derived from Monte Carlo calculations that allow real-time intrinsic fluorescence recovery. While this technique also needs a DRS measurement, it holds the promise to be applicable regardless of the measurement geometries for fluorescence and diffuse reflectance. This technique is validated on phantom measurements. For comparison, we also present a modification to the method developed by Zhang et al. and Müller et al. that allows intrinsic fluorescence reconstruction from fluorescence and DRS spectra measured with the same probe geometry but needs less prior information about probe-related parameters than the original technique. A comparison between the Monte Carlo look-up table and the photon-migration theory on various phantoms and tissues will be presented. Finally, we apply the modified photon-migration and the Monte Carlo look-up methods on biological tissue to evaluate their usefulness in practice.

When monochromatic light of intensity $Ix$ illuminates biological tissue via an optical fiber, fluorescence is emitted with the spectral distribution $fxm(λm)$. The indices $x$, $m$, and $xm$ are used for properties depending only on the excitation wavelength, only on the emission wavelength, and on both wavelengths, respectively. A fraction of the fluorescent light will be collected by a second fiber yielding the measured fluorescent light intensity $Fxm(λm)$. The spectral shape of $Fxm(λm)$ will be distorted relative to $fxm$ due to the influence of absorption and scattering. To determine the relative concentrations of the fluorophores present in the tissue, one first has to disentangle the effects of absorption and scattering from fluorescence by recovering the intrinsic fluorescence $fxm$ from $Fxm(λm)$. We also assume that a DRS spectrum is measured at the same location and that the spectral shapes $η(λm)$ of the various fluorophores are known.

A photon at the excitation wavelength $λx$ entering biological tissue typically will experience multiple scattering events until it is either absorbed or leaves the tissue. If it is absorbed by a fluorophore, then (with a probability of $ϕxm$) a new photon at the emission wavelength $λm$ is emitted. Again, this photon can be scattered multiple times before it is either absorbed or escapes the tissue so that it may be detected. A fluorescence photon may even be re-absorbed by a fluorophore resulting in secondary fluorescence, but this case is unlikely and shall be neglected here.

The scattering and absorption properties of a diffuse medium like biological tissue are characterized by the absorption coefficient $μa$, the scattering coefficient $μs$, and the anisotropy parameter $g$.⁹ The coefficients $μa$ and $μs$ are defined as the probability of photon absorption and scattering per unit path length, respectively, while the anisotropy parameter $g$ is defined as the average of the cosine of the scattering angle. In biological tissue, the anisotropy parameter is typically $g≈0.9$, which means predominantly forward scattering.¹⁰ Small changes of $g$ typically do not have a large influence on scattering as long as the reduced scattering coefficient $μs′=μs(1−g)$ (which gives the probability of equivalent isotropic photon scattering per unit length) remains constant. Therefore, we set $g=0.9$ for all calculations and simulations described in this paper. Scattering and absorption coefficients are strongly wavelength dependent, so the tissue properties at the excitation wavelength are described by $μax$, $μsx′$ and by $μam$, $μsm′$ at the emission wavelengths. The absorption coefficient from a given absorber is proportional to its concentration, $μa=∑ln(10)εici$, where $ci$ are the concentrations and $εi$ the extinction coefficients for the various absorbers present. In biological tissue, typically only a small fraction of the absorbed photons are actually absorbed by fluorophores.

Because fluorescence, absorption, and scattering are photon events, it is conceptually easiest to look at spectral probability distributions, i.e., the probability of an event for a photon with given wavelength. However, in practice one does not measure photon numbers or probabilities, but light intensities or energies. In this paper, spectral probability distributions are used when discussing principles and intensity spectra when showing actual measurement. Spectral probability distributions and intensity spectra are not proportional to each other due to the correction for the wavelength-dependent energy per photon $hc/λ$, where $h$ is the Planck constant and $c$ is the speed of light.

Following Zhang et al.,⁴ we define the intrinsic fluorescence $fxm$ as the total fluorescence intensity of a thin slab of thickness $l$ of nonabsorbing, nonscattering material containing the same concentration of fluorophores. The absorption coefficient due to the fluorophores only of this slab is denoted by $μafx$. If the slab is optically thin ($l≪1/μafx$) the intrinsic fluorescence intensity distribution can be written as Display Formula

To extract the intrinsic fluorescence $fxm$ from the extrinsic measured fluorescence $Fxm(λm)$, Zhang et al. and Müller et al. introduced a method based on photon-migration theory. Since, in this method the photon-migration due to scattering and absorption only will be used to disentangle these from the fluorescence measurement, it is essential that the diffuse reflectance and fluorescence are measured with the same geometry and at the same location. According to photon-migration theory, light traveling in the turbid medium is described in terms of photon propagation with discrete tissue interactions. These discrete events can either be scattering, absorption, or the emission of a fluorescence photon. A photon injected into the tissue from the delivery fiber can either be scattered a number of times before reaching the detection fiber, or be absorbed and be lost, or be absorbed and re-emitted as fluorescent photon. The probability that a photon enters the medium via the illumination fiber, is subject to $n$ scattering events, escapes the medium, and is collected at the detector fiber can be written as¹¹Display Formula

The measured diffuse reflectance of a homogenous medium is then given by Display Formula

The diffuse reflectance in absence of absorption ($a=1$), defined as $R0$, can be written as $R0=κ/ε$ with $ε=eβ−1$. Wu et al. showed that the parameter $β$ can be written as Display Formula

To describe the fluorescence, we first consider a dilute optical thin sample of thickness $l$, hence, $μfxl≪1$, having a homogeneous distribution of the fluorophore and no other absorbers and scatterers. The intrinsic fluorescence produced by this sample is given by Eq. (1). For fluorescence in turbid medium, the photon-migration can be divided into $i$ scattering events followed by a fluorescence event at $i+1$ and subsequent $n−i−1$ scattering events before reaching the detector fiber. The escape probability for these events can be approximated by⁴Display Formula

Similar as for the diffuse reflectance, the measured fluorescence intensity can then be written as Display Formula

The parameters $S=(1−g)/β$ and $l$ have to be determined for a given probe by validation on physical test phantoms as described by Müller et al. This method will be called the MüllerZhang method in this paper.

In practice we found that obtaining reliable values for $S$ and $l$ from phantom measurements is difficult. So to extract the intrinsic fluorescence in absence of prior knowledge of $l$ and $S$, we proceed as follows. For biological tissues, the anisotropy parameter $g$ is typically only slowly varying in the visible (VIS) and near infrared (NIR) part of the spectrum,¹⁰ so taking $β$ and $ε$ constant is a reasonable approximation. Furthermore, from Eq. (6) we deduce that the wavelength-dependent parameter $κ$ is determined by scattering only. Using the determined $R0$ and $am$ and substituting $κm=R0(eβm−1)$ in Eq. (6), we can deduce the parameter $βm$ and thus the parameter $εm$. From the measured $Rm$ and the determined diffuse reflectance in absence of scattering $R0m$ together with the parameter $εm$, we retrieve from Eq. (10) the intrinsic fluorescence from the measured fluorescence $Fxm$ apart from a scaling factor. Since for biological tissue we are interested primarily in the shape of the intrinsic fluorescence, we finally normalize the total retrieved intrinsic fluorescence to one. Throughout this paper this normalized retrieved intrinsic fluorescence shall be used for examples based on this photon-migration model. This variation of the MüllerZhang method will be referred to as the modified photon-migration (MPM) method.

The probability of an absorbed photon generating a fluorescent photon in the wavelength range between $λm$ and $λm+Δλm$ is given by the effective quantum yield $ϕxmeff$ and can be written as^8,12Display Formula

The intrinsic fluorescence $fxm$ is related to this effective quantum yield using Eq. (1) according to Display Formula

Furthermore, the probability $pxm$ that a photon from the excitation light is re-emitted at the wavelength $λm$ and then collected by the delivery fiber is given by Display Formula

If one ignores nonlinear effects (like photo-bleaching) $ϕxmeff$ and $pxm$ have to be proportional to each other: Display Formula

The factor $kxm$ represents the probability that a generated fluorescent photon is detected and depends strongly on the optical properties of the tissue at the emission wavelength ($μam$ and $μsm′$) and the probe properties (e.g., on the distance $ρ$ between the delivery and the collection fiber, the orientation, the diameters and the numerical apertures of the fibers, and the reflectivity of the probe material). It also depends on the tissue properties at the excitation wavelength, because $μax$ and $μsx′$ determine the location where the excitation photon is absorbed and therefore where the fluorescent photon is emitted.

For a fixed probe, $kxm$ depends only on the optical properties of the tissue at the excitation and the emission wavelengths: $kxm=kxm(μax,μsx′,μam,μsm′)$, hence we can write [using Eqs. (12) and (13)]: Display Formula

In principle, $kxm(μax,μsx′,μam,μsm′)$ can be determined experimentally, similar to an approach described by Rajaram et al. to determine the diffuse reflectance $R$ for a given probe at different $μa$ and $μs′$.¹³ In this approach, they made measurements on various phantoms with known optical properties to create a two-dimensional look-up table. In the case of fluorescence, $kxm$ depends on four variables, so the look-up table has to be four-dimensional. This makes determining a look-up table experimentally not only impractical but unfeasible. An analytical or numerical approach is required to obtain a reasonable sized look-up table for $kxm$. In this paper, a Monte Carlo approach will be presented.

An efficient Monte Carlo method to describe fluorescence has been developed in ⁸ based on earlier work presented in ¹². It considers fluorescence as a two-step process, with an excitation and an emission step. In the excitation step, the excitation light from the delivery fiber propagates through the tissue and part of it is absorbed by the fluorophore. In the emission step, the fluorophores emit light at the emission wavelength, which propagates through the tissue and is partly collected by the collection fiber. The results of both steps are then combined in a convolution procedure. Both steps are independent of each other and can be calculated separately with standard Monte Carlo routines as described by Wang et al.¹⁴ The excitation step depends only on the coefficient pair $μax$ and $μsx′$ while the emission step depends only on $μam$ and $μsm′$. The efficiency gain of the two-step model versus a single-step model (see Refs. ¹⁵–¹⁷) for creating a look-up table is enormous: it creates the four-dimensional look-up table out of two two-dimensional tables of Monte Carlo simulations. To create a 10 values each look-up table for $μax$, $μsx′$, $μam$, and $μsm′$ a single-step method needs to perform 10,000 Monte Carlo simulations, while the two-step method needs to do only 100 Monte Carlo simulations each for the emission and the absorption step.

We consider the case where an optical probe is in contact with semi-infinite tissue ($z>0$) with the circular delivery and collection fibers perpendicular to the tissue boundary at $z=0$. The origin of the coordinate system (0, 0) is set at the center of the exit facet of the delivery fiber. The light intensity then is rotational symmetric around the delivery fiber enabling the efficient use of cylindrical coordinates ($r,z$).

The excitation step is a standard Monte Carlo routine¹⁴ that calculates the absorbed energy probability density $A(r,z)$, which is the probability per unit volume that a photon from the source fiber will be absorbed at ($r,z$) for a given set of $μax$ and $μsx′$.

A photon emitted by a fluorophore has to escape the tissue before it can be detected. Due to the symmetry of the problem, the escape probability depends only on the depth $z$ at which the emitting fluorophore is located. The emission step calculates the fluorescence escape probability density $E(r,z)$, which is the probability per unit area that a photon emitted at ($0,z$) will exit the tissue at the coordinate ($r,0$) and is captured by a collection fiber taking the numerical aperture and diameter of the fiber into account. Fluorescence emission is isotropic, so each photon package starts in a random direction from ($0,z$). The photon package propagation is described by the absorption and scattering coefficients $μam$ and $μsm′$ of the tissue. Using an independent Monte Carlo step to calculate $E$ is a more general but also a more computationally expensive approach than using the principle of reciprocity to determine $E$ from $A$ as in Refs. ⁸ and ¹². And indeed in validation tests we obtained somewhat different values for $E$ from the two methods. We believe the reason for that is that the photon paths in the excitation and emission steps are not perfectly symmetric as assumed for reciprocity. For example, a fluorophore molecule close to the delivery fiber will have more photons coming from that direction of the fiber than from other directions, so the absorption can appear anisotropic while the emission of photons by the fluorophore is always perfectly isotropic. A similar argument can be made for the fibers: the light cone emitted from the delivery fiber will be symmetric but the collection fiber will accept any light within its numerical aperture (NA), whether symmetric or not.

$E$ and $A$ contain the spatial distribution of the fluorescence process. The probability that a photon of wavelength $λm$ is emitted at ($r,z$) is given by Display Formula

To obtain the probability of the generated fluorescent light being collected, this $pgen$ has to be convoluted with the fluorescence escape probability density $E$ which we shall do for each depth $z$ separately. The overall probability of fluorescence photons exiting at the surface at a given radius $r$ is then obtained by integrating over all depths:⁸Display Formula

This is Eq. (14) with an explicit expression for the correction factor $kxm$.

To calculate $kxm$, a convolution of $E$ and $A$ is required. For this paper, the convolution is done in Cartesian coordinates, so first $A(r,z)$ and $E(r,z)$ are coordinate transformed into $A(x,y,z)$ and $E(x,y,z)$. With the exit facet of the delivery fiber at (0, 0, 0) and the entry facet of the delivery fiber at ($ρ,0,0$) the equation for the correction factor becomes Display Formula

In a final step, $kxm$ is normalized by dividing it by the total absorption coefficient $μax$. This is necessary because the above derivation does not distinguish between absorption by fluorophores and by other chromophores.

The absorption, emission, and convolution steps described above were implemented in a Matlab (from MathWorks, Natick, Massachusetts) routine. All Monte Carlo calculations were done on 7 mm by 6 mm $r-z$ grids with a step-size of 0.02 mm. Test runs showed that photons that penetrate deeper than 6 mm or further than 7 mm from the source fiber have a negligible contribution to the measured fluorescence. The source and collection fiber had a numerical aperture of 0.22 and a core diameter of 200 μm. The look-up tables were calculated on a $28×20×28×20μsx′−μax−μsm′−μam$ grid. The grid values were 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1, 1.25, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 8, 10, and $12mm−1$ for $μsx′$ and $μsm′$ and 15, 8, 6, 5, 4, 3, 2.5, 2, 1.5, 1.25, 1, 0.75, 0.5, 0.25, 0.125, 0.1, 0.05, 0.025, 0.01, and $0mm−1$ for $μax$ and $μam$, respectively. This grid covers typical values encountered in tissue (a few $1mm−1$ for $μs′$ and of the order of few $0.1mm−1$ for $μa$). The Monte Carlo simulations were run for 250,000 photon packages both for the excitation and the absorption phase (except for $μam=0$, where the number was reduced to 100,000). Calculating the full look-up table took about three weeks on a standard office PC. Information about the validation of the Monte Carlo routine can be found in the ¹.

All measurements on tissue phantoms are performed with the probe shown in Fig. 1(a), while some of the biological tissue samples are measured with a probe geometry tip shown in Fig. 1(b). The probe is connected to an optical console which contains the light sources and the spectrometer. Fluorescence is excited from fiber A with light from a laser diode (Nichia NDU113E, $λx≈375nm$). The fluorescence light is collected by fiber B and sent to an optical spectrometer in the visible wavelength region (Andor DU420A-BR-DD). A filter in the spectrometer blocks light with a wavelength $λ<405nm$. For DRS measurements, an optical switch changes the light source of fiber A to a tungsten halogen broadband light source (Ocean Optics, HL-2000-HP). A second halogen light source is connected to fiber C, to allow the measurement of the DRS spectra at a second longer distance. The fluorescence and the DRS measurements use the same collection fiber B and the same spectrometer. Fiber D is connected to a second NIR spectrometer, which is not used for the fluorescence measurement, but extends the DRS spectra into the near-infrared region. The various light sources are switched on and off individually by computer-operated shutters. The setup is controlled by a laptop running software purpose-written in labview (from National Instruments, Austin, Texas). More details about the setup, including information about the calibration, can be found in ¹⁸.

During each measurement, three spectra are acquired: a fluorescence spectrum and diffuse reflectance spectra at the short fiber distance ($δ=0.3mm$) and the long fiber distance ($δ=1.8mm$). From the DRS measurement the absorption and the reduced scattering coefficients are determined for all wavelengths involved. This is done by fitting the spectra with the diffusion theory model of Refs. ¹³ and ¹⁹. The absorption coefficients $μam$ and $μax$ are determined from the extinction spectra of the individual chromophores using the concentrations as fitting parameters, while the reduced scattering coefficients $μsx′$ and $μsm′$ are approximated by a superposition of two exponential functions (one for Mie and one for Rayleigh scattering). A Levenberg-Marquardt nonlinear fitting algorithm is used to optimize the chromophores concentrations and scattering parameter. The algorithm is described and validated in detail in Refs. ¹⁸ and ²⁰. The algorithm generally yields very good agreement between measured and fitted DRS spectra. A drawback is that the model requires that the extinction spectra of all absorbers are known. Also, the algorithm becomes less reliable for short fiber distances ($ρ<1mm$) because the diffusion approximation breaks down. As a result, one typically obtains a different set of $μs′$ and $μa$ from the short distance (SD) DRS measurement than from the long distance (LD) with the long distance set being more trustworthy.

To determine the intrinsic fluorescence according to the MüllerZhang method, first the measured short distance DRS spectrum is fitted by the diffusion theory model as described above. From this fit the albedo $am$ can be determined. Furthermore, we can deduce the scattering in absence of absorption by setting the absorption coefficient from this fit equal to zero yielding $R0$. In the MüllerZhang method, the probe-specific parameters $S$ and $l$ are determined from measurements on known samples. From this $S$ we can then determine $β$ and subsequently $ε$. Inserting $S$, $l,ε$, and the measured $R$ and $Fxm$ in Eq. (10) yield the intrinsic fluorescence. In the MPM method, the parameters $S$ and $l$ are not a priori known. Here, we determine $κ$ using the relation $κ=R0(eβ−1)$. Inserting $am$, $R0$, and $κ$ in Eq. (6) and using for $R$ the measured short distance DRS spectrum we can deduce the parameter $β$ and thus the parameter $ϵ=(eβ−1)$. From Eq. (10) we can then deduce the intrinsic fluorescence from the measured fluorescence $Fxm$ apart from a scaling factor.

To determine the intrinsic fluorescence with the Monte Carlo look-up table (MCLUT) method, Eq. (15) is employed to calculate $ϕxmeff$, making use of the calculated previously look-up table for $kxm$. Then, Eq. (12) is used to determine the intrinsic fluorescence $fxm$, as defined above. This can be done with the $μs′−μa$-set from the short or the long distance DRS measurement, resulting in two different intrinsic fluorescence spectra, which we call short MCLUT and long MCLUT.

If the intrinsic fluorescence spectra $ηi$ of the individual fluorophores are known, then the (relative) concentrations $ci$ of the fluorophores can be determined by fitting $fxm$ according to Eq. (3) with the $ci$ as fit parameters. We are only interested in relative concentrations, so we set $ε$ and $ϕxm$ to 1 as scaling factors.

To validate the MCLUT method, phantoms with a single fluorophore are made by diluting latex balls (scatterer, diameter $diameter≈1μm$, Duke Scientific), red tissue marking dye (absorber, #63020-RD, Electron Microscopy Sciences), and a Lucifer Yellow preparation (fluorophore) in distilled water. These phantoms were also used to determine the constants $S$ and $l$ for the flat-head-probe according to ⁵.

To compare the different techniques of reconstructing the intrinsic fluorescence spectra in a quantitative way, a series of phantoms are prepared. They contain a mixture of two different fluorophores in various concentration ratios and an absorber in various concentrations. The diffuse reflectance and fluorescence spectra of the phantoms are measured and the intrinsic fluorescence spectrum of each phantom is recovered using all three techniques. From the recovered intrinsic fluorescence spectra the ratio of the two fluorophores in the phantom is determined. The deviation of the fluorophore ratio determined, thus from the true fluorophore ratio is taken as a measure for the quality of the recovery algorithm.

Nine phantoms were made, called N1 to N9. All phantoms contain the same amount of scatterers. The scatterers are Polybead® polystyrene microspheres from Polysciences, Inc. with a diameter of $1.025±0.01μm$. The initial suspension was diluted in demiwater to a polystyrene concentration of 2.4 V%. The absorbers used are a mixture of Amaranth and Tartrazine from SigmaAldrich. Amaranth and Tartrazine are both red food dyes with a high stability against photo-bleaching. Using a mixture creates a more characteristic absorption spectra than a single absorber. The extinction spectrum of the mixture was measured and it showed two distinct absorption maxima at 434 and at 522 nm. Phantoms N1, N2, and N3 are prepared with 5 V% of the absorber mix, phantoms N4, N5, and N6 with 1.25 V%, and phantoms N7 and N8 with 0.5 V%. Phantom N9 does not contain any absorber and is used as a reference.

The fluorophores used are preparations of Stilben-3 and brilliant Sulfaflavin in demiwater. Three different combinations of the two fluorophores were used. Both fluorophores have broad fluorescence peaks with a maximum around 425 nm (Stilben-3) and 525 nm (brilliant Sulfaflavin). These two peaks will typically still be visible in the intrinsic fluorescence spectra of a mixture, with the relative height of the peaks depending on the Stilben3/Sulfaflavin ratio. Phantoms N1, N4, and N7 were prepared with an (arbitrarily scaled) Stilben3/Sulfaflavin ratio of 1; N2, N5, and N8 with a Stilben3/Sulfaflavin ratio of 2.75; and N3 and N6 with Stilben3/Sulfaflavin ratio of 0.26. These ratios were determined by measuring the fluorescence spectra of the various mixtures and then fitting the spectra of the individual fluorophores. The intrinsic fluorescence spectrum of all fluorophore preparations are measured on the undiluted preparation without any additional absorbers or scatterer present.

To determine the applicability of the MPM and MCLUT techniques to biological tissue, we also measured various ex vivo human tissue samples. These measurements were conducted at The Netherlands Cancer Institute (NKI-AVL) under approval of the internal review board. Tissue samples from lung, liver, breast, and cervix are obtained from patients undergoing a surgical tumor resection. The freshly excised tissues are grossly inspected by a pathologist prior to our measurements. A different optical probe is used with a sharp step-like tip [see Fig. 1(b)] that can penetrate the tissue and measure inside the tissue volume. The measurements on cervical tissue are done on the surface of the ectocervix with the flat-head probe. The analysis is done as described above for the phantoms. No good values for the constants $l$ and $S$ are available for the step probe so the original MüllerZhang method is not used for the biological samples.

Figure 2 shows two typical examples of the look-up table $kxm$ calculated with our Monte Carlo routine. The images only show a two-dimensional subset where the values at the excitation wavelength are kept fixed at $μsx′=1.25mm−1$ and $μax=2.0mm−1$ in (a) and at $μsx′=5.0mm−1$ and $μax=1.5mm−1$ in (b) both for a source-detector fiber distance $ρ=0.3mm$. As one would expect, the correction factor $kxm$ decreases for increasing $μam$, i.e., the measured fluorescence intensity decreases with increased absorption at the emission wavelength. The influence of $μsm′$ is more complicated. For low values of $μam′$ there exists a maximum $kxm$, i.e., an increase of the scattering coefficient $μsm′$ increases the measured fluorescence intensity until a maximum is reached after which a further increase in $μsm′$ results in a sharp drop. The maximum grows less pronounced and finally disappears when the absorption coefficient $μam′$ increases further.

Figure 3(a) shows the dependence of $kxm$ when scattering and absorption properties are varied at the excitation wavelength for fixed $μsm′=10mm−1$, $μam=1mm−1$. The correction factor $kxm$ shows variations of the same order of magnitude as when varying the parameters at the detection wavelength.

For monochromatic excitation, there is only one value for $μsx′$ and $μax$ but a different ($μam′$, $μam$) duplet for every measured emission wavelength. Using those as input for the look-up table yields a correction factor spectrum $kxm(λm)$. Figure 3(b) depicts four such spectra, using the same ($μsm′$, $μam$) spectrum but with moderate differences in the $μsx′$, and $μax$ values. These differences cause significant deviations in the absolute values of the correction factor but also noticeable changes in the shape of the curves.

This section describes the recovery of the intrinsic fluorescence spectrum using the MCLUT method. As an example, we consider the phantom with one absorber, one fluorophore, and one type of scatters (see Sec. 2.5).

Figure 4 shows the DRS and fluorescence spectra measured on the single fluorophore phantom. The red paint has absorption peaks at 528 and 571 nm, which are clearly visible in the DRS spectrum and the latter also in the fluorescence spectrum. The absorption clearly has a major influence on the measured fluorescence spectra, which looks very different from the real intrinsic fluorescence of the fluorophore [Fig. 4(b)].

From the $μsx′$, $μax′$, $μsm$, and $μam$ from the short distance measurement the correction factor $kxm$ is determined by interpolation from the look-up table [Fig. 4(a)]. It is obvious that $kxm$ has similar features as the measured reflectance spectrum. Figure 4(b) shows the recovered intrinsic fluorescence spectrum using the MCLUT technique. Good agreement between the recovered intrinsic fluorescence spectra and the fluorescence spectra measured on the pure fluorophore preparation is obtained.

In this section, the recovery of the intrinsic fluorescence of the nine samples (N1 to N9) will be discussed. All samples are measured with the probe and setup described in Sec. 2.4, with DRS measurements being taken at both the short and the long distance. Figure 5 shows the DRS measurements and the resulting fits for two of the phantoms. Agreement between the measurement and the fits is limited by the uniform size of the scatterers, which causes ripples in the $μsm′$-spectra that cannot be reproduced by the exponential function used. Figure 6 depicts for the same two phantoms the measured fluorescence curves, the intrinsic fluorescence curve of the fluorophore mixture, and the intrinsic fluorescence curves recovered with the MüllerZhang, MPM, and MCLUT. The Monte Carlo look-up table method is used with both the short and the long DRS measurement.

Although not shown in this paper, the recovered intrinsic fluorescence of the other phantoms exhibit similar results: the strong absorption dips are removed from the fluorescence spectra for all the methods used. However, the relative heights of the two fluorescence peaks are rarely recovered correctly. Furthermore, the fluorescence intensity $>600nm$ is decreasing slower than observed in the real intrinsic fluorescence regardless of the method used. It is not clear whether this “tail” is caused by the algorithms or whether it is due to actual additional fluorescence (e.g., by the fiber probe or by the polystyrene microspheres).

To derive a measure for the quality of the fit, the reconstructed fluorescence curve is fitted as a superposition of the individual intrinsic fluorescence curves of Stilben 3 and brilliant Sulfaflavin. This fitting was done in the wavelength range from 405 to 575 nm to limit the influence of the low intensity, long wavelength “tail” on the result. From the fitting results the Stilben3/Sulfaflavin ratio is calculated. Table 1 lists the deviations (in %) from the expected Stilben3/Sulfaflavin ratio for each sample for the various methods used. The last two rows in Table 1 give the mean and the standard deviation for each technique.

In this section, different intrinsic fluorescence recovery methods shall be compared when applied to human breast, lung, liver, and cervix tissue samples. Typical examples of the DRS and fluorescence spectra obtained from these human tissue samples are shown in Fig. 7. To retrieve the scattering and absorption from the DRS spectrum, we apply the fitting as described above taking the chromophores hemoglobin, deoxyhemoglobin, water, lipid, beta-carotene, collagen, and bile (for liver samples) into account. In Fig. 7(a), 7(c), and 7(g) the double dip structure due to oxyhemoglobin near 550 nm is clearly visible and fitted well in these three cases. In Fig. 7(e), i.e., liver tissue, the double dip is less prominent and the blood is more deoxygenated. When we compare the recovered intrinsic fluorescence according to the MPM method with the MCLUT method, we observe that the intrinsic fluorescence reconstruction of the MPM and the short distance MCLUT methods are similar. In comparison, the long distance MCLUT methods seems to result in an “undercorrection,” meaning that the long distance MCLUT result is closer to the original measured spectra than those of the other methods. Apart from lung and liver where the absorption due to blood is large, the three models predict similar results for the recovered intrinsic fluorescence.

This paper reveals different ways to recover the intrinsic fluorescence spectrum from measured spectra of the fluorescence and the diffuse reflectance in turbid media. The MCLUT technique is, as far as the authors know, the first time that a look-up table has been used for this purpose.

The shape of the look-up table calculated by the Monte Carlo routine is complex but plausible: with increasing absorption, both at the excitation and the emission wavelength, the measured fluorescence decreases [Figs. 2 and 3(a)]. The influence of scattering is more complicated. The scattering coefficient at the emission wavelengths inversely correlates with the maximum distance from the fluorophore that, on average, the emitted photon will reach on its random walk through the tissue. Therefore, increasing $μsm′$ at first raises the measured fluorescence as the photon density at the collection fiber increases but then decreases because most photons cannot reach the collection fiber anymore. This explains the maximum in Fig. 2(a) and 2(b). Varying the scattering at the excitation wavelength also produces a maximum in the fluorescence intensity for medium $μsx′$. For low scattering, the excitation photons can enter deep into the tissue before being absorbed by a fluorophore, resulting in very few fluorescent photons exiting the tissue. For large $μsx′$, the excitation photons will be absorbed very close to the delivery fiber and many photons will be scattered out of the tissue without being absorbed. Therefore, maximum collection of fluorescent photons occurs for medium scattering at the excitation wavelength.

This complex dependency of $kxm$ on the optical properties at both the emission and the excitation wavelength justifies using a four-dimensional look-up table. In particular, the influence of $μsx′$ and $μax$ cannot be described by a simple scaling parameter. This is demonstrated in the examples in Fig. 3(b) where a comparison of, e.g., the solid red curve with the dash-dotted green curve shows that the shape of a correction spectrum can vary significantly if the optical properties at the excitation wavelengths change. This four-dimensionality is an important difference between the look-up table methods and the photon-migration based methods, since the latter does not consider the optical properties at the excitation wavelength except for using $μsx$ as a scaling factor.

On the tissue phantoms, the various techniques to reconstruct the intrinsic fluorescence all perform reasonably well. The strong absorption dips are filtered out regardless of the method used. However, the various methods produce subtly different shapes for the intrinsic fluorescence spectra that affect the determined relative fluorophore concentrations. These determined ratios of fluorophore concentrations were used as a measure of the goodness of intrinsic fluorescence reconstruction.

Table 1 shows that none of the methods works best for all samples. To get a figure of merit for the overall quality of each technique, the relative deviations are squared and summed up (see the row labeled $Σ2$ in Table 1). By this measure, the MPM and the MCLUT methods using the data from the long distance measurement yield the best results. The MCLUT method making use of the short distance diffuse reflectance measurement and the MüllerZhang method perform somewhat less well but the differences are marginal.

A common issue is caused by the fact that the source-detector fiber distance is of the same order of magnitude as the typical scattering length. This means that the light transport in tissue does not comply with the criteria necessary for the diffusion approximation to hold. In all the methods investigated in this paper, the step to retrieve the scattering and absorption parameter rely on this diffusion theory. In the MüllerZhang method and the MPM this is partly solved by using only the scattering information (i.e., the function $R0$) determined by this diffusion theory fit while the absorption is only taken into account by using the separately measured diffuse reflectance spectrum $Rm$. For the MCLUT methods, however, this is a serious issue when the short distance DRS measurement is used. The MCLUT using the long distance DRS measurement is affected least. The problem with the diffusion approximation for short distances therefore explains why the long MCLUT works slightly better than the short MCLUT for the phantom measurements. These fitting errors might be lessened by employing, instead of the diffusion approximation based fitting of the DRS spectrum, an inverse Monte Carlo based fitting of the DRS spectrum.

The slightly improved performance of the MPM compared to the original MüllerZhang method might be explained by the additional freedom in the MPM method, where the parameters $l$ and $S$ are determined in each fit separately, compared to the MüllerZhang method, where they are determined once based on various test phantoms.

Care must be taken with the ranking of the various methods to recover the intrinsic fluorescence used in the report. As measurements were only done on a small number of similar phantoms, there remains a large margin for statistical and systematic errors. For example, it appears as if all reconstruction methods produce increasing errors for phantoms with decreasing Stilben3/Sulfaflavin ratios, but it is unclear whether it is a real effect or a statistical artifact and, if it is a real effect, whether it is caused by our special choice of tissue phantom. Further investigation on this point is required. Also, the method chosen to compare the various techniques is rather arbitrary, even though recovering the intrinsic fluorescence spectra and determining the fluorophore concentrations from it is of large practical relevance. Still, we believe that the phantom measurements give a good first impression of the performance one can expect from the various intrinsic fluorescence recovery methods. To the best of our knowledge, this is the first attempt of a quantitative comparison of different techniques. In general, we can conclude that the various methods are able to reconstruct the intrinsic fluorescence but the error in the reconstruction can lead to errors in relative concentration on average of up to 25% (i.e., one standard deviation) as found in our examples.

From the examples presented on human breast, lung, liver, and cervix tissue we observe that the MPM and short distance MCLUT methods produce similar results for the recovered intrinsic fluorescence. For breast and cervix tissue these agree with the results from the long distance MCLUT method. Although no benchmark with the real intrinsic fluorescence can be made, these results are promising since the results are obtained with two highly independent reconstruction methods. When we compare the DRS spectra measured for this human tissue with the phantom studied, we observe that the amount of absorption and scattering present is of the same ballpark figure. Therefore, we expect that the recovery of the intrinsic fluorescence also has an accuracy similar to that of the recovery of the tissue phantoms. However, there may be additional unknown chromophores present in the biological tissues which may cause additional errors in the DRS fit and decrease the overall accuracy.

The long distance MCLUT method produced rather different reconstructions in tissues with strong absorption bands [Fig. 7(c) and 7(e)]. One also sees that traces of the blood absorption bands remain in the recovered intrinsic fluorescence spectra [Fig. 7(d) and 7(f)]. We attribute this to a low fidelity in the reconstruction of the $μa$ and $μs′$ values due to strong tissue inhomogeneities, high absorption, and a complicated probe geometry. In order for the MCLUT method to work, the optical properties in the volumes from the fluorescence and the DRS measurements have to be similar. This was the case for the phantoms, but biological tissues are much less homogenous. This is mainly an issue when using the long distance DRS measurement since it measures a much larger volume with less overlap with fluorescence volume. Also, for samples with high blood content (liver, lung), the long distance of the probe used was too long, resulting in excessive absorption. In Fig. 7(c) and 7(d), the intensity of the long distance DRS measurements in the absorption bands is practically zero. With the most characteristic part of the spectrum gone, fitting errors increase and the long-distance DRS fit becomes imprecise. This may attribute to the underfitting seen in Fig. 7(b) and 7(c) in the short wavelength regions. Last but not least, the complicated geometry of the fiber optical probe which was necessary to penetrate into tissue [Fig. 1(b)] might have added additional distortions to the long-distance DRS spectrum (even though the fiber distances are the same for both types of probe). For these reasons, we believe that for liver and long tissues the long-distance MCLUT results are less accurate than the results from the short-distance MCLUT and from MPM.

In general, the MCLUT method depends more on the quality of the DRS fit than the MPM method. This is because $μa$ and $μs′$ are more sensitive to fitting errors than $R0$. This reduced robustness may offset the advantage of not having to measure fluorescence and DRS in the same geometry. Since the computational effort is preloaded in calculating the look-up table, the recovery itself is fast and computationally cheap. Overall, the results show that it is possible to accurately recover the intrinsic fluorescence spectrum from a short-distance fluorescence measurement using the MCLUT method with only a long-distance DRS measurement. Results are best when the sample is homogenous on the scale of the long distance and has no absorption bands strong enough to introduce errors in the DRS fit.

A new method based on Monte Carlo look-up-tables to recover intrinsic fluorescence from measured fluorescence on turbid media has been presented. This new method requires an additional diffuse reflectance measurement taken at the same location as the fluorescence measurement but not necessarily with the same probe geometry. The method has been validated on various phantoms with known intrinsic fluorescence. The MCLUT was also benchmarked against the photon-migration method for intrinsic fluorescence recovery developed by Müller et al. and Zhang et al. (MüllerZhang method) and with a slightly modified photon-migration method that does not need prior knowledge of optical probe specific parameters. These methods require a diffuse reflectance spectrum measurement taken at the same location and with the same optical probe geometry as the fluorescence measurement. The MCLUT showed similar reconstruction accuracy as the methods based on the photon-migration method. Furthermore, also the MCLUT with different probe geometry for the DRS measurement resulted in similar intrinsic fluorescence reconstruction as the photon-migration method as long as the source-detector fiber distance is chosen such that the DRS spectrum in the relevant wavelength region is well detectable while, furthermore, the fiber distance is smaller than the average scale of the tissue inhomogenities. The modified photon-migration based method proposed in this paper resulted in a slightly improved result than found with the original MüllerZhang method.

In the current MCLUT method, the absorption and scattering parameters are determined from the DRS measurement by employing a diffusion theory model based fitting. The MCLUT method is therefore sensitive for errors in both absorption and scattering coefficient. The photon-migration based model uses only the scattering coefficient from this fitting step, making this method less sensitive for fitting errors than the MCLUT in the current form. In general, we can conclude that the MCLUT is a promising technique to recover the intrinsic fluorescence of biological tissues. It combines more flexibility in the probe design and fast reconstruction times with similar accuracy in reconstruction of the intrinsic fluorescence as found in the photon-migration based reconstruction methods.

We thank Prof. S. L. Jacques (Oregon Medical Laser Center) for providing the mfluor.c code to benchmark our Monte Carlo program code. Prof. Theo Ruers (Netherlands Cancer Institute) is acknowledged for enabling the measurements on the human tissue samples. We thank Rami Nachabé and Arnold van Keersop for their support in writing the MPM code, Walter Bierhoff for providing the optical probes, and Jeroen Horikx for providing the optical console.

To validate our Monte Carlo program, we benchmarked our Monte Carlo routines with existing peer-reviewed Monte Carlo routines. The basic Monte Carlo routine and the calculation of $A(r,z)$ were validated by comparing selected outputs with results produced by the program MCML by Wang et al.¹⁴ Both programs were in excellent agreement.

To validate the full fluorescence routine (called fluostitch.m), the program mcfluor.c from Prof. S. L. Jacques was employed.¹⁷Figure 8(a) shows the distribution of the fluorescence intensity $F(ρ)$ calculated both with mcfluor.c and with fluostitch.m for three different tissue parameter sets. To obtain the same accuracy required longer computation times for mclfuor.c program than our new code fluostitch.c. Figure 8(b) shows the results of another comparison between the two programs. Instead of showing the whole spectrum $F(ρ)$, only the values at two distances are compared. The tissue properties at the excitation wavelength are kept constant but are varied at the emission wavelength. The graph shows again, a good agreement. The deviations are of the same order of magnitude as the noise present in the data points.