2.1.

Monte Carlo Model

This study adapts the MC model of light transport in layered tissue code developed by Wang et al.²¹ implemented in parallel on a GPU using NVIDIA’s compute unified device architecture by Alerstam et al.^22,23 The MC model for modeling light transport is a stochastic method that simulates light transport in a scattering medium with the probabilities of scattering and absorption events determined by the user-specified optical properties of the medium and the geometry of the light source and measurement probe. Photons step sizes were selected from an exponential distribution that depended on the scattering coefficient, and scattering angles were determined by the scattering anisotropy (g) and the phase function. We used the Heyney–Greenstein phase function. Reflection and refraction due to index of refraction mismatch were calculated using the Fresnel equation and Snell’s law.

A two-layer model was used with reduced scattering in the bottom layer set to zero, and the absorption in the bottom layer was set to $1\times {10}^{15}{\mathrm{cm}}^{-1}$ so that photons reaching the bottom layer were terminated. Scattering anisotropy was held constant at 0.85. The sensitivity of photon path length and sampling depth to phase function and anisotropy (g) has been explored by Kanick et al. for single-fiber spectroscopy.¹⁶ They performed simulations with $g=[0.8,0.9,0.95]$ and with both the Henyey–Greenstein phase function and the modified Henyey–Greenstein phase function. The data showed that path lengths and the sampling depth are independent of anisotropy. The phase function was found to have an observable effect on path length, but the mean sampling depth remained relatively unchanged.

