Subdiffuse scattering model for single fiber reflectance spectroscopy

Abstract. To detect small-scale changes in tissue with optical techniques, small sampling volumes are required. Single fiber reflectance (SFR) spectroscopy has a sampling depth of a few hundred micrometers. SFR spectroscopy uses a single fiber to emit and collect light. The only available model to determine optical properties with SFR spectroscopy was derived for tissues with modified Henyey–Greenstein phase functions. Previously, we demonstrated that this model is inadequate for other tissue phase functions. We develop a model to relate SFR measurements to scattering properties for a range of phase functions, in the absence of absorption. Since the source and detector overlap, the reflectance cannot be accurately described by diffusion theory alone: SFR measurements are subdiffuse. Therefore, we describe the reflectance as a combination of a diffuse and a semiballistic component. We use the model of Farrell et al. for the diffuse component, solved for an overlapping source and detector fiber. For the semiballistic component, we derive a new parameter, psb, which incorporates the integrals of the phase function over 1 deg in the backward direction and 23 deg in the forward direction. Our model predicts the reflectance with a median error of 2.1%, compared to 9.0% for the currently available model.


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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 1 7 2 R½μ s ; pðθÞ; μ a ¼ Z ∞ 0 p½l; μ s ; pðθÞe −μ a l dl; (1) where μ s is the scattering coefficient, pðθÞ is the phase function (the probability distribution of scattering angles), and μ a is the absorption coefficient. Equation (1) describes the reflectance as the sum of photons with different pathlengths, according to the pathlength distribution pðlÞ. The pathlength distribution only depends on the scattering properties, μ s and pðθÞ. Absorption is accounted for by weighting each pathlength contribution according to the Beer-Lambert law (the exponential term in the integral). Equation (1) enables separate modeling of the effects of scattering and absorption on reflectance.
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 ) 7 and the second part includes absorption. 8 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 μ a and an effective pathlength hLi: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 2 3 2 R ¼ R 0 e −μ a hLi : (2) The model for R 0 describes the reflectance as a function of the product of the reduced scattering coefficient and the fiber diameter (μ 0 s d), the phase function parameter γ (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 sample ): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 1 3 3 · ð1 þ 0.632γ 2 e −2.308γ 2 ·μ 0 s d Þ · ðμ 0 s dÞ 0.574γ 2.308γ 2 þ ðμ 0 s dÞ 0.574γ : *Address all correspondence to Anouk L. Post, E-mail: a.l.post@amsterdamumc .nl 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. 9 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, 10 brain, 12 breast, 13 and muscle, 14 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∕μ 0 s . 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. 17 for spatially resolved diffuse reflectance, as solved by Faber et al. for an overlapping source and detector fiber. 18 The diffuse component of the reflectance depends only on the product μ 0 s d, where μ 0 s ¼ μ s ð1 − g 1 Þ, μ 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 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 (γ, 21 δ, 22 σ, 23 and R pNA 9 ). Therefore, we model the semiballistic contribution to the reflectance as a function of μ 0 s d and p 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. 7 for a range of phase functions, μ 0 s d values and NAs.

