1. Introduction

Reflectance spectroscopy is widely used for noninvasive measurement of the tissue absorption and scattering coefficients, which can provide diagnostic information about tissue function and structure, respectively [1–9]. Scattering in tissue originates from spatial heterogeneities of the optical refractive index that occur on size scales ranging from a few nanometers to a few millimetres [9]. Scattering properties of tissue include the scattering coefficient (μ_s), which describes the mean free path between scattering events, and the scattering phase function (p(θ)), which describes the angular distribution of scattering events. These properties are often described in terms of the scattering anisotropy g, which is the expectation value for the cosine of the scattering angle g = < cos(θ) >, and the reduced scattering coefficient ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{\mu }_s\left(1\hspace{0.17em}-\hspace{0.17em}g\right)$, which is a description of the combined effect of scattering coefficient and average scattering angle. Reflectance spectroscopy measurements with large source-detector separations (≈ 10 mean free paths [10]) can be modeled using diffusion theory, which utilizes ${\mu }_s^{\prime }$ to characterize the effect of scattering on reflectance. Diffuse measurements return information about bulk tissue properties, which is a limitation because many clinical diagnostic applications require more localized measurements [11–13]. Quantitative extraction of the absorption and reduced scattering coefficients from reflectance measurements with small source-detector separations is complicated because light transport in this non-diffuse regime, where collected photons may undergo few scattering events, is no longer dependent exclusively on ${\mu }_s^{\prime }$, but is also sensitive to large angle back scattering events, a characteristic of the scattering phase function.

The influence of phase function on reflectance intensity measured close to the source has been described in the literature. For example, Mourant et al. [14] reported significant differences in the Monte Carlo simulated reflectance intensities for a small source-detector probe for different phase functions (Henyey-Greenstein vs. Mie) but with the same anisotropy. Kienle et al. [15] also reported large deviations in the calculation of μ_a and ${\mu }_s^{\prime }$ after simulating the reflectances for a source-detector separation of < 1 mm for multiple phase functions by applying a standard solution for the diffusion equation [10] to extract the optical properties. Bevilacqua and Depeursinge [16, 17] systematically investigated the influence of higher order moments of the phase functions for several source-detector separations on the reflectance intensities. They introduced a phase function dependent coefficient γ which, together with μ_a and ${\mu }_s^{\prime }$ and refractive index n, would give a complete description of the reflectance for systems with measurements from multiple short source detector separations. Hull and Foster later developed an analytical description of light transport termed the P3 approximation, which models both diffuse and non-diffuse light scattering by including the first three moments of the phase function [18]; this approach has been shown to be valid for source-detector separations as low as ≈ 0.5 mm [13], but not directly applicable for devices with overlapping source and detector locations. Other empirical investigations of the relationship between reflectance intensity with small source-detector separations and ${\mu }_s^{\prime }$ [20–22] did not fully characterize the effect of phase function on the collected signal.

Recently, our group reported on the use of experimentally validated Monte Carlo (MC) simulations to investigate the effect that optical properties have on the propagation of photons that are collected during a single fiber reflectance (SFR) measurement [23]. SFR spectroscopy uses a single fiber for the delivery and detection of light; a schematic of the device setup is shown in Fig. 1. The advantages of an SFR device are the simplicity of the setup machinery and the small dimensions of the fiber-optic probes, which allows e.g. measurement through the lumen of FNA-needles [24, 25]. In the past, single fiber geometries have been used for the investigation of fluorescence in tissue and tissue simulation phantoms [26–29] and for particle size analysis [30, 31]. The relationship between single fiber reflectance intensity and optical properties was studied by Moffit and Prahl [19]; however, their analysis did not yield quantitative estimates of ${\mu }_s^{\prime }$ or characterize phase-function dependence of the signal. Further studies have focused on technical aspects of single fiber reflectance, such as the need to remove internal specular reflections from the collected signal [32] and the dependence of the collection efficiency of a single fiber on the optical properties of the sampled medium [33]. Our previous work has focused on utilizing SFR spectroscopy to quantitate the absorption coefficient (μ_a), and the corresponding chromophore concentrations, in turbid media. These studies introduced an empirical function that describes the dependence of the effective photon path length [23,34] on the optical properties of the sampled medium. A limitation of our method is that the path length relation requires knowledge of ${\mu }_s^{\prime }$, a parameter that is currently not quantitatively extracted from the measured spectra. The current spectral analysis algorithm uses estimated values for the wavelength-dependent ${\mu }_s^{\prime }$ and the uncertainty introduced has been documented [24]. A quantitative measurement of ${\mu }_s^{\prime }$ would improve the accuracy of our chromophore concentration estimates. Moreover, quantitative measurement of ${\mu }_s^{\prime }$ using an SFR device is inherently useful, because it may provide information about the tissue microarchitecture [9]; information that is complementary to the physiological information extracted from μ_a, and may have diagnostic value. Additionally, quantitative knowledge of ${\mu }_s^{\prime }$ and μ_a extracted from SFR measurements may be used to correct for the effect of optical properties on single fiber fluorescence measurements [26–29].

 

