1.

Introduction

Reflectance spectroscopy techniques are used to relate tissue optical properties to various types of diseases. Depending on the clinical question, an important consideration in choosing a spectroscopic technique is its sampling volume. Since optical properties are averaged over this volume, techniques with small sampling volumes are required to detect local, small-scale changes in tissue. Superficial sampling is relevant for early detection of cancer in general, and especially for the detection of (premalignant) lesions in the epithelium of, e.g., the esophagus or the colon, since the epithelium is only a few hundred micrometers thick. Single fiber reflectance (SFR) spectroscopy is a technique with a sampling depth in the order of a few hundred micrometers and has been used in a number of studies in the field of oncology.^1–6 SFR spectroscopy uses a single fiber, in contact with tissue, to emit and collect light, connected to a broadband light source and a spectrograph to detect the steady-state reflectance versus wavelength. From the measured spectra, tissue optical properties can be derived and related to the disease state of tissue.

Reflectance in terms of optical properties can be expressed in general as

Currently, a single model is available to derive optical properties from SFR measurements, developed by Kanick et al.^7,8 This model consists of two parts: one part describes reflectance as a function of tissue-scattering properties in the absence of absorption ($R_0$)⁷ and the second part includes absorption.⁸ Kanick et al. simplified the modeling by approximating Eq. (1) by Eq. (2), describing reflectance as reflectance in the absence of absorption ($R_0$) multiplied by a factor to account for absorption, featuring ${\mu }_a$ and an effective pathlength $⟨L⟩$:

The model for $R_0$ describes the reflectance as a function of the product of the reduced scattering coefficient and the fiber diameter (${\mu }_s^{\prime }d$), the phase function parameter $\gamma$ (which depends on the first and second Legendre moments of the phase function), the fiber numerical aperture (NA) and the refractive index of the sample ($n_{\text{sample}}$):

This model was derived using Monte Carlo (MC) simulations of tissues with modified Henyey–Greenstein (MHG) phase functions. We previously demonstrated that Eq. (3) does not predict the reflectance well for tissues with phase functions other than MHG.⁹

For many tissue types, it has been shown that scattering is more backward directed than can be described by the MHG. The two-term Henyey–Greenstein (TTHG) has been shown to describe the phase function of tissue in, e.g., the liver,^10,11 uterus,¹⁰ brain,¹² breast,¹³ and muscle,¹⁴ and for blood, it has been shown that the phase function is better represented by the Reynolds–McCormick (RMC, also known as Gegenbauer) phase function.^15,16 The phase function of a specific tissue under investigation is generally not known. For example, even though a TTHG phase function has been measured in the breast, this does not imply that breast tissue will always have a TTHG phase function. The phase function will depend on the type of breast tissue that is investigated (e.g., fat or ducts) and it might change during progression from healthy to diseased tissue. Therefore, to ensure an accurate determination of the optical properties from SFR measurements, a model that is valid for the wide range of phase functions that can be encountered in tissue is essential.

In this paper, as a first step toward a comprehensive SFR spectroscopy model that includes both scattering and absorption properties, we develop a model that accurately relates the reflectance in the absence of absorption ($R_0$) to tissue-scattering properties for a wide range of tissue phase functions. This model can be expanded in the future to also include absorption. We model $R_0$ as a combination of a diffuse and a semiballistic component. For SFR spectroscopy, diffusion theory alone is not appropriate to model reflectance as a function of tissue optical properties. Since a single fiber is both the source and detector, the distance between the location where photons enter the tissue and where they are detected is generally less than the transport mean free path $1/{\mu }_s^{\prime }$. SFR measurements are, therefore, in the so-called subdiffuse regime, where the measured reflectance is a combination of detected photons that underwent a large number of scattering events (diffuse photons) and detected photons that underwent only a few scattering events (semiballistic photons). We consider the detected photons semiballistic when they have experienced a single backscattering event in combination with an arbitrary number of forward-scattering events.

