Modeling subdiffusive light scattering by incorporating the tissue phase function and detector numerical aperture

To detect small-scale changes in tissue with optical techniques, small sampling volumes and, therefore, short source–detector separations are required. In this case, reflectance measurements are not adequately described by the diffusion approximation. Previous studies related subdiffusive reflectance to γ or σ, which parameterize the phase function. Recently, it was demonstrated that σ predicts subdiffusive reflectance better than γ, and that σ becomes less predictive for lower numerical apertures (NAs). We derive and evaluate the parameter RpNA, which incorporates the NA of the detector and the integral of the phase function over the NA in the backward and forward directions. Monte Carlo simulations are performed for overlapping source/detector geometries for a range of phase functions, reduced scattering coefficients, NAs, and source/detector diameters. RpNA improves prediction of the measured reflectance compared to γ and σ. It is, therefore, expected that RpNA will improve derivation of optical properties from subdiffusive measurements. © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JBO.22.5.050501]

Modeling subdiffusive light scattering by incorporating the tissue phase function and detector numerical aperture 1 Introduction Light that has traveled through tissue carries information about tissue properties, such as structure and biochemical composition. The distance that light travels determines the scale on which the tissue is investigated, as optical properties are averaged over the sampling volume. To detect localized, small-scale changes, techniques with small sampling volumes are, therefore, appropriate. Consequently, short source-detector separations are required to ensure most detected photons will have scattered only once or a few times. In this letter, we will investigate the reflectance measured in this so-called subdiffusive regime, specifically for overlapping source-detector geometries, relevant to, e.g., single fiber spectroscopy measurements.
In the diffusive regime, photon direction is randomized and reflectance is independent of the exact shape of the phase function [pðθÞ, the probability distribution of scattering angles]. In contrast, in the subdiffusive regime, the photon direction is not fully randomized; therefore, measurements are sensitive to the phase function. 1,2 Advances have been made to derive analytic models for subdiffusive reflectance. [3][4][5] However, for overlapping source-detector geometries, challenges remain regarding incorporation of the tissue phase function. 3,4 To model light transport, solutions to the radiative transport equation (RTE) involve expanding the radiance into a series of i spherical harmonics and the phase function into i Legendre polynomials. The latter are weighted with their moments g i , where g 1 is commonly referred to as the scattering anisotropy. For i ¼ 1, the diffusion approximation to the RTE is obtained, with the similarity relation μ s ð1 − g 1 Þ ¼ μ Ã s ð1 − g Ã 1 Þ, which expresses that tissues with a different scattering coefficient μ s and a different scattering anisotropy g 1 , but equal μ 0 s ¼ μ s ð1 − g 1 Þ yield the same reflectance. However, reflectance at short source-detector separations is inadequately described for i ¼ 1. Improvement is possible by increasing i, giving additional similarity relations. For i ¼ 2, Bevilacqua and Depeursinge 6 suggested an alternative form of these similarity relations using the parameter γ 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 ; 3 2 6 ; 4 1 1 Various studies have investigated describing subdiffusive reflectance using γ in combination with μ 0 s , μ a , and the detector diameter d det . 7-10 However, recent work showed that a range of γ values (for the same μ 0 s , μ a , and d det ) can result in the same reflectance. 11,12 Therefore, Bodenschatz et al. 13 incorporated more similarity relations into the parameter σ, employing all phase function moments 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 9 1 where c depends on both the measurement geometry and the scattering regime (empirically, c ¼ 0.5 was found to give the best results). The authors concluded that σ predicts reflectance better than γ; however, σ becomes less predictive for lower numerical apertures (NAs). We consider subdiffusive reflectance as the sum of a diffusive component and a semiballistic component: the diffusive component depends on μ 0 s d det , and, for the semiballistic component, only photons that have experienced a single backscattering event in combination with an arbitrary number of small-angle forward scattering events are considered. To model this contribution, not all details of the phase function are needed. Rather, we consider the phase function within the NA of the detector. The NA characterizes the range of angles over which a detector can accept incoming photons. We propose a theoretically derived parameter R pNA , which is related to the NA of the *Address all correspondence to: Anouk L. Post, E-mail: a.l.post@amc.uva.nl † Contributed equally Journal of Biomedical Optics 050501-1 May 2017 • Vol. 22 (5) detector and the integral of the phase function over the corresponding angular interval in the backward (p NA;backward ) and forward (p NA;forward ) directions, 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 3 ; 6 3 ; 7 1 9 R pNA ¼ We show the theoretical derivation of R pNA and the improved prediction of subdiffusive reflectance based on Monte Carlo simulations by R pNA compared to σ and γ. R pNA , therefore, has the potential to improve models used in subdiffusive measurements to extract optical properties.