Fig. 1 Schematic of single fiber reflectance probe machinery.

Download Full Size | PDF

There are two challenges associated with quantitative extraction of ${\mu }_s^{\prime }$ from an SFR measurement. First, a calibration procedure must be designed that allows absolute measurement of the reflectance (R_SF) of a turbid medium from an SFR measurement. Second, R_SF must be related to ${\mu }_s^{\prime }$. The first challenge will be addressed by designing a calibration phantom that returns a known percentage of incident light during measurement. This measurement allows the scaling of reflectance intensity measured on a turbid medium to a percentage of the absolute amount of incident measurement light. Although the concept of using calibration phantoms to scale the reflectance is not novel, for an SFR device (or any other small source-detector separation device) the scattering phase function of the calibration phantom is a critical factor that needs to be accounted for in the design of the phantom; the present study considers phase function effects in detail. The second challenge is addressed by investigating the relationship between R_SF and the scattering properties of the medium, and the influence of scattering phase function using an extensive set of Monte Carlo simulations. The simulations are experimentally validated using highly scattering ( ${\mu }_s^{\prime }\hspace{0.17em}\approx \hspace{0.17em}20\hspace{0.17em}-\hspace{0.17em}200\hspace{0.17em}{\text{mm}}^{-1}$) solid phantoms consisting of titanium dioxide (TiO₂) suspended in silicone, as well as low and medium scattering (0.05 < ${\mu }_s^{\prime }$ < 8 mm⁻¹) liquid phantoms consisting of a range of dilutions of Intralipid. A semi-empirical model is introduced that describes the relation between R_SF and dimensionless scattering ( ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$) for a wide range of phase functions.

Our methodology to experimentally extract ${\mu }_s^{\prime }$ of a turbid medium from SFR spectra, using a calibration phantom and a mathematical model of $R_{\mathit{\text{SF}}}\hspace{0.17em}=\hspace{0.17em}f\left({\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\right)$, is validated using SFR measurements on Intralipid phantoms. We show that with knowledge of the wavelength-dependent phase function, the model function can be used to accurately extract ${\mu }_s^{\prime }$ from SFR measurements for a wide range of reduced scattering coefficients ( ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.05\hspace{0.17em}-\hspace{0.17em}8\right]\hspace{0.17em}{\text{mm}}^{-1}$) and fiber diameters (d_fib = [0.2, 0.4, 0.6, 0.8, 1] mm).

2. Materials and methods

2.1. Monte Carlo simulation

The Monte Carlo (MC) code used to simulate single fiber reflectance spectroscopy in this study has been described in detail previously [23]. In short, our MC code is based on the MCML code [35] that stochastically simulates photon propagation within a turbid medium. During photon propagation, each photon step size was selected from an exponential distribution weighted by the scattering coefficient, and each scattering angle was selected from a user-specified phase function. Reflection and refraction due to the index of refraction mismatch at the medium/fiber and medium/air interface were calculated using the Fresnel equations and Snells law. The index of refraction (n) of the medium and fiber were specified at 1.37 and 1.5, respectively. The numerical aperture (NA) of the fiber was set as 0.22. Photons were initialized by selecting a location from a uniform distribution on the single fiber face in contact with the turbid medium (z = 0), and the launch direction was selected from a uniform distribution of angles within the fiber cone of acceptance, where the acceptance angle was given as ${\Theta }_a\hspace{0.17em}=\hspace{0.17em}\mathit{\text{asin}}\hspace{0.17em}\left(\frac{\mathit{\text{NA}}}{n_{\mathit{\text{medium}}}}\right)$. For photons propagating within the turbid medium that cross the medium interfacial boundary at z = 0, it was checked if they hit the fiber face; those in contact with the fiber face and traveling at an angle within the fiber cone of acceptance were collected, the rest were terminated. Photons propagating within the medium far from the fiber face do not contribute to the collected reflectance intensity and were terminated at a hemispherical limit from the fiber face of $10\frac{d_{\mathit{\text{fib}}}}{{\mu }_s^{\prime }}$; a limit that was confirmed to not influence model outputs for the range of optical properties investigated in this study. For each simulation, the MC code outputted the reflectance intensity as the ratio of the total number of photons collected (TPC) to the total number of photons launched (TPL) during the MC simulation, calculated as:

 

Fig. 2 Angular distribution of scattering events for selected scattering phase functions.

Download Full Size | PDF

2.2. Semi-empirical model of the single fiber reflectance intensity