We describe the diffuse component using the model of Farrell et al.¹⁷ for spatially resolved diffuse reflectance, as solved by Faber et al. for an overlapping source and detector fiber.¹⁸ The diffuse component of the reflectance depends only on the product ${\mu }_s^{\prime }d$, where ${\mu }_s^{\prime }={\mu }_s(1-g_1)$, ${\mu }_s$ is the scattering coefficient and $g_1$ is the scattering anisotropy (first Legendre moment of the phase function). After a few scattering events, photon direction is randomized. Therefore, diffuse reflectance does not depend on the details of the tissue phase function but only on the scattering anisotropy $g_1$. For semiballistic photons, the photon direction is not fully randomized and, thus, measurements are more sensitive to the shape of the phase function.^19,20 Therefore, we searched for a parameter that would optimally capture the phase function influence on the semiballistic contribution to the reflectance. We developed the parameter $p_{\mathrm{sb}}$ and showed that it improves the prediction of the semiballistic contribution to the reflectance compared to other parameters that have been proposed to model subdiffuse reflectance ($\gamma$,²¹ $\delta$,²² $\sigma$,²³ and $R_{p\mathrm{NA}}$⁹). Therefore, we model the semiballistic contribution to the reflectance as a function of ${\mu }_s^{\prime }d$ and $p_{\mathrm{sb}}$. Based on MC simulations with MHG, TTHG, and RMC phase functions, we show that our model (combining the diffuse and semiballistic contributions) predicts the reflectance more accurately, compared to the currently used model from Kanick et al.⁷ for a range of phase functions, ${\mu }_s^{\prime }d$ values and NAs.

2.

Methods

2.2.

Scattering Model

We model the reflectance ($R_0$) as the sum of a diffuse ($R_{\mathrm{SFR},\mathrm{dif}}$) reflectance and a semiballistic ($R_{\mathrm{SFR},\mathrm{sb}}$) reflectance:

Eq. (4)

$$R_0=R_{\mathrm{SFR},\mathrm{dif}}+R_{\mathrm{SFR},\mathrm{sb}.}$$

The ratio between the contribution of semiballistic and diffuse photons will depend on the scattering properties of the tissue. For example, for high values of ${\mu }_s^{\prime }d$, almost all detected photons will be diffuse, whereas for low values of ${\mu }_s^{\prime }d$, most detected photons will be semiballistic. Therefore, we rewrite Eq. (4) using this ratio

Eq. (5)

$$X=\frac{R_{\mathrm{SFR},\mathrm{sb}}}{R_{\mathrm{SFR},\mathrm{dif}}}$$

as

Eq. (6)

$$R_0=(1+X)·R_{\mathrm{SFR},\mathrm{dif}}.$$

2.2.1.

Diffuse component

We describe the diffuse component, $R_{\mathrm{dif}}$, using the model of Farrell et al.¹⁷ for spatially resolved reflectance as solved by Faber et al.¹⁸ for an overlapping source and detection fiber in the absence of absorption:

Eq. (7)

$$R_{\mathrm{dif}}({\mu }_s^{\prime }d)=\frac{\pi }{4}·d^2·\underset{0}{\overset{d}{\int }}{R}_{\mu a=0}(\rho ,{\mu }_s^{\prime })·p(\rho ,d)d\rho \mathrm{.}$$

Here, $R_{\mu a=0}(\rho ,{\mu }_s^{\prime })$ is the diffuse reflectance versus radial distance for a pencil beam illumination,

Eq. (8)

$$R_{\mu a=0}(\rho ,{\mu }_s^{\prime })=\frac{1}{4\pi }\{\frac{{{\mu }_s^{\prime }}^2}{{[1+{({\mu }_s^{\prime }·\rho )}^2]}^{\frac{3}{2}}}+\frac{(1+\frac{4A}{3})·{{\mu }_s^{\prime }}^2}{{[{(1+\frac{4A}{3})}^2+{({\mu }_s^{\prime }·\rho )}^2]}^{\frac{3}{2}}}\}$$

and $p(\rho ,d)$ is the distribution function of distances over the fiber face,

Eq. (9)

$$p(\rho ,d)=\frac{16\rho }{\pi d^2}{\cos }^{-1}\left(\frac{\rho }{d}\right)-\frac{16}{\pi d}{\left(\frac{\rho }{d}\right)}^2\sqrt{1-{\left(\frac{\rho }{d}\right)}^2.}$$

Equation (9) describes the distribution of distances between two randomly placed points on a disk with diameter $d$, which is a classic problem in the field of geometric probability.²⁷ The parameter $A$ depends on the refractive index mismatch between the tissue and the fiber.²⁸ For a tissue refractive index of 1.35 and a fiber refractive index of 1.45, $A$ is equal to 1.027.

Using Eqs. (7)–(9), $R_{\mathrm{dif}}$ models the fraction of photons that are diffuse and arrive at the fiber face. Photons are only detected if they arrive at the fiber at an angle smaller than or equal to the acceptance angle of the fiber. To account for the acceptance angle of the fiber, we need to multiply $R_{\mathrm{dif}}$ by the collection efficiency of the fiber (${\eta }_c$):

Eq. (10)

$$R_{\mathrm{SFR},\mathrm{dif}}={\eta }_c·R_{\mathrm{dif}}({\mu }_s^{\prime }d).$$

