2.1.

Modeling Reflectance

In the context of DRS, a forward model is any function that quantitatively describes the relationships between measured reflectance and the multiple input parameters that will influence reflectance, such as the tissue optical properties and the geometry specifications of the optical measurement. The primary motivation for understanding these reflectance relationships is not for use in forward direction calculations, but rather for use in solving the inverse problem (calculating tissue optical properties from measured reflectance). The relationship between forward and inverse models is expressed generically in Eq. (1), where $R$ represents measured reflectance, $p(\theta )$ is the scattering phase function, and $G$ represents all probe geometry parameters.

2.2.

Phase Functions

The phase function, $p(\theta )$, of a scattering sample specifies the probability distribution of scattering angles and is typically normalized to 1, integrated over all solid angles. To generalize the phase function, the anisotropy factor, $g$, is defined as the average cosine of the scattering angle, $\langle \cos (\theta )\rangle$. In tissue, $g$ ranges, typically, between 0.75 and 0.99.³¹ The scattering coefficient, ${\mu }_s$, and the anisotropy factor, $g$, are often combined into a single parameter, known as the reduced scattering coefficient: ${\mu }_s^{\prime }={\mu }_s(1-g)$. This similarity relationship is based on the assumed condition that observed optical measurements are equivalent for any combination of $g$ and ${\mu }_s$ that results in the same ${\mu }_s^{\prime }$.³²

Interest in the influence of the details of the phase function on reflectance measurements was originally reported by Mourant et al., Bevilacqua and Depeursinge, and Kienle et al., who demonstrated that reflectance is dependent on the specific form of phase function for reflectance measured close to the source ($<1\mathrm{mm}$ source-detector separation in tissue), thus negating the applicability of the reduced scattering coefficient similarity condition for short distances.^19,33,34 Bevilacqua and Depeursinge further explored the relationship between phase function and reflectance by considering higher moments of the phase function and by developing a new similarity relationship applicable for small source-detector separations.^3,19 More recently, Kanick et al. have also utilized this similarity function for single-fiber source-detector geometries.³⁵ In the derivation of the new similarity relationship, the phase function is expanded into a series of Legendre polynomials, where $g_n$ are the $n$’th-order Legendre moments. The first-order moment, $g_1$, is equivalent to the familiar anisotropy factor, $g$. To calculate its value from any arbitrary phase function, the following equation can be used:³⁶

The use of $\gamma$ is applicable to any phase function representation, including the commonly adopted Henyey-Greenstein (HG), modified HG (MHG), and more rigorous Mie theory formalisms, as described below.

2.2.2.

Mie theory phase function

A more rigorous estimation of the phase function in tissue can be calculated from Mie theory by approximating tissue as a distribution of discretely sized spherical particles. By combining the individual phase functions from each sphere size, weighted by relative concentrations of different sizes, an average phase function can be constructed to represent the bulk medium. One proposed approach to selecting the appropriate sizes and distribution of particles to best approximate tissue is to use a fractal distribution, based on observations that refractive index variations measured by phase-contrast microscopy scale according to a power law.^40–42 The number density, $N_i$, of each particle size $d_i$ is then defined as

Eq. (11)

$$N_i={(d_0/d_i)}^{\alpha },$$

where $d_0$ is the representative size of a spherical volume containing the particles, and $\alpha$ is the fractal volume dimension that relates to the ratio of small versus large particles.⁴⁰ A large $\alpha$ value corresponds to a greater ratio of small to large particles.

For each discrete particle size, Mie theory is used to calculate the scattering cross-section, ${\sigma }_s(d_i)$, and the anisotropy factor, $g(d_i)$, using values for the medium and particle indices of refraction that are representative of cytoplasm ($n_m$) and cellular organelles ($n_p$), respectively. Schmitt and Kumar calculated average values of $n_m=1.352$ and $n_p=1.39$ to 1.45.⁴² Using the assumption that waves scattered by individual particles add linearly, the anisotropy factor is defined as

Eq. (12)

$$g=\frac{\sum _{i=1}^m{d}_i^{-\alpha }{\sigma }_s(d_i)g(d_i)}{\sum _{i=1}^m{d}_i^{-\alpha }{\sigma }_s(d_i)},$$

where $m$ is the total number of discrete particle sizes in the distribution.^43,44 The combined phase function, $p(\theta )$, can be calculated from a linear combination of the scattering amplitude, $S_{11}$ (as defined by Bohren and Huffman in their formulation of Mie theory⁴⁵), scaled by the number density of spheres.

Eq. (13)

$$p(\theta )=\frac{\sum _{i=1}^m{d}_i^{-\alpha }{S}_{11}(\theta ,d_i)}{\sum _{i=1}^m{d}_i^{-\alpha }}.$$

Calculation of $g_1$ and $g_2$ can then be performed using Eqs. (5) and (6).

The use of Mie theory and the fractal distribution model of scattering sizes offers a unique opportunity to examine the relationship between $g$ and $\gamma$. Throughout this paper, calculations are performed using refractive index values and particle diameter values from the literature:^42,43 $n_p=1.42$, $n_m=1.352$, $d=[10,100,150,400,600,800\mathrm{nm},1,3,10,20\mu \text{m}]$, and a wavelength of $\lambda =600\mathrm{nm}$. The fractal dimension, $\alpha$, is varied between 2 and 6, each value of $\alpha$ corresponding to a unique pair of $g$ and $\gamma$. Different selections for $n_p$, $n_m$, $\lambda$, and the discrete particle sizes will yield slightly different numerical results for the relationship between $\gamma$ and $g$, but, in general, the trends remain the same.