Inspection of simulated reflectance data over the range of simulated scattering properties yielded an observable relationship between $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ and the product ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$; this product is referred to throughout this paper as dimensionless scattering. A semi-empirical model is utilize to describe this relationship, given by the following set of Equations,

The parameters [ρ ₁, ρ ₂, ρ ₃] were fitted by minimizing the weighted residual error between the simulated ( $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$) and model-estimated ( $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$) reflectance intensities, with each point weighted by the inverse of the simulated data point. It is important to note that neither Φ nor η_c are observable from our experiments, and while they were estimated by the MC model, the fitted parameter values [ρ ₁, ρ ₂, ρ ₃] were only fitted to reflectance intensity data using Eq. (2); i.e. we did not separately fit Eqs. (3) and (4) which represent unobservable quantities. Parameters were estimated using a Levenberg-Marquardt algorithm [38] written into Matlab code (version R2009a, MathWorks). Confidence intervals of the parameter estimates were calculated from the square root of the diagonal of the covariance matrix, as described in detail elsewhere [39].

2.3. Experimental setup for single fiber reflectance spectroscopy

Figure 1 shows a schematic of the single fiber reflectance spectroscopy device setup, which consists of a halogen light source (HL-2000-FHSA, Ocean Optics, Duiven, Netherlands), a spectrometer (SD2000, Ocean Optics, Duiven, Netherlands) and a solid core single fiber probe. The single fiber is connected to the light source and the spectrometer via a bifurcated fiber. During measurements, light travels from the light source through the fiber and exit into the probe tip that is in contact with the sampled medium. Photons that backscatter to the fiber face at an incident angle that is within the cone of acceptance enter the fiber core and travel to the spectrometer. To remove specular reflections within the collected reflectance intensity due to index of refraction mismatch at the probe/medium interface tip, the probe tip is polished at an angle of 15 degrees [34]. Variations in the output of the lamp, transmission characteristics of the fiber, and sensitivity of the spectrometer, as well as remaining specular reflectance and other internal reflections are being accounted for by performing a calibration procedure that we have described previously [34]. The calibration includes measurements of white and black spectralon standards (Labsphere SRS-99 and SRS-02, spectrally flat, with 99% and 2% reflectance, respectively) and a measurement of water in a black container. The single fiber reflectance signal can be calculated as

2.4. Solid phantom preparation

Solid phantoms were prepared according to the recipe reported by de Bruin et al. [40]. Silicone [Sylgard 184 Silicone Elastomer Kit (Dow Corning Europe SA, Seneffe, Belgium)] served as a stable matrix and Titanium dioxide (TiO₂) (Sigma Aldrich, Zwijndrecht, Netherlands) was used as a scatterer. All solid phantoms were prepared with a weight of 20 g, including different amounts of TiO₂ to create a series of different highly scattering phantoms (TiO₂ weight percentage= [0, 2.5, 5, 10, 20, 30] %). The silicone kit consisted of the silicone compound and a curing agent which was always mixed in a ratio of 9:1. To prepare the phantoms the silicone compound and TiO₂ were weighed, mixed with a Dremel 300 (Dremel, Leinfelden-Echterdingen, Germany) and placed into an ultrasound bath for 10 min, to assure that the TiO₂ particles were homogeneously distributed within the silicone. Then the curing agent was added and the composition was mixed again until the agent was completely mixed into the compound. The mass was poured into the mold and placed in a vacuum desiccator (Duran Group GmbH, Mainz, Germany) and held under vacuum until remaining air bubbles were removed from the phantom (1 – 3 hours). The phantoms were cured after approximately 12 hours at room temperature.

During the single fiber reflectance measurements the single fiber probe tip was brought into contact with the phantom surface, with no gap between fiber face and phantom surface. Three spectra were recorded at 3 different locations on the phantom surface and averaged. Each phantom was measured with the 5 different fiber probes with diameters of [0.2, 0.4, 0.6, 0.8, 1.0] mm. To test the homogeneity of the phantoms and the reproducibility of the preparation, 6 solid phantoms with a weight percentage of 30 % TiO₂ were prepared and single fiber spectra were measured at 10 different locations on each phantom surface. The maximum relative deviation between measurements on two phantoms was 0.66 %, while the intra-phantom variability was even lower.

2.5. Liquid phantom preparation