Monte Carlo Simulations
We performed absorption-free MC simulations to investigate the relationship between R 0 and tissue-scattering properties for SFR spectroscopy. We modified 9 the software of Prahl et al. 24 and Wang et al. 25 to enable simulation of an overlapping source/detector geometry and to allow the use of any phase function using the method of Zijp and ten Bosch. 26 Photons were launched from a location based on a uniform distribution across the source with an angle from a uniform angular distribution within the launch/acceptance angle of the fiber θ acc , where θ acc ¼ arcsinðNA∕n sample Þ. Photons were detected if they reached the fiber face at an angle within θ acc . We performed simulations for NAs of 0.10, 0.22, and 0.50. The refractive index was set to 1.35 for the tissue, 1.45 for the fiber face, and 1.00 for the medium above the sample. For all the MC simulations, we ran each simulation three times and ensured that enough photons had been launched such that the standard deviation over the mean of the reflectance for each set of three simulations was less than 2%. In conventional MC simulations that include absorption, the photon weight is reduced at the end of every step. There, to prevent endless tracking of photons with small weights, usually a "roulette" routine is implemented to terminate these photons when their weight drops below a certain value. In our absorption-free MC simulations, weights are maintained and another process for termination is needed. We terminated photons at a set distance from the fiber face (termination distance) and we checked that at least 99.9% of detected photons had traveled up to a distance from the fiber face less than 75% of the termination distance. We performed simulations using MHG, TTHG, and RMC phase functions, employing the parameters specified in Table 1 and applying the restrictions g 1 ≥ 0.5 and g 2 < 0.9 to exclude biologically unreasonable phase functions. This resulted in 207 phase functions (15 MHG, 146 TTHG, and 46 RMC) with fairly equally distributed g 1 values between , and 9 (μ 0 s ¼ 100 cm −1 , d ¼ 900 μm). The simulations with μ 0 s d ¼ 0.1 were used to derive a new parameter to describe the semiballistic contribution to the reflectance p sb . To confirm that, for μ 0 s d ¼ 0.1, the detected photons are indeed semiballistic, we determined the fraction of detected photons that underwent a single backscatter event in combination with an arbitrary number of forward-scattering events. The simulations with μ 0 s d ¼ 0.1 were used to investigate the relationship between R 0 and p sb , γ, δ, σ, and R pNA .
In addition, we performed simulations with 10 phase functions (selected such that they yielded reflectance values between 10 −5 and 4·10 −3 in 10 equal steps on a logarithmic scale, for an NA of 0.22 and μ 0 s d of 0.1) for 20 values of μ 0 s d between 0.1 and 100 (μ 0 s ¼ 10 to 10000 cm −1 , equally spaced in 20 steps on a logarithmic scale, d ¼ 100 μm). All the simulations combined were used to derive constants in our scattering model, as well as to test the accuracy of this new model.
We also investigated the limits of the new model for low values of μ 0 s d. Since part of our model is based on diffusion theory, it can be expected that there is a limit where μ 0 s d is too low for diffusion theory to accurately describe part of the reflectance. Therefore, we performed MC simulations for the same 10 phase functions, with μ 0 s d values of 0.005, and 0.01 to 0.09 in steps of 0.01 (d ¼ 100 μm), and we compared the simulated reflectance to the modeled reflectance.

Scattering Model
We model the reflectance (R 0 ) as the sum of a diffuse (R SFR;dif ) reflectance and a semiballistic (R SFR;sb ) reflectance: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 3 9 2 R 0 ¼ R SFR;dif þ R SFR;sb: (4) 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 μ 0 s d, almost all detected photons will be diffuse, whereas for low values of μ 0 s d, most detected photons will be semiballistic. Therefore, we rewrite Eq. (4) using this ratio E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 2 9 5 X ¼ as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 6 3 ;

Diffuse component
We describe the diffuse component, R dif , using the model of Farrell et al. 17 for spatially resolved reflectance as solved by Faber et al. 18 for an overlapping source and detection fiber in the absence of absorption: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; Here, R μa¼0 ðρ; μ 0 s Þ is the diffuse reflectance versus radial distance for a pencil beam illumination, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 3 2 6 ; 7 5 2 and pðρ; dÞ is the distribution function of distances over the fiber face, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 6 5 8 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. 27 The parameter A depends on the refractive index mismatch between the tissue and the fiber. 28 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 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 dif by the collection efficiency of the fiber (η c ): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 4 7 9 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 η c;L is equal to 29 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 3 2 6 ; 4 0 4 η c;L ¼ sin θ 2 acc : 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, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 3 1 6 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 μ 0 s d (1000; μ 0 s ¼ 10 5 cm −1 , d ¼ 100 μm) and compared these to the collection efficiency η c . For high values of μ 0 s d, the reflectance will be diffuse and, thus, R SFR;sb will approach 0 and R dif will approach 1. Therefore, the reflectance should equal to the collection efficiency for high values of μ 0 s d.

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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 1 2 3 pðcos θÞ ¼ 1 4π X n ð2n þ 1Þg n P n ðcos θÞ; where p is the probability of scattering at an angle θ, P n are the Legendre polynomials of order N, and g n are the Legendre Journal of Biomedical Optics 015001-3 January 2020 • Vol. 25 (1) moments. The scattering model from Kanick et al. 7 describes the influence of the phase function by incorporating the parameter Next to γ, three other parameters have been proposed that could be used to incorporate the phase function influence in models to describe subdiffuse reflectance: δ , 22 σ, 23 and R pNA . 9 The parameters δ and σ also incorporate moments of the phase function: We recently introduced the parameter R pNA , 9 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 scattering angles ≤ θ acc . The full derivation is described in our previous paper. 9 In short, based on this assumption, R pNA is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 5 5 0 R pNA ¼ Here, p NA;b is the probability of a scattering event to occur within θ acc in the backward direction, which equals the integral of the phase function over the acceptance angle in the backward direction: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 6 3 ; 4 6 8 p NA;b ¼ 2π Z π π−θ acc pðθÞ sinðθÞdθ: (15) Here p NA;f is the probability of a scattering event to occur within θ acc in the forward direction: We showed previously that for SFR spectroscopy the parameter R pNA predicts the reflectance better than γ or σ, 9 but we did not yet validate the assumptions behind the derivation of R pNA .
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 μ 0 s 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 pNA , which we named R pðθb;θfÞ : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 1 9 4 R pðθb;θfÞ ¼ E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 9 1p f ðθ f Þ ¼ 2π We investigated the relation between R 0 and R pðθb;θ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ðθb;θfÞ for a chosen reflectance, relative to the total range of log 10½R pðθb;θfÞ (Fig. 1). For each NA, we determined the relative dispersion for five different reflectance values (AE10%) equally spaced on a logarithmic scale. This relative dispersion is a measure for the suitability of R pðθb;θfÞ to model the reflectance, because the dispersion in the parameter will influence the uncertainty in all the estimated optical properties when R pðθb;θ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.
Based on these results, we defined the parameter p sb as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 3 2 6 ; 3 0 1 We compared p sb to γ, σ, δ, and R pNA by determining the relative dispersion for three different reflectance values (AE10%) per NA. Since we developed p sb for the semiballistic contribution to the reflectance, we compared these parameters for μ 0 s d ¼ 0.1, where most detected photons are semiballistic.