Figure 2 demonstrates that the relationship between $g$ and $\gamma$ is nonlinear and that $\gamma$ increases more rapidly for large $g$ values, especially when $g>0.8$. This is especially important since most tissues have $g$ values in this higher range. Further, the values of $\gamma$ in this region are generally $>2.0$ (dependent on the values of $n_p$, $n_m$, $d$, and $\lambda$), which is the upper limit of $\gamma$ for the HG and MHG phase functions, indicating that a phase function constructed from Mie theory is more appropriate when modeling tissue. Values of $\gamma$ in tissue have been found to range from 1.7 to 2.2, which highlights the importance of using a phase function capable of $\gamma$ values $>2$.³ Bevilacqua and Depeursinge have previously analyzed the space of possible combinations of values for $g_1$, $g_2$, and $\gamma$ for various phase functions, noting the larger range of possible combinations for Mie phase functions as compared to HG and MHG.¹⁹ The analysis here differs by limiting the Mie phase functions to those constructed using the fractal particle distribution model, providing a simplified understanding of the relationship between $g$ and $\gamma$ that is specifically representative of tissue.

Fig. 2

Relationship between the anisotropy values, $g$ and $\gamma$ (for the specific parameters $n_p=1.42$, $n_m=1.352$, $\lambda =600\mathrm{nm}$), as calculated from Mie theory using a fractal distribution of scattering particle sizes ($d=[10,100,150,400,600,800\mathrm{nm},1,3,10,20\mu \text{m}]$) with the fractal dimension, $\alpha$, varied between 2 and 6. Differences in refractive index values and/or particle sizes will provide different numerical results, but the generalized shape of the curve will remain consistent.

For the purposes of comparison, Fig. 3 illustrates constructed functions for the HG, MHG, and Mie phase function descriptions. All three functions are for $g=0.7$. However, the HG function has $\gamma =1.7$, while the MHG and Mie functions have values of $\gamma =1.4$, which can be seen in the differences in high-angle scattering between the HG function and the MHG and Mie functions. Note also that, despite having the same $\gamma$ value, the MHG and Mie phase functions are not identical. Differences between the two result from, and can be observed in variations in the higher-order Legendre moments ($g_n$).

Fig. 3

Comparison of a Henyey-Greenstein (HG), a modified Henyey-Greenstein (MHG), and a Mie phase function, all with $g=0.7$, but with different $\gamma$ values (HG: $\gamma =1.7$, MHG and Mie: $\gamma =1.4$).

It is important to note that these phase functions are presented using the standard representation that specifies the distribution of light intensity scattered as a function of angle [$p(\theta )$], i.e., the relative power per unit solid angle that is observed at a given angle in a steady-state system. Alternatively, the phase function can be expressed as a probability distribution of scattering angles [$p^{\prime }(\theta )$], i.e., the probability that a photon will scatter at a given angle. The difference reflects the inherent field versus photon representation issue common in optics; the two forms are related as $p^{\prime }({\theta }_j)=p({\theta }_j)\sin {\theta }_j$. In the literature, phase functions have been reported in both ways, without specifying which form is being used, which often leads to confusion.^46,47 The angular representation, $p^{\prime }(\theta )$, needs to be used when sampling scattering direction changes in Monte Carlo simulations.

3.1.

Monte Carlo Simulations

“Experimental” data for the work in this paper were obtained from Monte Carlo simulations. Computational simulation allows for a controlled and exact approach to the comparison of various models and input parameters, which would be more difficult experimentally using phantoms. A version of Monte Carlo for multilayered tissue⁴⁸ adapted for use on a GPU was used for increased speed.⁴⁹ This code was modified to launch and collect photons within the area and numerical aperture of user-defined fibers to model the two-fiber probe geometry (one illumination and one collection fiber), as illustrated in Fig. 1. Additionally, the code was modified to provide a choice in phase function, allowing the user to select $p_{\mathrm{hg}}$, $p_{\mathrm{mhg}}$, $p_{\mathrm{mie}}$, or any other user-defined phase function at run time.

Monte Carlo simulations were run in sets, each having identical geometry and phase function input parameters. Table 1 reports the values of all input parameters used for the simulations. In each set, simulations were run for a set of combinations of reduced scattering coefficient values, ${\mu }_s^{\prime }$, and absorption coefficient values, ${\mu }_a$. For all simulations, constant values for the refractive index of the tissue, the refractive index of the fibers, the numerical aperture of the fibers, and the fiber diameter were used. All simulations were terminated upon the collection of 10,000 photons, such that the variances of the simulations were equivalent at $\sim 1\%$ according to Poisson statistics. Because of the equivalence of variances, and their small 1% value, error bars are not included on plots of results.

Table 1

Input parameters for Monte Carlo simulations. Parameters listed with an asterisk (*) are constant for all simulations.

To examine the effects of source-detector separation, sets of simulations were run for a range of values, $\rho$. And for exploration of the effects of phase function, simulations were run assuming an MHG phase function form using various combinations of $g$ and $\gamma$. All combinations of the two variables were constructed, but due to the mathematical limitations of the form of the formula for MHG, phase function combinations with $\gamma =1.9$ and $g=0.75$ and 0.85 are not possible and, thus, could not be simulated and evaluated. Mie theory–based phase functions were also considered, using constant refractive index values, particle sizes $d$, and wavelength $\lambda$, while varying the fractal dimension, $\alpha$. Resulting anisotropy values of $g$ and $\gamma$ were calculated using Eqs. (5), (6), (7), and (13). The values were selected such that the values of $\gamma$ were spaced equally between 1.2 and 2.3; the range of $\gamma$ was chosen to correspond to physiologically relevant anisotropy values between $g=0.55$ and 0.96 (see Fig. 2). The particle sizes and refractive index values are representative of tissue.^42,43