The liquid phantoms were prepared by mixing 0.9 % NaCl (Baxter, Utrecht, Netherlands) with different amounts of Intralipid 20% (Fresenius Kabi, s-Hertogenbosch, Netherlands). The reduced scattering coefficient of the Intralipid stock solution was determined to be ${\mu }_s^{\prime }\left(800\text{nm}\right)\hspace{0.17em}=\hspace{0.17em}18.1\hspace{0.17em}{\text{mm}}^{-1}$ by using the formula given by Michels et al. [36]. The dependence of ${\mu }_s^{\prime }$ on the intralipid concentration was assumed to be linear and phantoms with ${\mu }_s^{\prime }\left(800\text{nm}\right)\hspace{0.17em}=\hspace{0.17em}\left[0.072,\hspace{0.17em}0.144,\hspace{0.17em}0.24,\hspace{0.17em}0.36,\hspace{0.17em}0.48,\hspace{0.17em}0.72,\hspace{0.17em}0.96,\hspace{0.17em}1.2,\hspace{0.17em}2.4,\hspace{0.17em}3.6\right]\hspace{0.17em}{\text{mm}}^{-1}$ were prepared to a volume of 20 ml. Reflectance measurements were performed by submerging the probe tip a few millimeters into the phantom and suspending it multiple centimeters from the container bottom or sides, such that the probe would not collect reflections from container surfaces. The measurement was repeated 3 times for each phantom, with these spectra averaged for each phantom. All measurements were performed with 5 fiber probes varying in diameter: [0.2, 0.4, 0.6, 0.8, 1] mm. The phantom was gently agitated before collecting each spectrum to ensure homogeneity.

2.6. Absolute calibration of experimental single fiber reflectance spectra

Measurement of single fiber reflectance intensity from the set of solid TiO₂ phantoms (outlined in Section 2.4) yielded $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{rel}}}$ vs. TiO₂ weight percentage (wt %) for a range of TiO₂ wt % and fiber diameters. TiO₂ wt % was converted to ${\mu }_s^{\prime }$ by estimating a linear conversion factor. In order to do this, measurements on the set of solid phantoms from all fibers were interpolated onto a curve that was normalized to the highest measured intensity in the set; this conversion corrected for different probe-spectralon distances, as noted in Section 2.3. Then a linear conversion factor ( ${\mu }_s^{\prime }\hspace{0.17em}/\hspace{0.17em}{\text{TiO}}_2\hspace{0.17em}\text{wt}\hspace{0.17em}\%$) was estimated by scaling the normalized reflectance vs. TiO₂ wt % data to fit the reflectance predicted by the mathematical model, presented in Eq. (2). This calculation yields a ${\mu }_s^{\prime }$ for each TiO₂ phantom.

Measurements of liquid Intralipid phantoms were converted from a value relative to spec-tralon distance, $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{rel}}}\hspace{0.17em}\left(\mathit{\text{IL}}\right)$, to an absolute percentage scale, $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}\hspace{0.17em}\left(\mathit{\text{IL}}\right)$, by,

3. Results

3.1. Mathematical characterization of $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ dependence on scattering properties

3.1.1. Influence of phase function on $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ vs. dimensionless scattering

Figure 3 shows the MC simulated reflectance intensity ( $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$) vs. dimensionless scattering ( ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$) with combinations of a wide range of ${\mu }_s^{\prime }\hspace{0.17em}\left[0.1\hspace{0.17em}-\hspace{0.17em}300\right]\hspace{0.17em}{\text{mm}}^{-1}$ and d_fib [0.2 – 1.0] mm. Here, $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ is expressed as a percentage of the incident photons. In Fig. 3, symbols differentiate the specified phase function, including $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ data from the HG (g values in the range of [0.5 – 0.9]) and MHG (effective g = 0.9) phase functions (see legend). Inspection of the data for each investigated phase function, shows that $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ increases monotonically vs. ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ until approaching an asymptotic limit for high dimensionless scattering values. For all phase functions investigated, the $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ data show a transition from a phase-function dependent regime at low ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ to a region insensitive to phase function at higher ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$, with the transition point between the two regimes occurring for all phase functions at approximately ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}>\hspace{0.17em}10$. In the low dimensionless scattering region, there is a clear stratification between reflectance intensities measured from different phase functions, with increased reflectance attributable to an increase in the backscattering component of the phase function. Therefore, for the same effective g, $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ is greater for the MHG than for the HG phase function; the difference is attributable to more high angle backscatter events in the MHG phase function, which is evident from the plot of the phase functions displayed in Fig. 2. Moreover, the $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ from HG simulations stratifies in order of decreasing backscatter component, with greater reflectance collected from phase functions containing a higher percentage of large angle backscatter events (e.g. g = 0.5 resulted in a higher reflectance than g = 0.9). Conversely, in the region of high dimensionless scattering ( ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}>\hspace{0.17em}10$) the $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ shows an insensitivity to backscattering features of the phase function, as data from all phase functions converge until they overlap and then asymptotically approach a limiting reflectance value. These results indicate that $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ is dependent on the phase function in a region of low-mid dimensionless scattering but approaches phase function independent behavior in a region of very high values of ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$.

 