The refractive tissue above the medium was set to 1.452 to match the refractive index of an optical fiber, and the refractive index of the medium was set to 1.33 to match the refractive index of water. The top layer absorption coefficient (${\mu }_{\mathrm{a}}$) ranged from 0 to $30{\mathrm{cm}}^{-1}$ in 20 increments, the top layer reduced scattering coefficient (${\mu }_{\mathrm{s}}^{\prime }$) ranged from 0 to $30{\mathrm{cm}}^{-1}$ in 20 increments, and the top layer thickness ranged from 0 to $3000\mu \mathrm{m}$ in 250 increments. This gave a total of 100,000 separate MC simulations with each using ${10}^7$ photons. The geometry for the simulations is shown in Fig. 1(a). Spatially resolved reflectance was calculated by convolving the impulse response using a Gaussian-shaped beam profile with radius $R_1$, and the reflectance signal was calculated by summing the reflectance values centered at the SDS with a collection fiber radius of $R_2$. For a given set of optical properties (${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ and probe geometry parameters (SDS, $R_1$, $R_2$, we plotted the percentage of photons that never reach depth $Z_0$ versus $Z_0$ [Fig. 1(b)]. If we model the curve in Fig. 1(b) as a sigmoid function, the greatest slope will occur at the depth that is reached by 50% of the photons, meaning that the measured reflectance is most sensitive to optical properties at that depth. Because of this, the sampling depth ($Z_{\mathrm{s}}$) of a probe for a given set of optical properties is defined as the depth reached by 50% of the photons.

Fig. 1

(a) A two-layer geometry was used for the Monte Carlo (MC) simulations. The bottom layer had an absorption coefficient of $1\times {10}^{15}{\mathrm{cm}}^{-1}$ and a scattering coefficient of zero so that photons reaching the bottom layer are terminated. Top layer thickness ($Z_0$) ranged from 0 to $3000\mu \mathrm{m}$ in 250 increments, the top layer absorption coefficient (${\mu }_{\mathrm{a}}$) ranged from 0 to $30{\mathrm{cm}}^{-1}$ in 20 increments. Reflectance measurements were recorded out to 1 cm from the source. (b) A plot of showing the percentage of photons that never reach a depth of $Z_0$ versus $Z_0$ with $\mathrm{SDS}=300\mu \mathrm{m}$, $R_1=100\mu \mathrm{m}$, $R_2=100\mu \mathrm{m}$, ${\mu }_{\mathrm{a}}=1.6{\mathrm{cm}}^{-1}$, and ${\mu }_{\mathrm{s}}^{\prime }=16{\mathrm{cm}}^{-1}$. Sampling depth ($Z_{\mathrm{S}}$) is defined as the depth reached by 50% of the photons.

2.2.

Empirical Measurements of Sampling Depth

To validate the computational results, 12 different phantoms were constructed in order to perform an experimental analysis of sampling depth for DRS of varying SDS. The phantoms were composed of 5-mL solutions of water, India ink (Salis International, Golden, Colorado), and scattering microbeads (Polyscienes, Warrington, Pennsylvania), which spanned absorption and scattering values across a range consistent with those normally found in human tissue. Mie theory was used to determine the scattering properties of the $0.99\text{-}\mu \mathrm{m}$ diameter beads. Mix ratios of water and microbeads were determined so that three different scattering spectra from 11 to $25{\mathrm{cm}}^{-1}$ were achieved at the reference wavelength of 630 nm. Each of these mix ratios was prepared with four different concentrations of India ink, so that the absorption coefficient for the samples ranged from 0 to $23{\mathrm{cm}}^{-1}$, resulting in 12 total phantoms with different scattering and absorption properties, as seen in Table 1.

Table 1

Optical properties of phantoms.

Each of these 12 phantoms was placed into a blackened beaker. Reflection measurements were taken while varying the distance between the probe and the bottom of the beaker from 0 to 3 mm in $50\text{-}\mu \mathrm{m}$ increments. Reflectance spectra were collected at wavelengths from 500 to 700 nm and at SDSs of 370, 740, and $1110\mu \mathrm{m}$. Using the known wavelength dependence of scattering and absorption, ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ were calculated at each wavelength, and for each set of ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ a plot of $P$ versus $Z_0$ was created. These plots were then used to calculate $Z_{\mathrm{s}}$ for each set of optical properties.

2.3.

Mathematical Model of Sampling Depth

The sampling depth for a DRS probe is dependent on the optical properties (${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ and probe geometry parameters (SDS, $R_1$, and $R_2$). The sampling depth data from the MC simulations were accurately described by the equation

Equation (1) is an empirical expression that accurately describes the MC sampling depth data. This expression was found by trying thousands of candidate functions with the help of TableCurve 2-D (Automated Curve Fitting and Equation Discovery Software, Systat, 2002).

Equation (1) has four free parameters ($a_1$, $a_2$, $a_3$, $a_4$) whose values must be determined by fitting the MC data. This was accomplished by minimizing the residual between the MC sampling depth results and the sampling depths calculated using Eq. (1) and a Levenberg–Marquardt algorithm scripted in MATLAB®. The dependence of the free parameters on SDS was then determined so that Eq. (1) could be used to determine sampling depth for a given probe geometry and set of optical properties. SDS and $Z_{\mathrm{S}}$ are in units of cm and ${\mu }_{\mathrm{a}}$ and ${\mu }_{\mathrm{s}}^{\prime }$ are in units of ${\mathrm{cm}}^{-1}$.

3.1.

Experimental Validation

The sampling depth results from the phantom experiments were used to validate the computational sampling depth results at SDSs of 370, 740, and $1100\mu \mathrm{m}$ with optical properties in the range ${\mu }_{\mathrm{a}}\in [0-25]{\mathrm{cm}}^{-1}$, and ${\mu }_{\mathrm{s}}^{\prime }\in [0-30]{\mathrm{cm}}^{-1}$. Figure 2 shows an overlay of the computational (transparent mesh) and the experimental (colored surface) results and provides a visual illustration of the good agreement between experimental and computational results. The root-mean-squared percent error for an SDS of $370\mu \mathrm{m}$ was 1.71%, for an SDS of $710\mu \mathrm{m}$ it was 1.27%, and for $1100\mu \mathrm{m}$ it was 1.24%. This agreement indicates that the MC model accurately returns sampling depth. The ripples in the data indicate that the agreement between the phantom measurements and the MC data is wavelength dependent. We believe this is due to the use of the inverse power law for the wavelength dependent description of scattering in the phantoms containing polystyrene microbeads, when in reality, the true scattering values of the phantoms as a function of wavelength contain “humps” in the curve due to the relatively narrow size distribution of the microspheres.

Fig. 2

Plots of $Z_{\mathrm{S}}$ predicted by Monte Carlo modeling versus the experimental values for $Z_{\mathrm{S}}$ at source-detector separations (SDS) of (a) 370, (b) 740, and (c) $1100\mu \mathrm{m}$. An overlay of two-dimensional (2-D) surfaces showing the relationship between scattering and absorption on sampling depth for both Monte Carlo and experimental results. These plots provide a visual illustration of the agreement between the computational (transparent mesh) and experimental (colored surface) results for source-detector separations of (d) 370, (e) 740 and (f) $1100\mu \mathrm{m}$.

3.2.

Analytical Model of Sampling Depth

The analytical model of sampling depth shown in Eq. (1) was fit to MC data for a probe with fiber diameters of 50, 100, 200, and $400\mu \mathrm{m}$. For each fiber diameter, SDS ranges from adjacent fibers to $1000\mu \mathrm{m}$, ${\mu }_{\mathrm{s}}^{\prime }$ ranges from 3 to $40{\mathrm{cm}}^{-1}$, and ${\mu }_{\mathrm{a}}$ ranges from 0 to $40{\mathrm{cm}}^{-1}$. Fitting parameters $a_1$ and $a_2$ were found to have a linear relationship with SDS, $a_3$ was found to have a quadratic relationship with SDS, and $a_4$ is a constant. Table 2 gives the fitting parameters for the four different common fiber diameters as a function of SDS. Table 2 allows the fitting parameters in Eq. (1) to be determined for a given fiber diameter and SDS so that Eq. (1) can be used to determine sampling depth for a specific probe geometry.

Table 2

Values for fitting parameters at various fiber diameters.

Figure 3 below shows the sampling depth predicted by the analytical model in Eq. (1) and Table 1 versus the MC sampling 3 depth. Model predictions were strongly correlated with the MC data with a mean residual error of 2.89%.

Fig. 3

Monte Carlo results for simulation of sampling depth versus sampling depth predictions from the analytical model for all four fiber diameters [Eq. (1)]. The line of unity is shown for comparative purposes. There is a 2.89% error between the Monte Carlo simulation results and the analytical model results.

3.3.

Effect of Anisotropy and Phase Function on Sampling Depth

Because scattering anisotropy and the choice of a phase function can impact reflectance at short SDSs,²⁴ a subset of MC simulations was performed to investigate the effect of anisotropy and phase function on sampling depth. The data showed no change in sampling depth for simulations of different anisotropy values ($g=[0.80,0.85,0.90,0.95]$) over a range of optical properties (${\mu }_{\mathrm{a}}\in [0–25]{\mathrm{cm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }\in [0–30]{\mathrm{cm}}^{-1}$) and probe geometries ($\mathrm{SDS}\in [50–800]\mu \mathrm{m}$), $R\in [50–400]\mu \mathrm{m}$. This is illustrated in Fig. 4(a), where sampling depth versus SDS is plotted for a probe with $50\text{-}\mu \mathrm{m}$ diameter fibers with a reduced scattering coefficient of $10{\mathrm{cm}}^{-1}$ and an absorption coefficient of $10{\mathrm{cm}}^{-1}$ for four different anisotropy values. The mean percent error across all anisotropy values for all probe geometries and optical properties was 3.47%. Additionally, the data showed no change in sampling depth for simulations performed with the Heyney–Greenstein (HG) phase function or the modified Heyney–Greenstein (MHG) phase function. This is illustrated in Fig. 4(b), where sampling depth versus SDS is plotted for a probe with $50\text{-}\mu \mathrm{m}$ diameter fibers with anisotropy values of 0.85, a reduced scattering coefficient of $10{\mathrm{cm}}^{-1}$ and an absorption coefficient of $10{\mathrm{cm}}^{-1}$ for both the HG and MHG phase functions. A change in g or the phase function did affect the reflectance values; however, there was no change in the sampling depth as it is defined in this study. These results agree with the findings by Kanick et al. for single-fiber reflectance spectroscopy¹⁶ that show sampling depth is unaffected by both the anisotropy value and the choice of phase function.

Fig. 4

(a) Sampling depth versus SDS for varying anisotropy values with ${\mu }_{\mathrm{a}}=10{\mathrm{cm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=10{\mathrm{cm}}^{-1}$, and fiber diameter at $50\mu \mathrm{m}$. (b) Sampling depth versus SDS for both HG and MHG phase functions with ${\mu }_{\mathrm{a}}=10{\mathrm{cm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=10{\mathrm{cm}}^{-1}$, $g=0.85$ and the fiber diameter at $50\mu \mathrm{m}$.