The collection efficiency of the fiber depends on the angular distribution of photons reaching the fiber face and the acceptance angle of the fiber. If the reflected light has a Lambertian profile, then ${\eta }_{c,L}$ is equal to²⁹

Eq. (11)

$${\eta }_{c,L}=\sin {\theta }_{\mathrm{acc}}^2.$$

Since it is not known what the angular profile is for SFR measurements, we account for a possible non-Lambertian profile of photons reaching the fiber by modeling the collection efficiency as the Lambertian collection efficiency multiplied by a constant,

Eq. (12)

$${\eta }_c=\sin {\theta }_{\mathrm{acc}}^2·a_1.$$

An important check of our model is whether the fit parameter $a_1$ accurately relates to the collection efficiency. Therefore, we determined the reflectance values for simulations with a high value of ${\mu }_s^{\prime }d$ (1000; ${\mu }_s^{\prime }={10}^5{\mathrm{cm}}^{-1}$, $d=100\mu \mathrm{m}$) and compared these to the collection efficiency ${\eta }_c$. For high values of ${\mu }_s^{\prime }d$, the reflectance will be diffuse and, thus, $R_{\mathrm{SFR},\mathrm{sb}}$ will approach 0 and $R_{\mathrm{dif}}$ will approach 1. Therefore, the reflectance should equal to the collection efficiency for high values of ${\mu }_s^{\prime }d$.

2.2.2.

Semiballistic component

For the semiballistic contribution to the reflectance, we searched for a parameter that would optimally capture the influence of the phase function. Usually, such attempts are based on the development of the phase function in Legendre polynomials:

Eq. (13)

$$p(\cos \theta )=\frac{1}{4\pi }\sum _{\mathrm{n}}(2n+1)g_n{P}_n(\cos \theta ),$$

where $p$ is the probability of scattering at an angle $\theta$, $P_n$ are the Legendre polynomials of order $N$, and $g_n$ are the Legendre moments. The scattering model from Kanick et al.⁷ describes the influence of the phase function by incorporating the parameter $\gamma =(1-g_2)/(1-g_1)$. Next to $\gamma$, three other parameters have been proposed that could be used to incorporate the phase function influence in models to describe subdiffuse reflectance: $\delta$ ,²² $\sigma$,²³ and $R_{p\mathrm{NA}}$.⁹ The parameters $\delta$ and $\sigma$ also incorporate moments of the phase function: $\delta =(1-g_3)/(1-g_1)$ and $\sigma =\sum _{i=2}^{\infty }{(-0.5)}^{i-2}\frac{1-g_i}{1-g_1}$.

We recently introduced the parameter $R_{p\mathrm{NA}}$,⁹ to describe the semiballistic contribution to the reflectance for SFR spectroscopy. Photons were considered semiballistic if they had experienced a single backscatter event in combination with an arbitrary number of forward-scattering events upon detection. For these semiballistic photons, it was assumed that they would only be detected if all scattering events occurred at $\text{scattering angles}\le {\theta }_{\mathrm{acc}}$. The full derivation is described in our previous paper.⁹ In short, based on this assumption, $R_{p\mathrm{NA}}$ is defined as

Eq. (14)

$$R_{p\mathrm{NA}}=\frac{p_{\mathrm{NA},\mathrm{b}}}{1-p_{\mathrm{NA},f}}.$$

Here, $p_{\mathrm{NA},b}$ is the probability of a scattering event to occur within ${\theta }_{\mathrm{acc}}$ in the backward direction, which equals the integral of the phase function over the acceptance angle in the backward direction:

Eq. (15)

$$p_{\mathrm{NA},b}=2\pi \underset{\pi -{\theta }_{\mathrm{acc}}}{\overset{\pi }{\int }}p(\theta )\sin (\theta )\mathit{d}\theta .$$

Here $p_{\mathrm{NA},f}$ is the probability of a scattering event to occur within ${\theta }_{\mathrm{acc}}$ in the forward direction:

Eq. (16)

$$p_{\mathrm{NA},f}=2\pi \underset{0}{\overset{{\theta }_{\mathrm{acc}}}{\int }}p(\theta )\sin (\theta )\mathit{d}\theta .$$

We showed previously that for SFR spectroscopy the parameter $R_{p\mathrm{NA}}$ predicts the reflectance better than $\gamma$ or $\sigma$,⁹ but we did not yet validate the assumptions behind the derivation of $R_{p\mathrm{NA}}$.