Fig. 3 Single fiber reflectance vs dimensionless scattering for HG phase function (g=[0.5,0.7,0.8,0.9] and MHG phase function (g=[0.9]).

Download Full Size | PDF

3.1.2. Mathematical model of $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ vs. dimensionless scattering

The observed relationship between $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ and ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ was described using the mathematical model introduced in Section 2.2. Figure 4(A) shows $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ for the HG phase function with g = 0.8 for ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.1\hspace{0.17em}-\hspace{0.17em}300\right]\hspace{0.17em}{\text{mm}}^{-1}$ and d_fib = [0.2 – 1.0] mm (blue cross marks), and the corresponding model predicted reflectance ( $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$) estimated from Eq. (2) (black line). The estimated parameter values, given in Table 1, result in predictions of $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ that are highly correlated with simulated reflectance intensities over the full dimensionless scattering range, as evidenced by a Pearson Correlation Coefficient of r = 0.999 and a mean residual error of < 3%.

 

Fig. 4 (A) Single fiber reflectance and mathematical model fit for HG phase function (g=0.8). Model estimates of incident photons contacting fiber face (B) and collection efficiency of fiber (C).

Download Full Size | PDF

Table 1. Estimated Parameter Values Resulting from Fits of Single Fiber Reflectance Model (Eqs. (2)–(4)) to R_SF Simulated by Monte Carlo Models for HG and MHG Phase Function (top) and Wavelength-Dependent Phase Functions [36] in Intralipid (bottom)

View Table

Equation (2) expresses $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ as the product of Φ, the fraction of incident photons that are remitted from the medium in contact with the probe face, and η_c, the fraction of photons that contact the fiber face within the cone of acceptance and are collected. Figures 4(B) and 4(C) show the MC model estimates (green squares and red diamonds) and mathematical model predictions (black lines) of Φ and η_c, respectively, expressed as percentages, as calculated from Eqs. (3) and (4). Figure 4(B) shows that Φ increases with ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ and exhibits a saturating behavior as it approaches the 100% upper limit, a feature well-described by the model. Figure 4(C) shows that η_c shows a decreasing trend vs. dimensionless scattering, until approximately ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}<\hspace{0.17em}1$, at which point it approximates a limit linked to the fiber opening angle of acceptance; this behavior is consistent with a previous theoretical analysis [33] and is discussed in detail in Section 4. It is important to note that the estimated model parameters in Eq. (2) were not fit to MC estimates of Φ and η_c, as they are not observable experimentally. Instead Φ and η_c were calculated from the parameters resulting from the fit of $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ from Eq. (2) to the data, as shown in Fig. 4(A). The results presented in Figs. 4(B) and 4(C) suggest that the mathematical expressions used to represent Φ and η_c as components of $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ are, in fact, representative of the true underlying factors contributing to collected reflectance.

Figure 5 shows $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ vs. $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ for ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.1\hspace{0.17em}-\hspace{0.17em}300\right]\hspace{0.17em}{\text{mm}}^{-1}$ and d_fib = [0.2 – 1.0] mm, with symbols distinguishing each phase function, including HG (g = [0.5 – 0.9]) and MHG (g = [0.9]). Model parameters for $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ were fit for each simulated phase function; estimated parameter values for each fit are shown in Table 1 (top). Model predictions show an excellent agreement with simulated reflectance intensities over the range of phase functions investigated, with all r values above 0.998; additionally, the mean residual between all model estimated and MC simulated data points was < 3%. These data indicate that the dependence of $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ on dimensionless scattering can be accurately estimated for a given phase function.

 

Fig. 5 MC simulated single fiber reflectance vs. model estimate.

Download Full Size | PDF

3.2. Experimental determination of ${\mu }_s^{\prime }$ from R_SF

3.2.1. $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ measured from solid TiO₂ phantoms

Reflectance measurements on TiO₂ phantoms yielded $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{rel}}}$ vs. TiO₂ wt % data in the range [2.5 – 30] wt %. As outlined in Section 2.6, these data were scaled to $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ using the semi-empirical reflectance model. Because the exact phase function of the TiO₂ phantoms is unknown, we used the averages of the fitted parameter values: ρ ₁ = 7.9, ρ ₂ = 1.08, ρ ₃ = 2.4. Using these parameter values, a linear conversion factor of ( ${\mu }_s^{\prime }\hspace{0.17em}/\hspace{0.17em}{\text{TiO}}_2\hspace{0.17em}\text{wt}\hspace{0.17em}\%$) = 7.42 ± 0.35 was found. Figure 6 shows the $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ (800nm) from all measurements on TiO₂ phantoms vs. ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ as measured with 5 probes with diameters in the range [0.2 – 1.0] mm, using the 7.42 ( ${\mu }_s^{\prime }\hspace{0.17em}/\hspace{0.17em}{\text{TiO}}_2\hspace{0.17em}\text{wt}\hspace{0.17em}\%$) ratio. Also displayed on the plot are $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ vs ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ from simulations of an HG phase function with g = 0.8, which shows good agreement between experimental and simulated data. Note that due to the high reduced scattering coefficients of the phantoms, these measurements are in the high dimensionless scattering region ( ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}>\hspace{0.17em}10$), which showed an insensitivity to phase function. Varying the parameters ρ ₁ – ρ ₃ in Eq. (2) to be representative of HG (g = [0.5 – 0.9]) and MHG (g = [0.9]) phase functions (see Table 1) induced only a 2.2% mean residual error between $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ (based on the average ρ ₁ – ρ ₃) and $R_{\mathit{\text{SF}}}^{\mathit{\text{Model}}}$ (based on the different ρ ₁ – ρ ₃ combinations listed in Table 1). These results led to the selection of the 30% TiO₂ phantom (which equated to ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}222.6\hspace{0.17em}\pm \hspace{0.17em}10.5\hspace{0.17em}{\text{mm}}^{-1}$) as the calibration phantom for subsequent measurements on Intralipid.

 