While only two phase function models (MHG and Mie) were chosen here for analysis, note that regardless of the mathematical form of the phase function, the $\gamma$ similarity relationship dictates that reflectance values will be equal for equal values of $\gamma$, and thus, the details of the higher-order moments of various models is not significant when considering reflectance model trends. Nonetheless, we note that only Mie theory can be invoked for values of $\gamma >2.0$.

3.2.

Inverse Modeling: A Different Approach to Error Analysis

As previously demonstrated in the literature,^19,26,35 reflectance models depend on phase function, and hence, the choice of a model to use in the inverse direction requires knowledge of the tissue phase function. However, in most cases, this information is not known a priori and, consequently, assumptions must be made. An incorrect assumption will affect the accuracy of extracted scattering and absorption properties from measured reflectance; we wish to quantify the degree to which an incorrect assumption affects these extracted values. Unfortunately, because there is currently no gold standard for measuring the optical and physiological parameters from point measurements in tissue, it is not possible to experimentally test the effects of model differences. Therefore, the error analysis presented in this work uses computationally constructed reflectance data, which allows for a more controlled and exact approach to the comparison of various models and input parameters.

This analysis first involves computationally constructing a set of reflectance spectra that are representative of the true model and tissue properties. Then the assumed reflectance model, which is different from the true model due to inaccurate assumptions about tissue phase function, is used in the inverse direction to extract the fitted optical property values. The differences between these extracted values and the actual property values, thus, provide a measure of error that represents the potential uncertainty in in vivo measurements. We note moreover, that this technique is novel in its approach to examining error in terms of extracted physiological values from reflectance spectra, as opposed to reporting error exclusively of the optical coefficient values.

It is difficult to define a single estimate of error because the errors in extracted quantities are directly related to ${\mu }_s^{\prime }$ and ${\mu }_a$ themselves. Hence, errors are calculated for multiple reflectance spectra, each with different scattering and absorption properties. A set of wavelength-dependent ${\mu }_a$ and ${\mu }_s^{\prime }$ spectra are first calculated using Eqs. (3) and (4), over a wavelength range from 400 to $800\mathrm{nm}$ at a resolution of 1 nm (for a total of 400 data points). Values for the five physiological fitting variables ($a$, $b$, $f_1$, $f_2$, and $r$) were selected such that the resulting ${\mu }_a$ and ${\mu }_s^{\prime }$ spectra would span the ranges characteristic of tissue. For scattering, four spectra were selected to represent overall low scattering ($a=8$, $b=1$), medium scattering ($a=16$, $b=1$), high scattering ($a=24$, $b=1$), plus one with a large slope representing a high density of smaller scatterers ($a=16$, $b=2$). The reference wavelength ${\lambda }_0$ was selected to be 600 nm. For absorption, spectra were selected for a physically representative range of absorption values (blood volume fraction, $f_1=0.5$, 2, and 5%), all with a moderate hemoglobin oxygen saturation level ($f_2=75\%$). For all cases, the average blood vessel radius is held constant at $r=10\mu \text{m}$ which is representative of capillaries. Using all combinations of the scattering and absorption conditions, a set of 12 unique scattering and absorption spectra were used for testing.

Computationally constructed reflectance spectra, representing measured data, are calculated using the test scattering (${\mu }_s^{\prime }$) and absorption (${\mu }_a$) spectra according to Eq. (1), where the forward reflectance model is specific to a single phase function [$p(\theta )$] and measurement geometry ($G$). Reflectance values are calculated using the forward model via a lookup table composed of Monte Carlo simulated reflectance data specific to the unique pair of $p(\theta )$ and $G$. A bilinear interpolation algorithm (available in MATLAB®) was used to calculate reflectance values for combinations of ${\mu }_s^{\prime }$ and ${\mu }_a$ not explicitly listed in the table. Noise of 1% was added to the spectra to better represent experimentally acquired data.

As previously mentioned, when using reflectance models in the inverse direction, the least-squares fitting algorithm is insensitive to the specific form of the forward model as long as the modeled reflectance is negligibly different over the set of ${\mu }_a$ and ${\mu }_s^{\prime }$ values. Hence, the choice of model format to use in the inverse direction can be based primarily on preference. For flexibility, we have chosen to use an interpolation lookup table method (with the same lookup table format used to construct the reflectance spectra). This is easily implemented in MATLAB® using the lsqcurvefit toolbox (using the Levenberg Marquardt optimization algorithm), as it can be constructed with any user-defined computational function for the forward model regardless of functional form.

When calculating optical property extraction error associated with inaccuracies in phase function assumptions, reflectance spectra were calculated for all combinations of ${\mu }_s^{\prime }$ and ${\mu }_a$ spectra listed above, yielding a set of 12 reflectance spectra. The assumed reflectance model, representing an inaccurate choice of phase function, is then used in the inverse direction to fit each of the 12 true reflectance spectra and extract the fitting coefficients ($a$, $b$, $f_1$, $f_2$, and $r$) and, consequently, the ${\mu }_s^{\prime }$ and ${\mu }_a$ spectra. Errors in extracted optical and physiological values are reported as an average across all 12 combinations on a percentage basis for the physiological parameters ($a$, $b$, $f_1$, $f_2$, and $r$), and as the mean percentage error across the entire wavelength range for the optical coefficients (${\mu }_a$ and ${\mu }_s^{\prime }$).