To test these assumptions here, we investigated the scattering angles of semiballistic photons that are detected. To this end, we examined the scattering angles of detected photons for the simulations with ${\mu }_s^{\prime }d=0.1$ and NAs of 0.10, 0.22, and 0.50. We extracted the normalized frequency distribution of scattering angles of detected photons from the MC simulations. We call this the “effective phase function.” Based on the results of the effective phase functions, we searched for a parameter that could better describe the influence of the phase function on the measured reflectance, by investigating the optimal integration limits for a parameter similar to $R_{p\mathrm{NA}}$, which we named $R_{p(\theta b,\theta f)}$:

Eq. (17)

$$R_{p(\theta b,\theta f)}=\frac{p_b({\theta }_b)}{1-p_f({\theta }_f)},$$

Eq. (18)

$$p_b({\theta }_b)=2\pi \underset{\pi -{\theta }_b}{\overset{\pi }{\int }}p(\theta )\sin (\theta )\mathit{d}\theta ,$$

Eq. (19)

$$p_f({\theta }_f)=2\pi \underset{0}{\overset{{\theta }_f}{\int }}p(\theta )\sin (\theta )\mathit{d}\theta .$$

We investigated the relation between $R_0$ and $R_{p(\theta b,\theta f)}$ for all possible backward and forward integration angles (in steps of 1 deg). Optimal integration angles were chosen based on a minimization of the dispersion in the parameter $\log 10[R_{p(\theta b,\theta f)}]$ for a chosen reflectance, relative to the total range of $\log 10[R_{p(\theta b,\theta f)}]$ (Fig. 1). For each NA, we determined the relative dispersion for five different reflectance values ($\pm 10\%$) equally spaced on a logarithmic scale. This relative dispersion is a measure for the suitability of $R_{p(\theta b,\theta f)}$ to model the reflectance, because the dispersion in the parameter will influence the uncertainty in all the estimated optical properties when $R_{p(\theta b,\theta f)}$ is used to model the reflectance. A relative measure is used because redefining the parameter by multiplying or subtracting a value would otherwise influence the measure.

Fig. 1

Concept of judging the suitability of the parameter $R_{p(\theta b,\theta f)}$ for a specific forward and backward integration angle. For a certain reflectance value ($R_0$), the dispersion of $R_{p(\theta b,\theta f)}$ for that reflectance (light blue) is divided by the entire range of $R_{p(\theta b,\theta f)}$ (green). This relative dispersion was determined for five reflectance values (depicted in red) for ${\mu }_s^{\prime }d=0.1$, per set of 207 simulations with the same NA.

Based on these results, we defined the parameter $p_{\mathrm{sb}}$ as

Eq. (20)

$$p_{\mathrm{sb}}=\frac{p_b(1\mathrm{deg})}{1-p_f(23\mathrm{deg})}.$$

We compared $p_{\mathrm{sb}}$ to $\gamma$, $\sigma$, $\delta$, and $R_{p\mathrm{NA}}$ by determining the relative dispersion for three different reflectance values ($\pm 10\%$) per NA. Since we developed $p_{\mathrm{sb}}$ for the semiballistic contribution to the reflectance, we compared these parameters for ${\mu }_s^{\prime }d=0.1$, where most detected photons are semiballistic.

2.2.3.

Ratio diffuse and semiballistic photons

We use $p_{\mathrm{sb}}$ in the model for the ratio between semiballistic and diffuse photons $X$. For low values of ${\mu }_s^{\prime }d$, the diffuse reflectance scales with ${({\mu }_s^{\prime }d)}^2$. Therefore, we include that factor in the denominator of Eq. (21). The semiballistic reflectance will scale differently with ${\mu }_s^{\prime }d$. We assume $R_{\mathrm{SFR},\mathrm{sb}}$ scales with $p_{\mathrm{sb}}{({\mu }_s^{\prime }d)}^{-B}$. This difference in scaling is accounted for through the parameter $a_3$.

Eq. (21)

$$X=a_2{\left[\frac{p_{sb}}{{({\mu }_s^{\prime }d)}^2}\right]}^{a_3}.$$

The parameters $a_1-a_3$ were determined by fitting the model [Eqs. (5)–(10), (12), (20)–(21)] to the MC results.

3.

Results

3.1.

Phase Function Parameter Semiballistic Contribution ($p_{\mathrm{sb}}$)

For ${\mu }_s^{\prime }d=0.1$ and NAs of 0.10, 0.22, and 0.50, the contribution of semiballistic photons (with only a single backscattering event and any number of forward-scattering events) to the total reflectance was 99%, 98%, and 98%, respectively. Therefore, we used these simulations to examine the semiballistic contribution to the reflectance.