Fig. 6 Single fiber reflectance measured on TiO₂ phantoms and simulated by MC model (HG phase function with g=0.8).

Download Full Size | PDF

3.2.2. $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ measured from Intralipid phantoms

Figure 7(A) shows $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ from measurements of Intralipid at multiple wavelengths (λ = [400 – 900] nm in steps of 100 nm) vs. dimensionless scattering. Here the $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ values represent the raw intensities scaled by measurements on the solid calibration phantom, allowing the conversion to an absolute scale. This calibration technique allows reflectance intensities from measurements with 5 fibers with different diameters to fall onto the same curve. Figure 7(B) shows $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ vs. ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ from MC simulations of the reflectance measurements in Intralipid over the same range of wavelengths. Here, the phase function for each wavelength was specified using an empirical expression (as in [36]), and the corresponding anisotropy values were in the range g = [0.477 – 0.818]. For both measurements and simulations, the R_SF values measured at different wavelengths (and therefore different phase functions) show stratification at low dimensionless scattering values, and gradually collapse towards a limiting value; these features are similar to the data shown in Fig. 3 for HG and MHG phase functions. Figure 8(A) displays the corresponding simulated $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ vs. measured $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ values. These data are highly correlated (r = 0.998) and show good agreement (with a mean residual of < 8%) over a range of measured reflectance values that span 2 orders of magnitude (range = [0.0158 – 1.57]%). It should be noted that relative deviation between simulation and experimental data increases for decreasing measured reflectance, and in turn for decreasing dimensionless scattering; this is discussed in Section 4. After accounting for the index of refraction mismatch between experimental and simulated phantoms (1.33 vs. 1.37, respectively), a linear fit of these data through the origin gives a slope of 1.04; potential reasons for the slight offset are given in Section 4.

 

Fig. 7 Single fiber reflectance measured (A) and simulated (B) from Intralipid optical phantoms at multiple wavelengths λ = [400, 500, 600, 800, 900] nm; other wavelengths measured in this range not shown to improve clarity.

Download Full Size | PDF

 

Fig. 8 Comparison of (A) single fiber reflectance intensity from MC simulations with experimental measurements from Intralipid phantoms λ = [400 – 900] nm; (B) ${\mu }_s^{\prime }$ estimated from $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ vs. known value.

Download Full Size | PDF

3.2.3. Experimental estimation of ${\mu }_s^{\prime }$ from $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$

The semi-empirical reflectance model given in Eq. (2) was fitted to $R_{\mathit{\text{SF}}}^{\mathit{\text{MC}}}$ simulated in Intralipid for each wavelength-dependent phase function. Table 1 (bottom) lists detailed data for each simulation, including the anisotropy, estimated parameters values, and Pearson correlation coefficients for each fit. These model fits were then used to analyze $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ measured experimentally from Intralipid phantoms and estimate ${\mu }_s^{\prime }$. Prior to this calculation, the $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ data are corrected for slight offset with MC simulations, as noted in Section 3.2.2. Figure 8(B) shows ${\mu }_s^{\prime }$ estimated using $R_{\mathit{\text{SF}}}^{\mathit{\text{meas}}–\mathit{\text{abs}}}$ and the mathematical model vs. the known experimental ${\mu }_s^{\prime }$ value; these data are highly correlated (r = 0.997) and show good agreement, with a mean residual of 9%. As observed in the reflectance data, the relative error in ${\mu }_s^{\prime }$ estimates increases as dimensionless scattering decreases; this is discussed in Section 4.

4. Discussion