4.1.

Phase Function Dependence

Many reflectance models, specifically empirical models describing reflectance at small source-detector separations, have assumed an HG phase function with an anisotropy value $g=0.9$.^17,18,20,22 As introduced in Sec. 2.2, the MHG and the Mie phase functions have been shown in the literature to be more representative of true scattering in tissue as a consequence of the higher probability for large-angle scattering.

To examine the impact of these phase functions that are more representative of tissue, specifically the influence of both $g$ and $\gamma$, reflectance relationships were simulated for a range of $g$ and $\gamma$ values. The MHG phase function is user-friendly for this purpose, since it conveniently enables independent variation of $g$ and $\gamma$. Phase functions constructed using values of $g$ and $\gamma$ described in Sec. 3.1 were examined using a source-detector separation of 250 μm. Figure 4 illustrates the resulting reflectance relationships when each of the phase functions is used for simulation.

Fig. 4

Reflectance relationship plots for a probe geometry with $\rho =250\mu \text{m}$, and different MHG phase functions. Colors and symbols correspond to different $\gamma$ values and line styles correspond to different g values. Note that plots for the combinations $g=0.95$ and $\gamma =1.7$ and 1.9 are not present due to the mathematical limitations of the MGH phase function with these values.

In much of what follows, reflectance models are reported graphically in a set of two plots, as illustrated in Fig. 4. The left plot illustrates the relationship between the relative reflectance and the reduced scattering coefficient, ${\mu }_s^{\prime }$, in the absence of absorption, and the right plot illustrates the relationship between reflectance and the absorption coefficient, ${\mu }_a$. In the plot of reflectance versus ${\mu }_a$, reflectance values are normalized to 1 at ${\mu }_a=0$. This highlights the shape differences rather than the overall amplitude differences in reflectance, as caused by differences in ${\mu }_s^{\prime }$. The two subplots correspond to the two components of the empirical model format inspired by Beer’s Law, presented in Eq. (2); the first represents Eq. (2) when ${\mu }_a=0$ such that $R^{\mathrm{REL}}=R_0({\mu }_s^{\prime })$, and the second represents Eq. (2) normalized such that $R^{\mathrm{NORM}}=R/R_0=e^{-{\mu }_a\langle L\rangle }$. This graphical representation is useful for identifying trends, as it helps to distinguish the dependences of scattering and absorption separately.

Most importantly, Fig. 4 reveals that when $\gamma$ is held constant, the influence of $g$ is remarkably small. Conversely, when $g$ is held constant, reflectance values vary significantly with $\gamma$. At ${\mu }_s^{\prime }=5{\mathrm{cm}}^{-1}$ and ${\mu }_a=0$, reflectance for $\gamma =1.3$ is more than three times larger than for $\gamma =1.9$. Differences due to $\gamma$ in the reflectance versus ${\mu }_a$ relationships are more moderate, but still show a significant dependence. This observation for a two-fiber geometry reflectance relationship follows from previous work by both Bevilacqua et al. and Kanick et al. (for radially resolved multifiber geometries and for single-fiber geometries) in that the critical phase function variable to consider when building reflectance models (for short source-detector distances) is $\gamma$, and not $g$.^3,19,35 We further emphasize not only that $\gamma$ is more influential than $g$, but also that the influence of $g$ is limited; phase functions with different $g$ values, but the same $\gamma$ value, have reflectance relationships that are essentially similar.

Because the MHG phase function is restricted to $\gamma$ values $<2$, and yet Mie phase functions constructed to represent tissue have possible $\gamma$ values $>2$, we also examined how Mie phase functions, which are more representative of tissue, influence reflectance relationships for an even wider range of $\gamma$. These results are presented in Fig. 5. As compared to the reflectance relationships built with the MHG phase function, the differences due to $\gamma$ in the reflectance versus ${\mu }_a$ relationships are qualitatively similar. The differences in reflectance for the reflectance versus ${\mu }_s^{\prime }$ relationships, however, have increased such that, at ${\mu }_s^{\prime }=5{\mathrm{cm}}^{-1}$, reflectance for $\gamma =1.2$ is nearly six times larger than for $\gamma =2.3$. Further, as $\gamma$ increases, the differences in reflectance become more substantial. For example, when considering $g$ values most representative of tissue ($g>0.85$), the differences in reflectance values of $\gamma$ between 1.75 and 2.3 are comparatively greater than those for values between $\gamma =1.2$ and 1.75, even though, in the first case, $g$ only ranges from 0.86 to 0.96, whereas in the second, $g$ ranges from 0.55 to 0.86. This indicates that the impact of $\gamma$ is even stronger for $g$ values most commonly expected in tissue.

Fig. 5

Reflectance relationship plots for a probe geometry with $\rho =250\mu \text{m}$, for a range of Mie phase functions, each with unique g and $\gamma$ values.

We offer an intuitive explanation for trends observed regarding the influence of $\gamma$ on reflectance, which can be explained by the relative number of large-angle scattering events that are expected to occur. For the reflectance versus ${\mu }_s^{\prime }$ relationships, reflectance is smaller when $\gamma$ is larger. At small source-detector-separations, such as that used in the above simulations (250 μm), for a photon to be collected, it can only scatter a modest number of times, and at least one of those events must be at a large angle such that the photon can turn around and reach the collector. Lower $\gamma$ phase functions have larger probabilities for these backward-scattering events (see Fig. 3) and, therefore, allow for more photons to be redirected back toward the fibers and collected. Conversely, the low probability of backscattering characterized by a high-$\gamma$ phase function will result in photons traveling far from the source and detector due to the predominance of forward scattering, and when the photon does experience a large-angle scattering event, it will likely be so far away from the detector that the probability for collection is low. This effect is amplified for small values of ${\mu }_s^{\prime }$, since average path lengths are longer, and therefore, photons will propagate much farther from the collection fiber with a fewer number of scattering events. With fewer opportunities for large-angle scattering, fewer photons will reach the collector.