Ratio diffuse and semiballistic photons
We use p sb in the model for the ratio between semiballistic and diffuse photons X. For low values of μ 0 s d, the diffuse reflectance scales with ðμ 0 s dÞ 2 . Therefore, we include that factor in the denominator of Eq. (21). The semiballistic reflectance will scale differently with μ 0 s d. We assume R SFR;sb scales with p sb ðμ 0 s dÞ −B . This difference in scaling is accounted for through the parameter a 3 . Fig. 1 Concept of judging the suitability of the parameter R pðθb;θf Þ for a specific forward and backward integration angle. For a certain reflectance value (R 0 ), the dispersion of R pðθb;θf Þ for that reflectance (light blue) is divided by the entire range of R pðθb;θf Þ (green). This relative dispersion was determined for five reflectance values (depicted in red) for μ 0 s d ¼ 0.1, per set of 207 simulations with the same NA.

Monte Carlo Simulations Single Fiber
Reflectance Spectra To demonstrate the use of our model to extract optical properties from measured spectra, we performed MC simulations for Intralipid 20% diluted to a volume fraction of 0.05, using the scattering coefficient and phase function from literature 30 and no absorption. To extract optical properties from SFR measurements, measurements with two fibers of different diameters are generally used (also referred to as multidiameter SFR or MDSFR). For robust fit results, the number of fit parameters should be considerably smaller than the number of data points.
For an SFR measurement, the number of data points (reflectance values) is equal to the number of wavelengths, λ n . Separate values for μ 0 s and p sb per wavelength cannot be determined directly using a fit on an SFR spectrum, since the number of fit parameters would equal 2λ n . Therefore, the number of fit parameters is reduced by modeling the reduced scattering coefficient as where a is the scattering amplitude, b is the scattering slope, and λ 0 is a reference wavelength. The wavelength dependence of tissue phase functions, in general, and p sb , specifically, is not well characterized. Therefore, a value of p sb will have to be fitted for each wavelength within the spectrum. For a single measured spectrum, the number of data points will equal λ n , but the number of fit parameters will equal 2 þ λ n (2 for μ 0 s and λ n for p sb ). Such an underdetermined system will not provide robust fit results. To overcome this issue, measurements can be performed using two different fiber diameters. In that case, the number of data points will equal 2λ n and the number of fit parameters will equal 2 þ λ n . Therefore, we performed MC simulations of spectra for two different fiber diameters: 300 and 600 μm. We modeled spectra from 400 to 900 nm in steps of 5 nm. For the fit of μ 0 s , we used a reference wavelength of 600 nm.