This paper presents a methodology to quantitatively extract ${\mu }_s^{\prime }$ from the reflectance intensity collected by a single fiber reflectance device (R_SF). Monte Carlo simulations were used to investigate the dependence of R_SF on the scattering properties of an optically sampled turbid medium. Simulated data were used to identify a relationship between R_SF and ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$, which showed that R_SF approached phase function independent behavior for ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}>\hspace{0.17em}10$. This theory motivated the construction of a solid optical phantom with a very high reduced scattering coefficient ( ${\mu }_s^{\prime }\hspace{0.17em}>\hspace{0.17em}200\hspace{0.17em}{\text{mm}}^{-1}$); such that measurements of R_SF from this phantom are independent of phase function and can be used to calibrate any measured reflectance intensity to an absolute scale. Experimental data show that this methodology can be used to accurately extract ${\mu }_s^{\prime }$ from R_SF measurements of optical phantoms over a wide range of ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.1\hspace{0.17em}-\hspace{0.17em}8\right]\hspace{0.17em}{\text{mm}}^{-1}$ and d_fib = [0.2 – 1] mm for a specified phase function. The work presented here provides an approach to quantitative estimation of both ${\mu }_s^{\prime }$ and μ_a from turbid media, such as tissue.

Monte Carlo simulations were used to investigate the effect of ${\mu }_s^{\prime }$, d_fib, and phase function on R_SF. Inspection of the data showed a relationship between R_SF and ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$, which was characterized by two regimes: (1) a low dimensionless scattering regime where R_SF depends on phase function, and (2) a high dimensionless scattering regime where R_SF becomes independent of phase function and asymptotically approaches a limit. The transition point between these regimes was observed to be approximately ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}\approx \hspace{0.17em}10$. The low dimensionless scattering regime contains reduced scattering values that are experienced during measurements of tissue (common range ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.5\hspace{0.17em}-\hspace{0.17em}5.0\right]\hspace{0.17em}{\text{mm}}^{-1}$) by single fibers with diameters commonly used clinically (common range d_fib = [0.2 – 1.0] mm). While the high dimensionless scattering regime will not be experienced in biological tissue, the insensitivity of R_SF to phase function in this region allows measurements of tissue (of any phase function) to be calibrated onto an absolute scale.

The experimental application of this theory motivated the construction of a solid optical phantom with a large concentration of scattering material (TiO₂) that corresponded to ${\mu }_s^{\prime }$ values that were within the high dimensionless scattering (and phase function insensitive) regime. These solid phantoms allow single fiber reflectance intensities measured from different fibers to be calibrated onto the same absolute scale; previously, this was difficult due to the fiber-diameter dependent effect of and variations in probe-spectralon distance utilized in our previous calibration technique (as noted in Section 2.3). In this study reflectance measurements on the set of phantoms investigated were scaled to a reflectance model to determine the relationship between TiO₂ wt % and ${\mu }_s^{\prime }$. We have assumed that this conversion factor is linear. While it is possible that the relation may not be linear, the data shown on Fig. 6 present a convincing argument that this assumption is appropriate. Uncertainty associated with the ${\mu }_s^{\prime }\hspace{0.17em}/\hspace{0.17em}{\text{TiO}}_2\hspace{0.17em}\text{wt}\hspace{0.17em}\%$ ratio is not expected to introduce error. This is because the highest scattering phantom (30 wt % TiO₂) was selected for scaling subsequent measurements (on liquid phantoms in this study, or on tissue in future work), and in this dimensionless scattering region, a mis-estimation of ${\mu }_s^{\prime }$ for this phantom are not expected be significant due to the relatively small changes in R_SF corresponding to changes in ${\mu }_s^{\prime }$ in this regime.