To investigate the assumption behind $R_{p\mathrm{NA}}$ that semiballistic photons are only detected if all the scattering events occurred at $\text{scattering angles}\le {\theta }_{\mathrm{acc}}$, we determined the effective phase functions (the frequency distributions of the scattering angles of detected photons) for each simulation. From these effective phase functions, it can be concluded that detected photons also underwent scattering angles outside ${\theta }_{\mathrm{acc}}$ (Fig. 2). As depicted in Fig. 2(b), the NA has almost no influence on the effective phase function for angles in the forward direction. The NA does have an influence on the percentage of detected scattering events for angles in the backward direction. The sharp increase corresponds to twice ${\theta }_{\mathrm{acc}}$, which can be explained by the fact that photons are launched with a random angle within the acceptance angle and that photons are only detected when they arrive at the fiber with an angle within the acceptance angle of the fiber. A photon that is launched at ${\theta }_{\mathrm{acc}}$ and undergoes only a single scattering event can be detected if that scattering angle is $2·{\theta }_{\mathrm{acc}}$.

Fig. 2

(a) Effective phase function: frequency distribution of scattering angles of detected photons for all 207 simulations with ${\mu }_s^{\prime }d=0.1$ (${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$, $d=100\mu \mathrm{m}$) and $\mathrm{NA}=0.22$. For each simulation, the frequency distribution is normalized to 1. The gray lines represent individual effective phase functions, the black line indicates the average effective phase function, and the red dashed lines indicate the acceptance angle of the fiber, ${\theta }_{\mathrm{acc}}$. (b) Average effective phase functions for ${\mu }_s^{\prime }d=0.1$ (${\mu }_s^{\prime }=10{\mathrm{cm}}^{-1}$, $d=100\mu \mathrm{m}$) and $\mathrm{NA}=0.10$, 0.22, or 0.50. The dashed lines indicate ${\theta }_{\mathrm{acc}}$ for each NA. Detected photons also underwent scattering angles outside ${\theta }_{\mathrm{acc}}$.

Since detected photons also underwent scattering events at angles larger than ${\theta }_{\mathrm{acc}}$, we searched for a parameter that could better describe the influence of the phase function on the measured reflectance, by investigating the optimal integration limits for $R_{p(\theta b,\theta f)}$ [Eqs. (17)–(19)]. The optimization of the integration angles was based on minimizing the relative dispersion in the parameter $\log 10[R_{p(\theta b,\theta f)}]$. Figure 3(a) depicts the relative dispersion versus backward integration angle, for the optimal forward integration angle per NA. The relative dispersion is fairly constant for integration angles up to 20 deg and the results are similar for different NAs. Figure 3(b) depicts the relative dispersion versus forward integration angle, for the optimal backward integration angle per NA. There is a sharp optimum for the forward integration angle and the optimum integration angle varies with the NA. For each NA, the optimal integration angles can be found in Table 2.

Fig. 3

Relative dispersion versus integration angle for simulations with NAs of 0.10, 0.22, and 0.50, ${\mu }_s^{\prime }d=0.1$ (a) Relative dispersion versus backward integration angle, for the optimal forward integration angles (16 deg, 23 deg, and 23 deg, respectively). (b) Relative dispersion versus forward integration angle, for the optimal backward integration angle (1 deg for all three NAs). The relative dispersion is much more sensitive to the forward integration angle than the backward integration angle and is slightly influenced by the NA.

Table 2

Optimal integration angles backward (θb) and forward (θf) for different values of the NA.

NA	θb	θf
0.10	1	16
0.22	1	23
0.50	1	23

Since fibers with an NA of 0.22 are generally used for SFR measurements, we optimized our parameter for this NA. The optimal integration angle that results in the lowest relative dispersion is 1 deg backward and 23 deg forward. Therefore, we will use $p_{\mathrm{sb}}$ [Eq. (20)] in our model for the semiballistic contribution to the reflectance.

Figure 4 shows the relation between the reflectance and $\gamma$, $\sigma$, $\delta$, $R_{p\mathrm{NA}}$, and $p_{\mathrm{sb}}$ for simulations with an NA of 0.22 and ${\mu }_s^{\prime }d=0.1$. Overall, $p_{\mathrm{sb}}$ has a lower relative dispersion. A quantitative measure of the dispersion was calculated by determining the relative dispersion of each parameter for three reflectance values per NA (Table 3). For an NA of 0.10, $p_{\mathrm{sb}}$ outperforms the other parameters for all three reflectance values. For an NA of 0.22, $p_{\mathrm{sb}}$ outperforms the other parameters for the lower three reflectance values but has a higher relative dispersion for the highest reflectance value. For an NA of 0.50, the results for $R_{p\mathrm{NA}}$ and $p_{\mathrm{sb}}$ are similar (and both outperform the other parameters), which is expected since the integration angles for $R_{p\mathrm{NA}}$ are 22 deg for an NA of 0.50.

Fig. 4

Simulated reflectance ($R_0$) versus $\gamma$, $\sigma$, $\delta$, $R_{p\mathrm{NA}}$, and $p_{\mathrm{sb}}$, for an NA of 0.22 and ${\mu }_s^{\prime }d=0.1$; colors indicate phase function types. Note the log scales for both the reflectance, $R_{p\mathrm{NA}}$, and $p_{\mathrm{sb}}$. Here $p_{\mathrm{sb}}$ is most closely related to the semiballistic reflectance.

Table 3

Variability of psb, RpNA, σ, δ, and γ for μs′d=0.1, defined as the dispersion in psb, RpNA, σ, and γ values for a chosen reflectance (+/- 10%) relative to the total range of each parameter.

NA	Reflectance	Variability
		psb	RpNA	σ	δ	γ
0.10	0.0001	0.013	0.030	0.042	0.178	0.074
	0.0002	0.045	0.104	0.092	0.344	0.183
	0.0004	0.087	0.153	0.138	0.322	0.226
0.22	0.0005	0.016	0.029	0.050	0.231	0.111
	0.001	0.045	0.060	0.114	0.272	0.185
	0.003	0.217	0.130	0.096	0.095	0.086
0.50	0.001	0.002	0.001	0.015	0.131	0.056
	0.003	0.030	0.035	0.089	0.344	0.185
	0.005	0.088	0.115	0.141	0.322	0.226

3.2.

Scattering Model

Table 4 lists the fit parameters obtained from fitting Eqs. (5)–(10), (12), and (20)–(21) to all the MC simulations. For all three NAs, the value of $a_1$ is similar. Therefore, we chose a single value for $a_1$ for all three NAs and repeated the fit to obtain values for $a_2$ and $a_3$ for each NA. This resulted in a new set of parameter values to be used in the SFR scattering model (Table 5).

Table 4

Resulting fit parameters per NA using Eqs. (5)–(10), (12), (20), and (21) on simulated SFR measurements. The 95% confidence intervals on these fit parameters are indicated.

	NA=0.10		NA=0.22		NA=0.50
	Value	95% CI	Value	95% CI	Value	95% CI
$a_1$	1.130	($\pm 0.006$)	1.119	($\pm 0.005$)	1.098	($\pm 0.005$)
$a_2$	4427	($\pm 87$)	3065	($\pm 50$)	1461	($\pm 22$)
$a_3$	0.785	($\pm 0.003$)	0.750	($\pm 0.003$)	0.684	($\pm 0.002$)

Table 5

Resulting fit parameters per NA using Eqs. (5)–(10), (12), (20), and (21) and a1=1.11 on simulated SFR measurements. The 95% confidence intervals on these fit parameters are indicated.

	NA=0.10		NA=0.22		NA=0.50
	Value	95% CI	Value	95% CI	Value	95% CI
$a_2$	4370	($\pm 87$)	3046	($\pm 49$)	1475	($\pm 23$)
$a_3$	0.780	($\pm 0.003$)	0.748	($\pm 0.003$)	0.688	($\pm 0.003$)

To determine whether the fit parameter $a_1$ accurately relates to the collection efficiency, we compared the reflectance values for simulations with a high value of ${\mu }_s^{\prime }d$ (1000; ${\mu }_s^{\prime }={10}^5{\mathrm{cm}}^{-1}$, $d=0.01\mathrm{cm}$) to the collection efficiency using Eq. (12) in combination with the $a_1$ values from Table 4. For high values of ${\mu }_s^{\prime }$, the reflectance should equal the collection efficiency. The differences are 1.1%, 1.6%, and 0.7% for NAs of 0.10, 0.22, and 0.50, respectively.

Using the parameter $a_1$, we can now calculate the diffuse contribution to the reflectance ($R_{\mathrm{SFR},\mathrm{dif}}$) and show that the reflectance ($R_0$) is indeed a sum of a diffuse and a semiballistic component (Fig. 5). For high values of ${\mu }_s^{\prime }d$, the total reflectance equals the diffuse reflectance (dashed black line), but for lower values of ${\mu }_s^{\prime }d$ there is an additional semiballistic contribution to the reflectance. With increasing $p_{\mathrm{sb}}$, the fraction of semiballistic photons increases.

Fig. 5

Simulated $R_0$ as a function of ${\mu }_s^{\prime }d$. For lower values of ${\mu }_s^{\prime }d$, an additional semiballistic reflectance is added to the diffuse reflectance ($R_{\mathrm{SFR},dif}$). With increasing $p_{\mathrm{sb}}$, the fraction of semiballistic photons increases.

3.3.

Accuracy Scattering Model

To determine the accuracy of the scattering model, we used Eqs. (5)–(10), (12), and (20)–(21), in combination with the parameters from Table 5, to determine the difference between the reflectance predicted by the model to the reflectance from the MC simulations. First, we investigated the limit with respect to ${\mu }_s^{\prime }d$ for our model, by determining the difference for all the simulations with 10 different phase functions ($0.005<{\mu }_s^{\prime }d<{10}^3$). For values of ${\mu }_s^{\prime }d$ that are lower than 0.1, the difference increases significantly (Fig. 6).

Fig. 6

Absolute value of the difference (%) between the model and the MC simulations versus ${\mu }_s^{\prime }d$. For low values of ${\mu }_s^{\prime }d$, the difference increases significantly.

Next, we determined the difference between the reflectance predicted by the model and the reflectance from the MC simulations for all simulations where ${\mu }_s^{\prime }d\ge 0.1$ (Fig. 7). The median error is 2.1%, with a standard deviation of 3.0% and a maximum error of 16%. The error in the predicted reflectance is similar for all three NA values. The errors are much lower compared to the model from Kanick et al.⁷ [Fig. 7(b)], where the median error is 9.0%, with a standard deviation of 36.8% and a maximum error of 303%.

Fig. 7

Reflectance as predicted by the model ($R_{0,\text{model}}$) versus the reflectance obtained from the MC simulations ($R_{0,\text{simulations}}$). The black line depicts the perfect prediction. (a) For our model, the median error is 2.1% with a standard deviation of 3.0%. (b) For the model from Kanick et al.,⁷ the median error is 9.0% with a standard deviation of 36.8%.

3.4.

Monte Carlo Simulations Single Fiber Reflectance Spectra

Figure 8 shows the simulated spectra for Intralipid 20%, diluted to a volume fraction of 0.05. The residual of the fit is below 5% for the $300\text{-}\mu \mathrm{m}$ fiber and fluctuates around 0% for the $600\text{-}\mu \mathrm{m}$ fiber. The obtained ${\mu }_s^{\prime }$ and $p_{\mathrm{sb}}$ values are close to the values used in the simulations.

Fig. 8

(a) Fit on simulated spectra for Intralipid 20%, diluted to a volume fraction of 0.05. (b) Residual of fit. (c) Fit result for the reduced scattering coefficient. (d) Fit result for $p_{\mathrm{sb}}$.

4.

Discussion

SFR spectroscopy is a technique that has a sampling depth of a few hundred micrometers and has, therefore, the potential to detect local, small-scale changes in tissue—such as small or superficial tumors. However, the model currently used to obtain optical scattering properties from SFR measurements is limited to tissues with MHG phase functions. Here, we developed a new model that accurately relates the reflectance in the absence of absorption ($R_0$) to scattering properties of tissue for a wide range of phase functions, scattering coefficients, fiber diameters, and NAs. Our model predicts $R_0$ substantially better compared to the model from Kanick et al.⁷ (2.1% versus 9.0% median error). To determine optical properties of tissues with SFR spectroscopy, the next step will be to add absorption to the model.

We demonstrated the use of our model to extract optical properties from SFR spectra, using simulated spectra for two fibers of different diameters (MDSFR). Since there is currently no model available for the wavelength dependence of $p_{\mathrm{sb}}$, a measurement with two fibers of different diameters is required. The fitted values we obtained for ${\mu }_s^{\prime }$ are very close to the simulated ${\mu }_s^{\prime }$. The difference between the fitted and the simulated values of $p_{\mathrm{sb}}$ are slightly larger. A model for the wavelength-dependent behavior of $p_{\mathrm{sb}}$ could further simplify the measurements, since it would enable the use of only a single fiber.

We model $R_0$ as a sum of a diffuse and a semiballistic contribution. For the semiballistic contribution to the reflectance, we proposed a new parameter, $p_{\mathrm{sb}}$, to capture the influence of the phase function. Our findings imply that the semiballistic contribution is mainly influenced by the probability of photons scattering forward within 23 deg. A physical reason for this angular limit is yet to be found. The suitability of $R_{p(\theta b,\theta f)}$ to model the semiballistic contribution depends much more on the forward integration angle than the backward integration angle, even though this might not be expected based on our results for the effective phase functions. The large influence of the forward integration angle might be explained by the fact that semiballistic photons undergo only a single backscatter event, but generally undergo multiple forward-scattering events. Even though we had chosen $p_{\mathrm{sb}}$ to optimize the parameter for an NA of 0.22, the model performed similarly well for all three NAs. Previous work on subdiffuse scattering presented the subdiffuse parameters $\gamma$, $\delta$, $\sigma$, and $R_{p\mathrm{NA}}$. The parameters $\gamma$, $\delta$, and $\sigma$ are all based on the Legendre moments of the phase function. It can be concluded that the Legendre moments are not the most suitable for describing the semiballistic contribution to SFR measurements.

We model the diffuse contribution to the reflectance by solving the model for spatially resolved reflectance of Farrell et al.¹⁷ for an overlapping source and detection fiber.¹⁸ For low values of ${\mu }_s^{\prime }d$, the diffuse reflectance scales with ${({\mu }_s^{\prime }d)}^2$, which justifies the inclusion of that factor in the denominator of the ratio between semiballistic and diffuse photons. The semiballistic reflectance scaled differently with ${\mu }_s^{\prime }d$, which is accounted for in the parameter $a_3$.

As expected, the fraction of diffuse photons increases for higher values of ${\mu }_s^{\prime }d$ and lower values of $p_{\mathrm{sb}}$. Thus, the fiber diameter will influence the ratio of semiballistic to diffuse photons that are detected.

In our model, we model the collection efficiency as $\sim 1.11$ times the collection efficiency that would be expected if the light that reached the fiber face would have had a Lambertian profile. We determined that this indeed still represents the collection efficiency by showing that, for very high values of ${\mu }_s^{\prime }d$, the difference between the collection efficiency and the reflectance was around 1%. Another assumption for the collection efficiency could be that light reaching the fiber face has a flat angular distribution. For an NA of 0.22, this would result in a collection efficiency of 3.8%, compared to 2.7% for a Lambertian profile.³¹ In our model, we use 3% for this NA. This suggests that the angular profile of light reaching the fiber face is in between a Lambertian and a flat angular distribution.

Since part of our model is based on diffusion theory, it can be expected that there is a limit where ${\mu }_s^{\prime }d$ is too low for diffusion theory to accurately describe the diffuse contribution to the reflectance. For values of ${\mu }_s^{\prime }d$ below 0.1, our model becomes inaccurate. For tissues with low values of ${\mu }_s^{\prime }$, it is, therefore, advisable to use fibers with a larger diameter. For example, with a $500\text{-}\mu \mathrm{m}$ fiber, our model accurately predicts the reflectance for ${\mu }_s^{\prime }\ge 2{\mathrm{cm}}^{-1}$, and for most tissues ${\mu }_s^{\prime }$ will be larger than $2{\mathrm{cm}}^{-1}$.³²

For the derivation of our model, we have included a range of phase functions that have been measured in tissue. As far as we know, no other tissue phase functions have been measured yet, but new phase functions might be discovered in the future. In that case, it will have to be investigated whether the model is still accurate. Nevertheless, the tested phase functions are widely variable in shape, most importantly within the forward and backward integration limits of the phase function.

Our findings could have implications for other subdiffuse techniques. For other techniques, such as spatially resolved reflectance, it will be interesting to see whether $R_{p(\theta b,\theta f)}$ can also be used to model the semiballistic contribution to the reflectance and, if so, which integration angles will have to be used.

Furthermore, since semiballistic photons have shorter pathlengths than diffuse photons, time-resolved or low coherence-based subdiffuse techniques could be developed to separate the semiballistic and diffuse contributions to the reflectance.³³ Ultimately, this would lead to a understanding of the pathlength distribution ($l=ct$), and possible application of Eq. (1) to derive models for SFR spectroscopy.

5.

Conclusion

We developed a model for the reflectance in the absence of absorption ($R_0$) measured with SFR spectroscopy that provides accurate results for a wide range of tissue phase functions. Since the phase function of a specific tissue under investigation is generally not known, a model that is valid for the wide range of phase functions that can be encountered in tissue is essential. We modeled the reflectance as the sum of a diffuse and a semiballistic component. We used the model of Farrell et al.¹⁷ for the diffuse component, solved for an overlapping source and detector fiber. We proposed the parameter $p_{\mathrm{sb}}$ to model the influence of the phase function on the semiballistic contribution to the reflectance. For SFR measurements, $p_{\mathrm{sb}}$ outperforms the subdiffuse parameters $\gamma$, $\delta$, $\sigma$, and $R_{p\mathrm{NA}}$. The new model predicts the measured reflectance substantially better, compared to the commonly used model from Kanick et al.⁷ To determine the tissue optical properties accurately with SFR spectroscopy, the next step will be to include absorption in the model.