The influence of backscattering can also explain why reflectance is greater for high $\gamma$ values in the reflectance versus ${\mu }_a$ relationships. As we have established, with a larger $\gamma$, the probability of high-angle scattering is lower, and the probability of forward scattering is larger. Many forward-scattering events will move the photon away from the detector very quickly. Therefore, the most likely scenario for a photon to be collected would then be to experience one high-angle scattering event in conjunction with a small number of forward-scattering events. When $\gamma$ is lower, the photon experiences more scattering events at moderate scattering angles, which allows it to scatter more times while still staying near the collector, thus enhancing its probability for being collected. This larger number of scattering events corresponds to a longer average path length, and, as determined by Beer’s law, a longer traveled path length results in higher absorption and, thus, a lower reflectance.

4.1.1.

Error analysis

Under practical conditions, when reflectance from a single collection fiber is being used for an inverse analysis to extract optical property values, assumptions about the value of $\gamma$ and its corresponding reflectance model must be made since the exact value of $\gamma$ is typically unknown. With most tissues of interest, even if the anisotropy value, $g$, is assumed to be $>0.85$, the range of $\gamma$ is still quite large ($\gamma =1.75$ to 2.3), as is demonstrated in Fig. 2. Within this range, it is then logical to assume an average value of $\gamma$ near 2.0. To estimate the effect that an inaccurate assumption will have on the accuracy of the extracted optical coefficient values, the novel error analysis method, described in Sec. 3.2, was employed. A set of 12 reflectance spectra were constructed from the ${\mu }_s^{\prime }$ and ${\mu }_a$ test spectra, using the reflectance model for $\gamma =1.75$. The reflectance model for $\gamma =2.05$ was then used in the inverse direction to extract the optical property values from the reflectance spectra. This represents a situation where the tissue being measured has a phase function with $\gamma =1.75$, but the inverse model incorrectly assumes a phase function with $\gamma =2.05$. This analysis was also performed assuming a large $\gamma$ value, $\gamma =2.3$. Figures 6 and 7 illustrate the errors in optical coefficient values that result from this analysis for an example case with tissue physiological properties $a=16$, $b=2$, $f_1=5\%$, $f_2=75\%$, and $r=10\mu \text{m}$.

Fig. 6

(a) Fitting results of forward model to reflectance data built using a phase function with $\gamma =1.75$, but fit using a model for a phase function with $\gamma =2.05$. (Optical property spectra built assuming $a=16$, $b=2$, $f_1=5\%$, $f_2=75\%$, and $r=10\mu \text{m}$.) (b) The real and extracted optical coefficient spectra.

Fig. 7