The semi-empirical model function used to describe the relation between R_SF and dimensionless scattering for each investigated phase function contains 3 fitted parameters, given by Eqs. (2)–(4). The development of this mathematical model provided insight into the factors contributing to the collected signal. The model represents the light collected during measurement as a product of the photons contacting the fiber surface (Φ) and the fraction of those photons within the fiber cone of acceptance (η_c). Figure 4(B) shows that Φ increases with dimensionless scattering, until it saturates as it approaches the 100% limit. This behavior was described using a saturating mathematical expression, which included a phase-function specific fitted power (ρ ₂) on the dimensionless scattering term; while the fitted values for this parameter were relatively close to unity (range=[1.018–1.152], see Table 1) the model fit quality was significantly increased by fitting the parameter. Figure 4(C) shows that η_c is large for small dimensionless scattering values and decreases as it approaches a limiting value. This transition can be rationalized by considering the effect of the optical properties of the turbid medium on the angular distribution of remitted light; this concept has been described in detail previously by Bargo et al. [33]. In brief, small dimensionless scattering values indicate that photons travel relatively large distances between scattering events. Therefore photons remitted from the medium at the fiber face have a high probability of being scattered from (relatively) deep locations in the medium, and in turn are traveling at angles near normal to the tissue interface. As the dimensionless scattering increases, photons are more likely to be remitted from shallow distances into the medium with remitted light exiting diffusely at all angles. Theoretically, η_c in this regime approaches a diffuse limit defined by the projection of the solid angle of the acceptance cone onto the medium/fiber plane divided by all possible solid angles for light remitted across that plane; the mathematics underlying this theory has been described in detail by Horn et al. [41], and applied to single fiber reflectance spectroscopy by Bargo et al. [33]. It should be noted that in the highly scattering region ( ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}\hspace{0.17em}>\hspace{0.17em}10$), the angular distribution of remitted light calculated by the MC model is not perfectly diffuse; instead, a subcomponent of the light collected by the single fiber device has undergone only a few scattering events, which causes slight variations (approximately 5%) in η_c for different phase functions. This phenomenon is known to occur for optical devices with overlapping source and detection areas on a medium interface, and has been described in detail previously by Snabre et al. [42]. As noted in Section 2.2, Eq. (4) is dependent on the fiber NA. Simulations investigating variations of NA on R_SF showed that the effect could be characterized by adjusting η_limit exclusively (without refitting parameters [ρ ₁, ρ ₂, ρ ₃]), with ≈ 5% mean error between estimates of R_SF measured by fibers of NA= [0.22] and NA= [0.1, 0.4] in the biologically relevant scattering range ( ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.44\hspace{0.17em}-\hspace{0.17em}4\right]\hspace{0.17em}{\text{mm}}^{-1}$) over a range of phase functions, with increasing deviations associated with decreasing ${\mu }_s^{\prime }{d}_{\mathit{\text{fib}}}$ (data not shown).

This study includes an experimental proof of concept for the extraction of ${\mu }_s^{\prime }$ from experimental measurements of optical phantoms containing Intralipid. As part of that investigation, R_SF was simulated by a MC model that emulated the experimental measurements in optical phantoms. These techniques yielded R_SF values that were highly correlated, as observed in Fig. 8(A) and 8(B). However, as noted in the Results, there was a slight discrepancy between the simulated and the measured R_SF from Intralipid phantoms, with a 4% increase in measured value vs. the corresponding simulated values. One explanation for this slight mismatch is that the probes used for the experimental measurements have the fiber tip polished at an angle of 15 degrees, altering the fiber face from a circle to an elliptical cross-section, and increasing the area by approximately 3%; this mismatch is expected to represent the bulk of the 4% offset. Additionally, it was observed that the relative error in R_SF estimates increased with decreasing dimensionless scattering. This is most likely caused by slight mismatches between the actual phase function of the liquid phantoms and the model approximation utilized in the MC model [36]. Another source of mismatch may be in the specification of the index of refraction of both fiber and sample in the MC simulation, which was assumed wavelength to be independent. Slight mismatches between experiment and simulation would be observed especially in regions of small dimensionless scattering values where the reflectance signal is most sensitive to phase function.

The goal for performing the absolute R_SF calibration was to extract an absolute value for ${\mu }_s^{\prime }$ within the sampled turbid medium from the measured single fiber reflectance intensity. This was achieved in this study by performing measurements on Intralipid phantoms, calibrating the measured reflectance with measurements on the highly scattering solid calibration phantom, and analyzing the data using a phase-function specific mathematical model describing the dependence of reflectance on dimensionless scattering. The work in this study calculated ${\mu }_s^{\prime }$ for 5 fiber diameters, on Intralipid phantoms of 10 different concentrations at 11 different wavelengths (range [400–900] nm), representing a wide range of phase functions and ${\mu }_s^{\prime }\hspace{0.17em}=\hspace{0.17em}\left[0.05\hspace{0.17em}-\hspace{0.17em}8\right]\hspace{0.17em}{\text{mm}}^{-1}$. Because this calculation utilized the model fit to MC simulations to analyze measured reflectance data, the measurements were corrected for the observed 4% offset prior to ${\mu }_s^{\prime }$ estimation, as described in detail earlier in this Section. After this correction we were able to accurately extract the reduced scattering coefficient from Intralipid measurements, with a mean residual of 9% over the entire range of ${\mu }_s^{\prime }$; this ${\mu }_s^{\prime }$ range is wider than experienced in tissue, and therefore error associated with biologically relevant ${\mu }_s^{\prime }$ values is expected to be < 9%. In this study the semi-empirical model was used in conjunction with measured reflectance on phantoms to provide absolute estimates of ${\mu }_s^{\prime }$, estimates that have been demonstrated to be accurate for a given phase function. In tissue, the phase function is generally unknown and assumptions about the values for ρ ₁ – ρ ₃ that feed into the model expression must be made. Future studies will investigate how the assumed parameter values affect the estimates of ${\mu }_s^{\prime }$, and whether the parameter values can be expressed in terms of metrics extracted from a known wavelength-dependent phase function profile. Furthermore, we are currently investigating methods to measure the parameter values experimentally for biological tissues.