Phase Function Parameter Semiballistic Contribution (p sb )
For μ 0 s 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 pNA that semiballistic photons are only detected if all the scattering events occurred at scattering angles ≤ θ 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 θ 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 θ 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 θ acc and undergoes only a single scattering event can be detected if that scattering angle is 2 · θ acc .
Since detected photons also underwent scattering events at angles larger than θ 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ðθb;θfÞ [Eqs. (17)- (19)]. The optimization of the integration angles was based on minimizing the relative dispersion in the parameter log 10½R pðθb;θ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.
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 sb [Eq. (20)] in our model for the semiballistic contribution to the reflectance. Figure 4 shows the relation between the reflectance and γ, σ, δ, R pNA , and p sb for simulations with an NA of 0.22 and Overall, p 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 sb outperforms the other parameters for all three reflectance values. For an NA of 0.22, p 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 pNA and p sb are similar (and both outperform the other parameters), which is expected since the integration angles for R pNA are 22 deg for an NA of 0.50. 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).

Scattering Model
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 μ 0 s d (1000; μ 0 s ¼ 10 5 cm −1 , d ¼ 0.01 cm) to the collection efficiency using Eq. (12) in combination with the a 1 values from Table 4. For high values of μ 0 s , 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 SFR;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 μ 0 s d, the total reflectance equals the diffuse reflectance (dashed black line), but for lower values of μ 0 s d there is an additional semiballistic contribution to the reflectance. With increasing p sb , the fraction of semiballistic photons increases.

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 μ 0 s d for our model, by determining the difference for all the simulations with 10 different phase functions (0.005 < μ 0 s d < 10 3 ). For values of μ 0 s d that are lower than 0.1, the difference increases significantly (Fig. 6).
Next, we determined the difference between the reflectance predicted by the model and the reflectance from the MC Table 3 Variability of p sb , R pNA , σ, δ, and γ for μ 0 s d ¼ 0.1, defined as the dispersion in p sb , R pNA , σ, and γ values for a chosen reflectance (+/-10%) relative to the total range of each parameter.  Table 5 Resulting fit parameters per NA using Eqs. (5)-(10), (12), (20), and (21) and a 1 ¼ 1.11 on simulated SFR measurements. The 95% confidence intervals on these fit parameters are indicated.  simulations for all simulations where μ 0 s d ≥ 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. 7 [ Fig. 7(b)], where the median error is 9.0%, with a standard deviation of 36.8% and a maximum error of 303%.

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-μm fiber and fluctuates around 0% for the 600-μm fiber. The obtained μ 0 s and p sb values are close to the values used in the simulations.

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. 7 (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 sb , a measurement with two fibers of different diameters is required. The fitted values we obtained for μ 0 s are very close to the simulated μ 0 s . The difference between the fitted and the simulated values of p sb are slightly larger. A model for the wavelengthdependent behavior of p 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 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ðθb;θ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 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 γ, δ, σ, and R pNA . The parameters γ, δ, and σ 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. 17 for an overlapping source and detection fiber. 18 For low values of μ 0 s d, the diffuse reflectance scales with ðμ 0 s 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 μ 0 s d, which is accounted for in the parameter a 3 .
As expected, the fraction of diffuse photons increases for higher values of μ 0 s d and lower values of p 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 ∼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 μ 0 s 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. 31 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 μ 0 s d is too low for diffusion theory to accurately describe the diffuse contribution to the reflectance. For values of μ 0 s d below 0.1, our model becomes inaccurate. For tissues with low values of μ 0 s , it is, therefore, advisable to use fibers with a larger diameter. For example, with a 500-μm fiber, our model accurately predicts the reflectance for μ 0 s ≥ 2 cm −1 , and for most tissues μ 0 s will be larger than 2 cm −1 . 32 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ðθb;θfÞ can also be used to model the semiballistic contribution to the reflectance and, if so, which integration angles will have to be used. Journal of Biomedical Optics 015001-9 January 2020 • Vol. 25 (1) Furthermore, since semiballistic photons have shorter pathlengths than diffuse photons, time-resolved or low coherencebased subdiffuse techniques could be developed to separate the semiballistic and diffuse contributions to the reflectance. 33 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.

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. 17 for the diffuse component, solved for an overlapping source and detector fiber. We proposed the parameter p sb to model the influence of the phase function on the semiballistic contribution to the reflectance. For SFR measurements, p sb outperforms the subdiffuse parameters γ, δ, σ, and R pNA . The new model predicts the measured reflectance substantially better, compared to the commonly used model from Kanick et al. 7 To determine the tissue optical properties accurately with SFR spectroscopy, the next step will be to include absorption in the model.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.