(a) Forward model fitting results to reflectance data built using a phase function with $\gamma =2.3$, but fit using a model for a phase function with $\gamma =2.05$. (Optical coefficient spectra built assuming $a=16$, $b=2$, $f_1=5\%$, $f_2=75\%$, and $r=10\mu \text{m}$. (b) The real and extracted optical coefficient spectra.

When $\gamma$ is underestimated, the ${\mu }_s^{\prime }$ and ${\mu }_a$ extracted values are higher than expected, and when γ is overestimated, the ${\mu }_s^{\prime }$ and ${\mu }_a$ extracted spectra have values lower than expected. The mean errors in optical and physiological values, averaged over both cases, are presented in Table 2. Note that these errors represent comparison of $\gamma$ values with differences $\le 0.3$. In the extreme case, when comparing forward models with differences in $\gamma >1.0$, extraction errors in ${\mu }_a$ and ${\mu }_s^{\prime }$ were found to exceed 100%. But in realistic situations, assuming an average $\gamma$ for the inverse fit is much more reasonable.

Table 2

Average errors in the extracted physiological and optical property values by incorrectly assuming a model with γ=2.05, when the true tissue has γ=1.75 or γ=2.3.

True $\gamma$	$a$	$b$	$f_1$	$f_2$	$r$	${\mu }_{s,\mathrm{mean}}^{\prime }$	${\mu }_{a,\mathrm{mean}}$
$\gamma =1.75$	15.6%	$-\mathbf{20.1}\%$	4.9%	$-1.6\%$	1.2%	15.3%	6.0%
$\gamma =2.3$	$-\mathbf{14.9}\%$	27.6%	$-4.5\%$	0.2%	$-3.1\%$	$-\mathbf{14.5}\%$	$-5.6\%$

Note: Largest errors are indicated in bold.

We remind the reader that the errors reported here are the averages over the 12 combinations of ${\mu }_s^{\prime }$ and ${\mu }_a$ test spectra. However, the errors are also dependent on the specific ${\mu }_s^{\prime }$ and ${\mu }_a$ values being used. For scattering in particular, the low-scattering spectra ($a=8{\mathrm{cm}}^{-1}$, $b=1$) will result in extraction errors in ${\mu }_s^{\prime }$ that are two to three times larger than for the high-scattering spectra ($a=24{\mathrm{cm}}^{-1}$, $b=1$). Although partially due to the fact that percentage errors are greater for smaller values, this is primarily due to the greater differences in reflectance at smaller ${\mu }_s^{\prime }$ values from variations in $\gamma$ (see Figs. 4 and 5). Because differences in reflectance as a function of ${\mu }_a$, which are due to variations in $\gamma$, are much smaller than differences in reflectance as a function of ${\mu }_s^{\prime }$, the errors in extracted ${\mu }_a$, $f_1$, $f_2$, and $r$ are also significantly smaller.

Note on goodness-of-fit

A new and remarkable observation from Figs. 6 and 7 is that, despite significant errors in the extracted optical properties, the reflectance spectra for the forward model fit are almost perfect! As compared to the fit of the reflectance spectra when the true and assumed values of $\gamma$ match, the fits of the reflectance spectra when the true and assumed values of $\gamma$ are different result in fittings of the reflectance model that are only slightly inferior. This is quantitatively verified by comparing the mean normalized residual of the fits: $<R_{\mathrm{fit}}-R_{\mathrm{true}}|/R_{\mathrm{true}}>$ over all wavelengths. For matching $\gamma$, the average of the mean normalized residuals (over the 12 test spectra) is 0.0045, and for nonmatching $\gamma$, the average value is 0.0056, both of which are within the 1% noise that was added to the reflectance spectra.

The implication of this novel observation is that accuracy in the extracted optical coefficient values cannot be associated with the goodness-of-fit of the forward reflectance model to the data. This is of particular significance because the validity and suitability of a given reflectance model is often reported as being verified by the goodness-of-fit that it demonstrates when used to fit experimental reflectance data. In reality, the goodness-of-fit measure can only be used to verify a model if that model accounts for all sources of potential variability (phase function, source-detector separation, fiber diameter, fiber and probe material indices of refraction, etc.), and if the wavelength-dependent assumptions concerning the spectral shapes of ${\mu }_s^{\prime }$ and ${\mu }_a$ (when using wavelength-dependent fitting as is employed in this work) are well representative of the tissue being measured. In practice, it is prudent not to assume a significant level of certainty in these factors to make conclusive verifications based on goodness-of-fit alone.

While this observation regarding goodness-of-fit as a deceptive measure of model accuracy is most critically applicable to uncertainties in phase function, it is also applicable to conditions for which other model variables are inaccurate. For example, when the reflectance model being used to fit a measured spectrum is built using incorrect values for probe geometry details (e.g., numerical aperture, fiber diameter, etc.), similar inaccuracies in extracted physiological and optical properties will result despite excellent spectral fits.⁵⁰ This suggests that building reflectance models that are unique to the specifications of individual measurement probes is important for maximal accuracy; the use of GPU-enabled Monte Carlo codes makes this much more feasible, especially when the data can be used directly in lookup tables for inverse model fitting.

4.2.

Effects of Source-Detector Separation

To explore the influence of fiber separation, i.e., the distance between the source and detector fibers, reflectance relationships were constructed from Monte Carlo simulations for a number of source-detector separations, $\rho =[250,500,1000,3000,2000,5000,\mathrm{10,000}]\mu \text{m}$. For each separation distance, several phase functions constructed from Mie theory were used. The separation distances were chosen to represent the nondiffuse regime ($\le 500\mu \text{m}$), the diffuse regime ($\ge 5\mathrm{mm}$), and the intermediate region (1, 2, and 3 mm). The intermediate region is of particular interest because reflectance models have not been extensively investigated for these source-detector separations. Thus, by examining all three regimes simultaneously, a more comprehensive understanding of how tissue optical properties influence reflectance can be attained.

As has been previously utilized by others,^25,26 we employ a convenient method of investigating the combined effects of both source-detector separation and optical coefficient value in a single parameter for both scattering (${\mu }_s^{\prime }\rho$) and absorption (${\mu }_a\rho$). These values are referred to as dimensionless scattering and dimensionless absorption, respectively. Correspondingly, relative reflectance is presented in a dimensionless form as dimensionless reflectance ($R{\rho }^2$). This makes it possible to investigate trends in reflectance relationships in terms of dimensionless variables without any loss of generality.

Figure 8(a) illustrates the dimensionless relationship between reflectance and scattering for all simulated values of $\rho$ between 250 μm and 1 cm. Figure 8(b) illustrates the relationship between normalized reflectance (which is not explicitly presented as dimensionless since it represents the ratio of two dimensionless reflectance values such that the ${\rho }^2$ values cancel) and dimensionless absorption; in this plot, data are provided for the case where dimensionless scattering (${\mu }_s^{\prime }\rho$) is equal to 0.25. (Note that the relationship between normalized reflectance and dimensionless absorption is unique for each value of dimensionless scattering.) Given the values of ${\mu }_s^{\prime }$ and $\rho$ used to generate the simulation data, only a limited number of combinations of ${\mu }_s^{\prime }$ and $\rho$ yield a dimensionless scattering value of 0.25. Thus, the data in Fig. 8(b) are provided for only $\rho =250\mu \text{m}$, 500 μm, 1 mm, and 5 mm. The choice of presenting reflectance, ${\mu }_s^{\prime }$, and ${\mu }_a$ in dimensionless units is for simplicity of normalization and to show consistency of trends over a large range of $\rho$ values. It is important to note that when plotted for a fixed value of $\rho$ using nondimensionless units, the shape of the curve remains the same; the dimensionless normalization simply scales the reflectance values.

Fig. 8

Dimensionless reflectance relationships for a range of source-detector separations. (a) Dimensionless reflectance versus dimensionless scattering. (b) Normalized reflectance versus dimensionless absorption for a dimensionless scattering value (${\mu }_s^{\prime }\rho$) of 0.25. All data are for a phase function gamma value ($\gamma$) of 2.3.

The results presented in Fig. 8 also enable several important observations about reflectance relationships, illustrating the significant variations in reflectance over large dimensionless scattering and absorption ranges, especially for the reflectance versus scattering curve. In contrast to the nearly linear relation of relative reflectance versus scattering illustrated in Figs. 4 and 5, the plot in Fig. 8(a) demonstrates considerable nonlinearity. For low dimensionless scattering, the dimensionless reflectance increases as both scattering and source detector separation increase. But at a given distance, dimensionless reflectance reverses its relationship with ${\mu }_s^{\prime }\rho$ and begins to decrease as ${\mu }_s^{\prime }\rho$ increases.

The observation of this scattering-dependent reflectance relationship is a key component of the generalized reflectance formalism. Several authors have observed this behavior experimentally and it is also predicted by diffusion theory.^25,51,52 Further, Kumar and Schmitt related the curve peak to an ideal distance at which reflectance is only slightly dependent on ${\mu }_s^{\prime }$, which they reported to be between 2 and 5 mm, depending on the range of ${\mu }_s^{\prime }$ values of interest.⁵¹ For a typical tissue value of ${\mu }_s^{\prime }$ of $10{\mathrm{cm}}^{-1}$, this range of separation distances corresponds to dimensionless scattering values between 2 and 5, which correlates well with the peak of the dimensionless reflectance versus dimensionless scattering curve. It is at this peak value that moderate changes in ${\mu }_s^{\prime }$ will result in only minimal change in reflectance, suggesting that reflectance is insensitive to ${\mu }_s^{\prime }$ at source detector separations between 2 and 5 mm.

A novel and significant observation results from extending this analysis to consider variations in phase function. Figure 9 illustrates the dimensionless reflectance versus dimensionless scattering relationships for a range of phase function $\gamma$ values (1.2, 1.75, and 2.3). The critical observation in this figure is the existence of an isosbestic point, where reflectance is not dependent on $\gamma$ at a dimensionless scattering value of $\sim 0.7$. (Note that we adopt the term “isosbestic” loosely, with the intent of describing a condition for which one measurement parameter is independent of another.) While reflectance is known to be insensitive to phase function in the diffusion regime (${\mu }_s^{\prime }\rho >10$), it is an unexpected revelation that there is an additional singular dimensionless scattering distance, well below the diffusion regime, at which this insensitivity to phase function also exists. Because this isosbestic point is part of the generalized reflectance model, it is applicable to any turbid medium (i.e., independent of tissue model), and is not dependent on the fiber probe (including fiber separation, which is inherently part of dimensionless scattering).

Fig. 9

Dimensionless reflectance versus dimensionless scattering for a range of phase function $\gamma$ values, plus diffusion theory. The three curves with varying $\gamma$ values all intersect at a ${\mu }_s^{\prime }\rho$ value of 0.7, although diffusion theory does not.

For low dimensionless scattering values, a high $\gamma$ corresponds to lower reflectance, as was observed in Figs. 4 and 5 (where ${\mu }_s^{\prime }\rho$ ranged from 0.125 to 0.75). As dimensionless scattering increases, however, a point is reached at which reflectance is insensitive to the value of $\gamma$. And at yet larger dimensionless scattering values, a larger value of $\gamma$ results in a larger reflectance. As noted earlier, we hypothesize that this is because, in tissues with low $\gamma$, photons with more backward scattering and less forward scattering are more likely to remain close to the source. With high $\gamma$, more forward scattering and fewer backscattering events allow a photon to travel farther from the source before experiencing the high-angle scattering event that will direct it back toward the surface for collection. At large scattering and source-detector separations, this means that the proportion of photons collected at larger source-detector separations to photons collected at smaller separations is higher for large $\gamma$, leading to higher reflectance.

For comparison, the diffusion theory equation was used to calculate dimensionless reflectance as a function of dimensionless scattering, and is also included in Fig. 9 (dashed line). Recalling that diffusion theory is not dependent on phase function, above the isosbestic point the diffusion theory plot best approximates the $\gamma =1.2$ curve, and below the isosbestic point diffusion theory best approximates the $\gamma =2.3$ curve. This same relationship between diffusion theory and Monte ○Carlo simulations with different anisotropy $g$ values was recently observed by Kim et al. for values of ${\mu }_s^{\prime }$ and $\rho$ corresponding to dimensionless scattering from 0 to 1,⁵² which is consistent with our findings. At large dimensionless scattering values, the similarity between diffusion theory and the $\gamma =1.2$ plot is consistent with the fact that diffusion theory is descriptive of isotropic scattering; since the value $\gamma =1.2$ corresponds to $g=0.55$, the low $\gamma$ curve (in this range) is closest to the isotropic case ($g=0$ and $\gamma =1$).

We now plot normalized reflectance versus dimensionless scattering, which provides different observations compared with the previously illustrated dimensionless reflectance versus dimensionless scattering relationships in Figs. 8 and 9. Figure 10(a) presents the relationships over a large range of dimensionless absorption values from 0 to 10, and for five values of dimensionless scattering between 0.1 and 8.0. We note that as dimensionless scattering increases, the plots relating normalized reflectance and dimensionless absorption do not monotonically progress. For values of dimensionless scattering below $\sim 2.5$, normalized reflectance is larger for larger values of ${\mu }_s^{\prime }\rho$, as was seen in the relationships illustrated in Figs. 4 and 5. However, for dimensionless scattering values $>2.5$, normalized reflectance is smaller for larger values of ${\mu }_s^{\prime }\rho$. To better visualize this nonmonotonic relation, Fig. 10(b) plots normalized reflectance as a function of dimensionless scattering for a single value of dimensionless absorption, ${\mu }_a\rho =4.0$.

Fig. 10

(a) Normalized reflectance versus dimensionless absorption for five dimensionless scattering values between 0.1 and 8.0. (b) Normalized reflectance versus dimensionless scattering for a constant value of dimensionless absorption of 4.0.

This relationship was explored by Mourant et al., who reported on an optimized source-detector separation for minimizing the effects of scattering.⁵³ In contrast to the analysis of Kumar and Schmitt,⁵¹ who identified a source-detector separation that reduced variance in reflectance, Mourant et al. examined the source-detector separation at which variation in photon path length is minimized (enabling application of Beer’s law for determining the absorption coefficient), which they experimentally determined to be $\sim 1.7\mathrm{mm}$. That observation is supported by the results presented in Fig. 10(b), which illustrates that for dimensionless scattering values between $\sim 2$ and 3, the variation in normalized reflectance (which is directly related to path length) is minimal. When considering average values of ${\mu }_s^{\prime }$ in tissue, this range of dimensionless scattering values is consistent with the reported 1.7 mm source-detector separation distance for which normalized reflectance is independent of scattering (at a given ${\mu }_a$).

Like the dimensionless reflectance versus dimensionless scattering relationships, the normalized reflectance versus dimensionless absorption relationships are also dependent on phase function. Figure 11 presents the percent difference between normalized reflectance values with a phase function value $\gamma =1$.2 and normalized reflectance values with a phase function value $\gamma =2.3$. The percent difference is illustrated for a range of dimensionless scattering values between 0.1 and 8.0. As expected, when absorption is zero, the error in normalized reflectance is zero, since, by definition, normalized reflectance is equal to 1.0 when ${\mu }_a=0$. When dimensionless scattering is low (corresponding to low ${\mu }_s^{\prime }$ and/or $\rho$), the error in normalized reflectance due to phase function is maximal, reaching as much as 110%. As dimensionless scattering increases and approaches the diffuse regime, the error in reflectance due to phase function variation reduces significantly. However, even at a dimensionless scattering value of 8.0, the maximal reflectance error due to phase function is still over 20%, which is consistent with the known breakdown of the diffusion approximation for large values of absorption. The importance of this observation is that the phase function is an important variable for reflectance even at relatively large source-detector separations.

Fig. 11

Percent difference in normalized reflectance values between simulations when the phase function parameter $\gamma =1.2$ and $\gamma =2.3$. Percent differences are illustrated for a range of dimensionless scattering values between 0.1 and 8.0.

4.3.

Extracting Phase Function Information

While extracting ${\mu }_s^{\prime }$ and ${\mu }_a$ values from tissue has generally been of primary interest, the variation in reflectance as a function of $\gamma$ presents an encouraging opportunity to extract additional information about the tissue using DRS. Given that $\gamma$ relates to the relative backscattering contributions of a tissue phase function, which is dictated by the size distribution of scatterers, the ability to quantify $\gamma$ can provide knowledge about the tissue microstructure. A handful of other authors have also examined this opportunity^3,13,54,55 for other measurement geometries (e.g., single-fiber probes that both illuminate and collect light) and have observed that there is not enough information in a single reflectance spectrum to uniquely identify ${\mu }_s^{\prime }$, ${\mu }_a$, and phase function information simultaneously, and therefore, multiple measurements are necessary. In particular, Kanick et.al. have demonstrated that the value of $\gamma$ can be extracted from multiple reflectance measurements using multiple single-fiber probes of different diameters.³⁵ Extending this to the two-fiber collection geometry, collecting measurements at various source-detector separations is a natural choice, especially since the reflectance relationships vary significantly as a function of $\rho$.

The implementation of this approach involves measuring reflectance spectra at multiple distances and simultaneously fitting the multiple reflectance spectra to a reflectance model in two dimensions, with both $\rho$ and $\lambda$ being independent variables. In addition to the original inverse fitting coefficients described in Sec. 2.1.2, $\gamma$ then becomes the sixth fitting coefficient. The reflectance model to which the data are fit is a set of three-dimensional lookup tables for which relative reflectance is defined as a function of ${\mu }_a$, ${\mu }_s^{\prime }$, and $\gamma$. Each three-dimensional lookup table is unique to a single source-detector separation geometry. Experimental reflectance spectra are fit to their respective probe geometry reflectance lookup tables, as would be done with only one reflectance measurement being analyzed, but a least-squares algorithm performs iterative calculations on all data simultaneously, optimizing the physiological fitting coefficients, including $\gamma$, for the entire set of data.

Using the ${\mu }_s^{\prime }$ and ${\mu }_a$ spectra described in Sec. 3.2, reflectance spectra were constructed for a sample value of $\gamma$ and for a number of different source-detector separations. Various pairs of reflectance spectra, representing different pairs of source-detector separations, were then analyzed in the inverse direction using the above lookup table approach. In all cases, the extracted optical coefficient values, including the value of $\gamma$, were accurately retrieved with errors $<0.25\%$. This indicates that with only two separate reflectance spectra and with the differences in source-detector separations as small as 250 μm (when the two reflectance spectra are measured using $\rho =250$ and 500 μm), not only can errors in the extracted optical properties due to uncertainty in phase function be significantly reduced, but an additional optical parameter can be gained (that of the phase function parameter gamma, $\gamma$), which can serve to better characterize and differentiate tissue types.