Theoretical Derivation of R pNA
We consider semiballistic photons that are backscattered once in combination with an arbitrary number of forward directed scattering events and assume that these photons will only be detected if subsequent scattering events occur at angles smaller or equal to the acceptance angle θ NA (Fig. 1), where θ NA ¼ arcsinðNA∕n sample Þ and n sample is the refractive index of the sample. We introduce the probabilities p NA;backward and p NA;forward 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 4 ; 6 3 ; 5 2 2 p NA;backward ¼ 2π 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 ; 4 6 6 p NA;forward ¼ 2π The contribution of semiballistic photons to subdiffusive reflectance scales with the probability of a single backscattering event in the pathlength interval dl; this probability equals μ s · p NA;backward · dl. To determine the effect of forward semiballistic scattering, we use the result obtained by Wang and Wilson, 14 which gives the probability of a photon traveling a pathlength l while experiencing N forward scattering events 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 ; 3 2 6 ; 7 5 2 prob NA;forward ðl; NÞ ¼ ðp NA;forward · μ s · lÞ N N! e −μ s l : Summing Eq. (6) over all possible numbers of forward scattering events, N from 0 to ∞, we obtain 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 ; 3 2 6 ; 6 9 7 prob NA;forward ðlÞ ¼ e p NA;forward ·μ s l · e −μ s l ¼ e −μ s lð1−p NA;forward Þ Thus, the combined probability of one backscatter event and an arbitrary number of forward scattering events, each occurring within θ NA , is given by integrating over all possible pathlengths, 0 to ∞ 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 ; 6 0 7

Methods
For our Monte Carlo simulations, we modified the software of Prahl et al. 15 (which was the core programming later used in Monte Carlo model of steady-state light transport in multilayered tissues (MCML) 16 software), to allow the use of arbitrary phase functions using the method of Zijp and ten Bosch. 17 The modified code was benchmarked against standard MCML using the Henyey-Greenstein phase function. We modeled geometries where the source and detector geometry were equal and overlapping, where photons were launched from a location based on a uniform distribution across the source with an angle from a uniform angular distribution within the NA. We performed simulations using 15 modified Henyey-Greenstein (mHG), 6 144 double Henyey-Greenstein (dbHG), 18 8 modified power of cosines (MPC), 6 and 46 Reynolds McCormick 19 phase functions (RMC, which is equivalent to the Gegenbauer kernel phase function) employing the parameters specified in Table 1 Fig. 1 Semiballistic scattering: a combination of any number of forward scattering events N in combination with a single backscatter event. We assume photons will only be detected if all scattering angles θ are smaller than or equal to the acceptance angle θ NA of the detector. The scattering angle θ is determined with respect to the previous scattering direction as indicated by the dashed lines. and applying the restrictions g 1 ≥ 0.5 and g 2 < 0.9 to exclude biologically unreasonable phase functions. For each set of phase functions, simulations were performed for three values of μ 0 s d det ∶0.1 (μ 0 s ¼ 10 cm −1 , d det ¼ 100 μm), 1 (μ 0 s ¼ 100 cm −1 , d det ¼ 100 μm), and 9 (μ 0 s ¼ 100 cm −1 , d det ¼ 900 μm) in combination with an NA of 0.22 or 0.5. An absorption coefficient of 0.1 cm −1 and refractive indices inside and outside the sample of 1.35 and 1, respectively, were used. Figure 2 shows the simulated reflectance versus R pNA , σ, and γ for an NA of 0.22 and 0.5, and three different values of μ 0 s d det (0.1, 1, and 9). Reflectance correlates with all three parameters for μ 0 s d det of 0.1 and 1. This correlation is far less pronounced for μ 0 s d det ¼ 9. To compare R pNA , σ, and γ, we determined the spread in each of these three parameters for a chosen reflectance (AE10%) relative to the total range of the parameter ( Table 2). For example, for μ 0 s d det ¼ 0.1, NA = 0.22 and a reflectance of 0.001 (AE10%), γ ranges from 1.43 to 1.66. The total range of γ is 0.68 to 2.43, so the variability is calculated as 0.23/1.75 = 0.13. For all combinations of NA (0.22 and 0.5) and μ 0 s d det values (0.1 and 1), we determined these variability values for three reflectance values and found the variability to be lowest for R pNA compared to σ and γ.

Discussion
We have theoretically derived the parameter R pNA to model subdiffusive light scattering for overlapping source-detector geometries, and we have shown with Monte Carlo simulations that the reflectance depends on R pNA . In comparison to σ and γ, R pNA improves prediction of the reflectance. Our findings indicate that the reflectance does not depend on the details of the phase function-knowledge of all Legendre moments is not required-but only on the magnitude of the phase function within the acceptance angle of the detector (in the backward and forward directions).  Currently, the parameter γ is widely used to model subdiffusive reflectance. Since R pNA improves prediction of the measured reflectance, incorporating R pNA into subdiffusive models is expected to improve their reliability and thereby the estimation of other optical properties, such as μ 0 s . This improvement is significant because subdiffusive measurements have the potential to detect small-scale tissue changes.
For higher values of μ 0 s d det , the reflectance becomes more diffusive and, therefore, depends less on R pNA , σ, and γ. The contribution of diffuse photons is a remaining challenge for describing subdiffusive scattering. In our simulations, we kept the absorption coefficient constant. For subdiffusive reflectance, path lengths are short; therefore, the effect of μ a will be minor. While the absorption affects the diffuse contribution to the reflectance, it has a minor effect on semiballistic photons. To implement R pNA for subdiffusive measurements, models have to be developed incorporating R pNA , the diffusive reflectance, and absorption for specific measurement geometries, such as single fiber reflectance spectroscopy.

Conclusion
In conclusion, the theoretically derived parameter R pNA predicts subdiffusive light scattering for overlapping source-detector geometries. Consequently, the reflectance does not depend on the details of the entire phase function, but on the phase function within the acceptance angles of the detector. Since R pNA improves prediction of the measured reflectance compared to σ and γ, the use of R pNA is expected to improve derivation of optical properties from subdiffusive measurements.

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