A predictive model of backscattering at subdiffusion length scales

Vladimir Turzhitsky; Andrew Radosevich; Jeremy D. Rogers; Allen Taflove; Vadim Backman

doi:10.1364/BOE.1.001034

1. Introduction

Diffusion approximations are often utilized to allow fast predictions of reflectance signals. These approximations involve a simplification of the transport equation and are typically validated with a more exact numerical method such as Monte Carlo [1]. Generally, diffusion approximations can be accurate for predicting the reflectance when the observed lateral distance (r) is much greater than the transport mean free path (ls*). Diffusion approximations are not accurate at small distances of light transport (e.g. source-detector separation), r, because they do not take the shape of the phase function into account. While at diffusion length scales (r/ls*>>1), light transport is primarily governed by the value of the transport mean free path (ls*), at sub-diffusion distances (r/ls*<1), the shape of the phase function may significantly affect the radial reflectance distribution. Foster and others have shown that the accuracy can be improved by accounting for higher order moments of the phase function in the P3 approximation [2,3]. However, even the P3 approximation becomes inaccurate for reflectance closer than ½ ls* [2].

For certain applications, it is important to be able to model and predict the backscattering signal from turbid media at source-detector separation distances that are significantly smaller than the transport mean free path ls*. For example, cancer detection often requires the isolation of a signal from superficial tissue such as the epithelium or mucosa. In many tissues, the thickness of the epithelium is much smaller than ls* therefore requiring a source detector separation much smaller than ls*. For this reason, several groups have developed fiber probes that sample small source-detector separations [4–6]. Other methods for sampling small radial transport distances include polarization gating or coherence based methods such as Low-coherence Enhanced Backscattering [7,8]. Thus far, these methods have relied on time intensive computational solutions of the radiative transport equation (RTE), typically with Monte Carlo simulations that predict the backscattering signal as a function of optical properties [9,10]. These computational solutions have established that variations in the anisotropy factor and the shape of the phase function result in substantial variations to the reflectance at small r/ls* [2,11]. A fast and accurate predictive model of scattering at these small length scales is therefore of great interest for measuring properties of the scattering phase function as well as predicting epithelial tissue scattering.

In this paper, we introduce a simple approach that allows for accurate prediction of the backscattering signal down to length scales that are several orders of magnitude smaller than ls*. The approach involves the construction of a simple model that predicts an infinitely narrow normally incident illumination beam response, termed p(r), from a turbid scattering medium. The response p(r) is a fundamental property of the turbid medium and is the objective for predictive modeling of most diffusion approximation models and Monte Carlo methods. If p(r) is known, both the effects of a finite source and a numerical aperture can be modeled. We consider two types of phase functions in order to construct the models: the commonly used Henyey-Greenstein phase function, and a more general two parameter phase function which encompasses the Henyey-Greenstein phase function and uses parameters that quantify the sample refractive index correlation function. As an experimental example, we measure the reflectance distribution from a tissue phantom composed of a mixture of polystyrene microspheres using Low-coherence Enhanced Backscattering and compare the measured distribution at small length scales to the newly developed model.

1. Monte Carlo Simulation of Reflectance

One of the main determinants of the accuracy of a Monte Carlo simulation at small length scales is the choice of the phase function. Although the Henyey-Greenstein phase function has been described as a sufficiently accurate choice for prediction of backscattering at intermediate length scales (r~ls*), accurate modeling of small length scales (r<<ls*) requires a more general choice of phase function. For this purpose, we will follow a recently developed model that is based on the Whittle-Matérn correlation function [12,13].

The model implements the Born approximation in order to obtain the phase function from the refractive index correlation function. The Born (i.e. weakly scattering) approximation is valid in the regime relevant for soft biological tissue [14]. In the Born approximation, the differential cross section and, thus, the phase function, are completely defined through the Fourier transform of the refractive index correlation function [15]. In turn, the refractive index is a linear function of the local density of tissue macromolecules and the refractive index correlation function is proportional to that of the mass density. Thus, one can calculate the phase function if the mass density correlation is known (and vice versa). There have been several hypotheses on the best functional form that can model the refractive index correlation function in tissue, but one convenient expression can encompass nearly all of these possibilities. The refractive index fluctuations of biological tissue can be modeled with the Whittle-Matérn correlation function [16]: $B_{n} (r) \propto Δ n^{2} {(r / l_{c})}^{m - 3 / 2} K_{m - 3 / 2} (r / l_{c})$ , where Δn² is the variance of the refractive index fluctuations, l_c is the correlation length, and m is a parameter that determines the form of the function. The function K_m-3/2 denotes the modified Bessel function of the second kind of order m-3/2. When m < 1.5, B_n(r) is a power law, thus corresponding to a mass fractal medium with mass fractal dimension D_mf = 2m. 1.5< m <2 corresponds to a stretched exponential function, m = 2 corresponds to an exponential function, and as m becomes much larger than 2, B_n(r) approaches a Gaussian function. The correlation length l_c has different physical meaning depending on the type of the correlation function. For m = 2, $B_{n} (r) \propto Δ n^{2} \exp (- r / l_{c})$ while in a mass fractal case of m < 1.5, l_c represents the upper length scale at which the correlation function loses its fractal behavior. The differential scattering cross section can be derived by applying the Born approximation to the Whittle-Matérn correlation function [13]. For the scalar wave case, the phase function assumes the following form:

\begin{array}{l} F (\cos θ) = \frac{2 \hat{g} (m - 1) \cdot {[1 - 2 \hat{g} \cos θ + {\hat{g}}^{2}]}^{- m}}{{(1 - \hat{g})}^{2 - 2 m} - {(1 + \hat{g})}^{2 - 2 m}} \\ \hat{g} = 1 - \frac{\sqrt{1 + 4 {(k l_{c})}^{2}} - 1}{2 {(k l_{c})}^{2}} \\ k l_{c} = \frac{\sqrt{\hat{g}}}{1 - \hat{g}}, \end{array}

where k = 2π/λ and the normalization is such that

\int F (\cos θ) \cdot d \cos θ = 1

. We will refer to this phase function as the Whittle-Matérn phase function. The phase function can also be expressed as a function of m and klc without any change to the normalization:

F (θ) = \frac{2 {(k l_{c})}^{2} (m - 1) \cdot {[1 + {(2 k l_{c} \sin (θ / 2))}^{2}]}^{- m}}{1 - {(1 + {(2 k l_{c})}^{2})}^{1 - m}} .

It is important to note that a special case is observed when the value of m = 1.5. In this special case the correlation function is that of the space filling random field and the phase function simplifies to the commonly used form known as the Henyey-Greenstein phase function. The parameter ĝ then becomes the average cos(θ), also known as the anisotropy factor g. For other values of m, g is given by taking the forward moment:

g = 〈 \cos θ 〉 = {\begin{cases} \frac{{(1 - \hat{g})}^{2 - 2 m} \cdot (1 + {\hat{g}}^{2} - 2 \hat{g} (m - 1)) - {(1 + \hat{g})}^{2 - 2 m} \cdot (1 + {\hat{g}}^{2} + 2 \hat{g} (m - 1))}{2 \hat{g} (m - 2) [{(1 + \hat{g})}^{2 - 2 m} - {(1 - \hat{g})}^{2 - 2 m}]} m \neq 2 \\ \frac{1 + {\hat{g}}^{2}}{2 \hat{g}} + \frac{\ln (1 + \hat{g}) - \ln (1 - \hat{g})}{g [{(1 - \hat{g})}^{- 2} - {(1 + \hat{g})}^{- 2}]} m = 2 \end{cases}

Fig. 1(a) shows examples of the Whittle-Matérn phase functions with the same value of m, and varying values of g, while Fig. 1(b) shows examples of phase functions with the same value of g and varying values of m. The parameter g influences the width of the phase function while m influences the shape of the phase function independently of the width. There are two cases in the generalized phase function which are removable discontinuities: m = 1 and g = 0. We can evaluate the phase function for these cases by employing L’Hospital’s rule:

Fig. 1 Example phase functions. (a) The Henyey-Greenstein case (m = 1.5) for various values of g. (b) The generalized Whittle-Matérn phase function for various values of m (g = 0.9).

Download Full Size | PDF

F (\cos θ) = {\begin{cases} \frac{\hat{g} \cdot {[1 - 2 \hat{g} \cos θ + {\hat{g}}^{2}]}^{- 1}}{\ln (1 + \hat{g}) - \ln (1 - \hat{g})} m = 1 \\ 1 / 2 \hat{g} = 0 \end{cases}

In these two cases, g becomes:

g = {\begin{cases} \frac{1 + \hat{g}}{2 \hat{g} [\ln (1 + \hat{g}) - \ln (1 - \hat{g})]} m = 1, \hat{g} \neq 0 \\ 0 \hat{g} = 0 \end{cases} .

The reflectance distribution was calculated with the Whittle-Matérn phase function by implementing the Monte Carlo method. An existing Monte Carlo code was modified to implement the generalized phase function [9,17]. The code was validated by comparing the results for m = 1.5 (Henyey-Greenstein case) with existing codes that implement the Henyey-Greenstein phase function [10]. Simulations with values of g varying from 0 to 0.98 and values of m varying from 1.01 to 1.9 were obtained for the backscattering direction (0-10°). We found that the variations of backscattering probability distributions were small within the 10° angular collection range when the backscattering probability distribution was stored as a function of the position of the final scattering event. Therefore, all reflectance distributions were stored as a function of the position of the last scattering event. ls* was maintained at 100μm with a scattering slab thickness of 1cm, resulting in a scattering medium that approaches semi-infinite, with less than 2% of the intensity transmitting through the entire thickness of the slab. The infinitely narrow illumination beam was oriented orthogonally to the scattering medium. The boundary at the interface of the scattering medium was assumed to be index-matched and absorption was not present. The scattering angle, θ, in the Monte Carlo simulation was chosen by expressing the probability of a selected angle as a function of the random variable ξ:

\cos θ = \frac{1 + {\hat{g}}^{2} - {[ξ ({(1 - g)}^{2 - 2 m} - {(1 + g)}^{2 - 2 m}) + {(1 + g)}^{2 - 2 m}]}^{\frac{1}{1 - m}}}{2 g}

where ξ is uniformly distributed between 0 and 1 [10]. The azimuthal angle was chosen from a uniform random distribution: ψ = 2πξ. The remaining Monte Carlo simulation elements were identical to previously developed methodology for light propagation in turbid media [10].

2. Model of Reflectance: Henyey-Greenstein Phase Function

The backscattering distributions were obtained from Monte Carlo simulations by collecting rays within 10° of the backward direction and normalizing the reflectance such that $\int P (s) d s = 1$ , where s = r/ls* and $P (r) = \int_{0}^{2 π} p (r, θ) r d θ = 2 π \cdot r \cdot p (r)$ , r and θ being the polar coordinates in a plane perpendicular to the illumination beam. In a Monte Carlo simulation, P(r) is the obtained reflectance distribution that is collected with azimuthally integrated radial storage (θ is the azimuth angle). Figure 2(a) shows example P(r/ls*) curves obtained from Monte Carlo simulations using the Henyey-Greenstein phase function for four different values of g. All of the length scales in the Monte Carlo simulation are determined by ls, the mean free path. Additionally, it is known that the determining length scale in the diffusion regime is ls*. Therefore, the axes in Fig. 2 are normalized with respect to ls* in order to be scalable for any value of ls as well as observe the convergence of the results in the diffusion regime. All of the curves can be translated into units of P(r) by multiplying the abscissa axis by ls* and dividing the ordinate axis by ls*. Another words, P(r) = P(ls*·s)/ls*. The division by ls* is required due to the change of variable in the normalization integral (ds = dr/ls*) such that $\int P (r) d r = 1$ .

Fig. 2 Scaling relationships of P(r/ls*) from the Henyey-Greenstein phase function. (a) P(r/ls*) curves for varying values of g. The curves have variations at r/ls*<1 and converge for larger values of r/ls*. (b) Probability difference obtained by subtracting P(r/ls*) for the isotropic case (g = 0) from P(r/ls*) of a given g. (c) Probability difference curves for varying values of g with the amplitude rescaled by a coefficient that depends only on g.

Download Full Size | PDF

Most diffusion approximations make the simplifying assumptions of isotropic scattering. We can evaluate the effect of anisotropy on P(r/ls*) at subdiffusion length scales by subtracting P(r/ls*) curves for isotropic scattering from P(r/ls*) for non-isotropic cases (i.e. g > 0). In Fig. 2(b), three difference curves are plotted for g values of 0.9, 0.8, and 0.7. Note that the integral of each difference curve is 0 because the integral of P(r/ls*) is always 1. The curves in Fig. 2(b) have very similar shapes but varying amplitudes. When each of these curves is rescaled by a constant that depends on g, they closely overlap [Fig. 2(c)].

We can therefore employ a predictive model of P(r/ls*) that depends on just two simulation results: P(r/ls*) for g = 0 and P(r/ls*) for a particular g>0. While any value of g>0 can be used, we use g = 0.9 in the following analysis for convenience (this results in accurate prediction within the range of tissue anisotropy):

\begin{array}{l} P_{g} = P_{0} + c (g) [P_{0.9} - P_{0}] \\ c (g) = a g^{b}, \end{array}

where c(g) is an empirical model for the coefficients that multiply the difference term. The values of the constants a and b are approximately 1.244 and 2.338 respectively. The shortened notation P_g represents P(r/ls*) for a given value of g (e.g. P_0.9 = P(r/ls*) for g = 0.9). The values of c(g) were determined by fitting Monte Carlo results for a particular g to the expression for P_g in Eq. (7). The values of c(g) and the empirical model for c(g) are plotted in Fig. 3(a) . We can understand the difference between P_0.9 and P₀ as the alteration in the backscattering due to anisotropy. As g increases, the anisotropy contribution increases in amplitude but retains a very similar radial shape. This allows for a predictability of P(r) for any value of g with only two reference P(r) distributions. Fig. 3(b) shows a comparison of the Monte Carlo simulations and the model based on the difference relationship. Fig. 3(c) further illustrates the details of the model fit at small values of r/ls*. Note that the fits for g = 0 and g = 0.9 are not shown because the model and the Monte Carlo result are identical for those two cases (c = 1 when g = 0.9 and c = 0 when g = 0). The model has excellent agreement for values of g that are close to 0.9, but begins to deviate slightly at g = 0.7. As r/ls* becomes large, all of the curves converge and the backscattering can be predicted with an isotropic scattering model of equivalent ls*.

Fig. 3 Backscattering model for P(r/ls*) based on the Henyey-Greenstein phase function utilizing the difference method. (a) Value of amplitude coefficient as a function of g, c(g). The dots represent a least square fit to the Monte Carlo data and the solid line represents the model for c(g) from Eq. (7). The amplitude increases as a function of g follow a power law. (b) Comparison of Monte Carlo simulations (green, red, and purple curves) and model (black curves) for three values of g. g = 0.9 is not shown because the model and simulation are, by definition, identical for that case. (c) is a rescaled version of (b) shown for a smaller range of r/ls*.

Download Full Size | PDF

Another possible approach for predicting the backscattering signal at small radial distances is through an implementation of principle component analysis (PCA). PCA is a variance reduction technique that is often used to obtain dependencies when a large number of variables are present. The analysis typically involves mean-centering the data (i.e. subtracting the mean from each variable), followed by a transformation which decomposes the data into orthogonal components that explain the largest proportions of the variance within the data set. In order to apply this method to build a model of P(r/ls*), we used each value in r/ls* as an input variable. Instead of mean centering, we subtracted the P(r/ls*) curve for g = 0. The effect of subtracting the isotropic P(r/ls*) is similar to that of mean-centering, but results in a more predictable model that is independent of the particular reflectance distributions used in the PCA analysis. We then obtained a series of principle components and found that when the first three components are used, P(r,ls*) can be predicted more accurately than the single-component difference model described above. The first three principle components (PC1-PC3) predicted 99.966%, 0.027%, 0.002% of the variance in the data, respectively. Based on this model, P(r/ls*) can be predicted according to:

P_{g} = P_{g = 0} + c_{1} (g) P C 1 + c_{2} (g) P C 2 + c_{3} (g) P C 3,

where c₁, c₂, and c₃ are the weights of the principle components. We utilized a polynomial equation to fit the weights with the order of the polynomial chosen such that the R² coefficient is greater than 0.99. Fig. 4(a) shows a comparison of the principle component model with Monte Carlo simulation results for three tissue-relevant values of g. The r/ls* axis is in log scale, showing that the model is in excellent agreement with the Monte Carlo simulations for the entire simulated range of 0.001< r/ls*<10. Fig. 4(b) shows the same comparison in linear scale for the subdiffusion range of r/ls* < 1, again, showing excellent agreement. Fig. 4(c) shows the distributions of the three principle components that were used in the model. Note that the contribution of each successive component decreases, with higher components being noisier. Fig. 4(d) is a plot of the weights of the three components along with the polynomial fits that are used for the predictive model.

Fig. 4 Backscattering model for the Henyey-Greenstein case utilizing principle component analysis method. (a) Semi-log plot of Monte Carlo simulations and the PCA model for three values of g. (b) Comparison of PCA model and Monte Carlo simulation for r/ls* < 1. (c) Semi-log plot of the three principle components. Note that higher principle components have less amplitude and contribute less to P(r/ls*) prediction. (d) Plot of coefficients that multiply the principle components and their polynomial fits used to obtain the predictive model.

Download Full Size | PDF

3. Model of Reflectance: Whittle-Matérn Phase Function

To obtain a generalized model for reflectance based on the Whittle-Matérn phase function, we simulated P(r/ls*) distributions for varying values of g and m. In Fig. 5(a) , P(r/ls*) curves for four values of m are shown with a constant anisotropy factor of g = 0.9. The isotropic component is subtracted from these curves in Fig. 5(b). From Fig. 5(b), it is apparent that a simple scaling in amplitude cannot account for the difference between these curves. There is an m-dependent alteration in the shape of the non-diffuse component of the curves. However, for a given value of m, changes in g only alter the amplitude of the non-diffuse component [Fig. 5(c)]. Therefore, it is clear that another component needs to be introduced that can account for the alterations in the shape of P(r/ls*) due to varying m. We can extend the difference model developed for the Henyey-Greenstein phase function (m = 1.5) discussed in the previous section by defining a second difference component that is calculated by subtracting the isotropic probability from P(r/ls*) for g = 0.9 and a particular m. In our analysis, we chose m = 1.01. This value of m was chosen because the shape of P(r/ls*) becomes dramatically altered as m approaches 1. The shapes of the two difference components are compared in Fig. 5(d). P(r/ls*) can then be predicted according to a two-component model:

Fig. 5 P(r/ls*) distributions for the generalized Whittle-Matérn phase function. (a) P(r/ls*) dependence on m for g = 0.9. The shape of the P(r/ls*) gradually changes with varying m. (b) Probability difference between P(r/ls*) curves and the isotropic P(r/ls*) for g = 0.9 and varying m. Unlike changes in g, alterations in m cannot be accounted for by a simple scaling of the amplitude of the probability difference curves. (d) Probability difference curves for various values of g and m = 1.8. The probability difference only changes in amplitude for a constant m and varying g. (d) The two probability difference curves used to model P(r/ls*) from Eq. (9).

Download Full Size | PDF

\begin{array}{l} P_{g, m} = P_{g = 0} + c_{1} (g, m) P_{1} Δ + c_{2} (g, m) P_{2} Δ \\ P_{1} Δ = P_{g = 0.9, m = 1.5} - P_{g = 0} \\ P_{2} Δ = P_{g = 0.9, m = 1.01} - P_{g = 0} \end{array}

The coefficients c₁ and c₂ vary smoothly and continuously with g and m. These coefficients can be fit to a variety of functions, depending on the desired simplicity and accuracy of the model. We fit these coefficients to a third order polynomial in two dimensions described by Eq. (10).

c_{i} = a_{i} + b_{i} x + c_{i} y + d_{i} x^{2} + e_{i} y^{2} + f_{i} x^{3} + g_{i} y^{3} + h_{i} x y + i_{i} x^{2} y + j_{i} x y^{2},

where x = ln(m) and y = 1/(1-g). The constants, a – j, are supplied in Table 1 . These constants were optimized to obtain a minimized error for g≥0.6.

Table 1. Values of constants from coefficients for PΔ model [Eq. (9)]

View Table

Alternatively, a principle component model can also be adopted similar to the one described by Eq. (8), except that the coefficients c₁, c₂, and c₃ each vary as a function of g and m in the generalized model. The variation of these coefficients is also smooth and continuous and can be fit to a polynomial equation based model such as the one in Eq. (10). Fig. 6(a) shows a comparison of the difference model (PΔ model) for the Whittle-Matérn phase function with Monte Carlo results for varying values of m and a constant g of 0.9. The agreement is excellent, although the error slightly increases for larger values of m. The agreement is improved for the PCA based model, shown in Fig. 6(b). The error is quantified for the entire range of g and m for the PΔ and PCA models in Fig. 6(c) and Fig. 6(d), respectively. Although the equations were optimized for g≥0.6, the average error for r/ls* between 0 and 1 is less than 2% for the entire range of g and m values that were evaluated. The error was less than 1% for all values of m and biologically relevant anisotropy factors (g≥0.6). The PCA model had less error than the PΔ model for this biologically relevant anisotropy range.

Fig. 6 Evaluation of P(r/ls*) predictive model accuracy. (a) PΔ model for P(r/ls*) (lines) compared to Monte Carlo data (points) for varying values of m and g = 0.9. (b) PCA model for P(r/ls*) (lines) compared to Monte Carlo data (points) for varying values of m and g = 0.9. (c) Average error of PΔ model as a function of g for r/ls*<1. (d) Average error of PCA model as a function of g for r/ls*<1.

Download Full Size | PDF

It is important to point out that the presented model does not incorporate absorption. However, the absorption is expected to be much smaller at small radial distances because the path length and the exit radius of a ray traversing a turbid medium are known to be well correlated. To estimate the range of absorption under which the presented models are accurate, we performed three Monte Carlo simulations with identical optical properties but a varying absorption coefficient. In this case, each of the backscattering distributions was normalized by the total backscattering reflectance that would be present in the absence of absorption. The resulting reflectance distributions are shown in Fig. 7 . Absorption primarily alters the intensity of backscattering at larger radial distances and has a minimal effect at r << l_a, where l_a = 1/µ_a. This can be reasoned from the Beer-Lambert law and the fact that the path lengths contributing to backscattering at small radial distances are small. For larger radial distances, the path length becomes longer and the effect of absorption becomes substantial. Therefore, one can also incorporate traditional diffusion approximation models of absorption in order to quantify the backscattering contribution. In this case, the isotropic scattering portion [P_{g = 0} from Eq. (7) to (9)] can be modeled with standard diffusion approximation equations for reflectance [18,19].

Fig. 7 P(r) for varying µ_a. The scattering properties of all three simulations were maintained at µ_s* = 100 cm⁻¹, g = 0.9, and m = 1.6. Absorption has a minimal effect at small radial distances because of the short path lengths involved.

Download Full Size | PDF

4. Experimental Measurement of P(r) from Tissue Phantom

As an illustrative example of an application of the presented models, we utilized Low-Coherence Enhanced Backscattering (LEBS) to measure P(r) for r<ls*. The experimental system used for the measurements has been described elsewhere and validated with P(r) measurements and simulations of mono-dispersed polystyrene microspheres in water [20]. A phantom was constructed from a mixture of polystyrene microsphere suspensions in order to mimic the Whittle-Matérn phase function. The proportions of microsphere sizes were determined by fitting three commonly available sizes such that the final phase function matched closely to the Whittle-Matérn phase function. Microspheres that were sized 0.26, 0.52, and 1.0µm were mixed in respective volume proportions of 1:0.2:0.7. When the suspension is dilute and the scattering is weak, the weighted average of the phase functions (with the weights being proportional to the concentration and the scattering cross section of the microspheres) results in the effective phase function of the scattering medium. The phase functions for microsphere particles and their scattering cross sections were calculated with Mie theory. The phase function fit was further improved by raising the refractive index of the medium to 1.34 with glycerol. A Whittle-Matérn phase function with an m of 1.6 and g of 0.8 is compared to the resulting phase function obtained from the microsphere mixture in Fig. 8(a) . In principle, the agreement may further be improved by using a larger variety of sphere sizes. We then obtained measurements for the unpolarized P(r) by measuring the unpolarized LEBS signal and dividing the Fourier transform of the signal by the coherence function [20]. In this example, we approximated the coherence function to be the shape of an idealized first order Bessel function of the first kind with a spatial coherence length of 110µm, as calculated from the van-Cittert-Zernike theorem. We accounted for the decrease in the enhancement factor due to the contribution of orthogonally polarized light by assuming that unpolarized light is an equal mixture of the two orthogonal polarization states. This results in an enhancement of 1.75 when dividing by the incoherent baseline, or 0.75 when the baseline is subtracted [21]. We also accounted for the discrepancy between Monte Carlo and experimental measurement that was observed in earlier work with a factor of 0.5 [20]. These constants were determined to be accurate by measuring P(r) from mono-disperse microsphere suspensions with varying sizes and comparing to Monte Carlo simulations that track polarization and utilized the Mie phase function. The experimentally measured P(r) and the P(r) obtained from the PCA model were in excellent agreement for r << ls* [Fig. 8(b)], with an average error of 7.1% for the range shown. It is interesting to note that the noise level from the experimental measurement in Fig. 8(b) varies with the radial position. The LEBS method measures the two-dimensional quantity p(x,y) via a two-dimensional Fourier transform. The higher noise level at small radii is a result of the decreased number of pixel elements that contribute to the summation over the polar angle. As r increases, the noise level decreases until the coherence function begins to approach small values, at which point the noise level quickly rises due to the division of two numbers that have values near 0.

Fig. 8 Tissue phantom for experimental measurement of P(r). (a) Comparison of the Whittle-Matérn phase function (m = 1.6 and g = 0.8) with the fit obtained from a mixture of microspheres of three sizes. (b) Resulting P(r) experimentally measured with LEBS, compared with P(r) obtained with the PCA based model and a Monte Carlo simulation. An ls* of 800μm, determined from the Mie theory calculation, was used for the model and Monte Carlo curve. The experimental LEBS measurement was scaled by the 0.75 factor to account for unpolarized light and the empirical 0.5 factor, determined from microsphere studies in prior work [20], to account for the difference observed between LEBS and Monte Carlo (see text). The average error between the experimentally measured P(r) and the model was 7.1%.

Download Full Size | PDF

5. Conclusions

We have presented two models capable of predicting P(r/ls*) for the entire range of r/ls*. Each model was applied to predicting reflectance from random media with a Henyey-Greenstein phase function and a generalized phase function that is derived from the Whittle-Matérn correlation function. The presented models consist of an empirical equation that quantifies the difference between the isotropic case of g = 0 and non-isotropic cases of g>0. This difference is predicted with two methods: the PΔ method and the PCA method. The PΔ method is based on the observation that these differences only scale in amplitude as g is varied [Fig. 2(c)] and do not change their shape. Although the shape of the radial dependence of this anisotropy contribution depends on m, it can be predicted with a linear combination of two difference terms [Eq. (9)]. The PΔ model results in accurate prediction for the Henyey-Greenstein case [Fig. 3(c)] and the generalized Whittle-Matérn phase case [Fig. 6(a)]. The PΔ(r/ls*) terms used in this model can be thought of as the anisotropic contribution to P(r/ls*), with the amplitude of these terms being related to the anisotropy factor.

The PCA method uses a similar idea as the PΔ method in that the anisotropy contribution is quantified. The advantage of the PCA method is that it decomposes the radial contributions into orthogonal principle components that explain the largest variations in the Monte Carlo data. Although, in principle, using more principle components results in a more accurate prediction of the data, in practice, accurate prediction also relies on being able to model the coefficients associated with the principle components. The dependence of these coefficients with the phase function properties become increasingly complicated for higher principle components [Fig. 4(d)], and therefore the model that fits these coefficients requires increasing complexity. We used a polynomial model for the Henyey-Greenstein phase function, and a two-dimensional cubic [Eq. (10)] for the Whittle-Matérn phase function. This resulted in excellent agreement for the biological tissue regime where g>0.6 [Fig. 4(a) and Fig. 6(d)]. The results of the fits can potentially be improved by implementing more principle components, increasing the accuracy of the model fit to the coefficients, or fitting to a smaller range of optical properties. That said, the error for the model applied to biologically relevant optical properties (g≥0.6 and 1<m<2) was less than 1%. The same procedure that was presented here can be used for modeling other ranges of g and m in order to obtain improved accuracy.

As mentioned in section 3, absorption is not included in the models described in this manuscript. However, in biological tissue, absorption varies dramatically with wavelength and is typically small for λ>600nm. Therefore, the technique described here can be utilized to measure scattering properties in the non-absorbing wavelength regions. Absorption can then be characterized by understanding the path length distribution for varying optical properties and measuring the backscattering for varying wavelengths. From Fig. 7, we can conclude that absorption primarily alters the intensity of backscattering at larger radial distances and has a minimal effect at r << l_a. This is due to shorter path lengths at smaller radial distances resulting in less attenuation of the scattered rays (The Beer-Lambert law). In cases where absorption cannot be neglected, a traditional diffusion approximation model of absorption in order to quantify the backscattering contribution can be used. In this case, the isotropic scattering portion [P_{g = 0} from Eq. (7) to (9)] can be modeled with standard diffusion approximation equations for reflectance [18,19].

In conclusion, the models presented in this work allow for accurate prediction of the impulse response function, P(r), to a random medium with a tissue-relevant range of optical properties and without the need for performing a large number of Monte Carlo simulations. Only three simulations are required including a simulation for isotropic scattering and two simulations for anisotropic scattering (g = 0.9 with m of 1.5 and 1.01). A Henyey-Greenstein based P(r) model is simpler in that it only requires two Monte Carlo simulations; however, it may not be as comprehensive of a model for tissue characterization. Finally, we presented a methodology for obtaining phantoms that have the potential to closely mimic optical properties of tissue, including the backscattering at small length-scales. The ability to predict the backscattering distribution at subdiffusion length scales holds promise for using techniques such as LEBS to measure optical properties of tissue (such as g, m and ls*) by measuring P(r). These results may also allow for faster, simpler and more accurate solutions to the inverse problem of measuring optical properties from tissue by providing an alternative for existing inverse Monte Carlo methods [5,6,11,22,23]. The three simulation and coefficient equations necessary for predicting P(r) will be made available online for public use. Furthermore, there are currently no existing empirical or theoretical models that allow for the prediction of the backscattered light at subdiffusion length scales without the need for performing repetitive and time intensive Monte Carlo simulations. The high degree of accuracy of the presented models and experimental illustration of a P(r) measurement from the Whittle-Matérn phase function at r<ls* indicate that the presented models and experimental phantom will be useful for characterizing the optical properties of biological samples.

Acknowledgements

This work was supported by National Institutes of Health (NIH) grants R01 CA128641, R01 CA109861, R01 EB003682, and R01 CA118794.

References and links

1. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef] [PubMed]

2. E. L. Hull and T. H. Foster, “Steady-state reflectance spectroscopy in the P-3 approximation,” J. Opt. Soc. Am. A 18(3), 584–599 (2001). [CrossRef]

3. I. Seo, C. K. Hayakawa, and V. Venugopalan, “Radiative transport in the delta-P1 approximation for semi-infinite turbid media,” Med. Phys. 35(2), 681–693 (2008). [CrossRef] [PubMed]

4. M. C. Skala, G. M. Palmer, C. F. Zhu, Q. Liu, K. M. Vrotsos, C. L. Marshek-Stone, A. Gendron-Fitzpatrick, and N. Ramanujam, “Investigation of fiber-optic probe designs for optical spectroscopic diagnosis of epithelial pre-cancers,” Lasers Surg. Med. 34(1), 25–38 (2004). [CrossRef] [PubMed]

5. A. Amelink, H. J. C. M. Sterenborg, M. P. L. Bard, and S. A. Burgers, “In vivo measurement of the local optical properties of tissue by use of differential path-length spectroscopy,” Opt. Lett. 29(10), 1087–1089 (2004). [CrossRef] [PubMed]

6. R. Reif, O. A’Amar, and I. J. Bigio, “Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media,” Appl. Opt. 46(29), 7317–7328 (2007). [CrossRef] [PubMed]

7. Y. L. Kim, Y. Liu, V. M. Turzhitsky, H. K. Roy, R. K. Wali, and V. Backman, “Coherent backscattering spectroscopy,” Opt. Lett. 29(16), 1906–1908 (2004). [CrossRef] [PubMed]

8. V. M. Turzhitsky, A. J. Gomes, Y. L. Kim, Y. Liu, A. Kromine, J. D. Rogers, M. Jameel, H. K. Roy, and V. Backman, “Measuring mucosal blood supply in vivo with a polarization-gating probe,” Appl. Opt. 47(32), 6046–6057 (2008). [CrossRef] [PubMed]

9. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13(12), 4420–4438 (2005). [CrossRef] [PubMed]

10. L. H. Wang, S. L. Jacques, and L. Q. Zheng, “Mcml - Monte-Carlo Modeling of Light Transport in Multilayered Tissues,” Comput Meth. Prog Biol. 47(2), 131–146 (1995). [CrossRef]

11. F. Bevilacqua and C. Depeursinge, “Monte Carlo study of diffuse reflectance at source-detector separations close to one transport mean free path,” J. Opt. Soc. Am. A 16(12), 2935–2945 (1999). [CrossRef]

12. C. J. R. Sheppard, “Fractal model of light scattering in biological tissue and cells,” Opt. Lett. 32(2), 142–144 (2007). [CrossRef] [PubMed]

13. J. D. Rogers, I. R. Capoğlu, and V. Backman, “Nonscalar elastic light scattering from continuous random media in the Born approximation,” Opt. Lett. 34(12), 1891–1893 (2009). [CrossRef] [PubMed]

14. I. R. Çapoğlu, J. D. Rogers, A. Taflove, and V. Backman, “Accuracy of the Born approximation in calculating the scattering coefficient of biological continuous random media,” Opt. Lett. 34(17), 2679–2681 (2009). [CrossRef] [PubMed]

15. A. Ishimaru, Wave propagation and scattering in random media (Academic Press, New York, 1978).

16. P. Guttorp and T. Gneiting, “Studies in the history of probability and statistics XLIX On the Matern correlation family,” Biometrika 93(4), 989–995 (2006). [CrossRef]

17. J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part II,” Opt. Express 13(25), 10392–10405 (2005). [CrossRef] [PubMed]

18. A. Kienle, L. Lilge, M. S. Patterson, R. Hibst, R. Steiner, and B. C. Wilson, “Spatially resolved absolute diffuse reflectance measurements for noninvasive determination of the optical scattering and absorption coefficients of biological tissue,” Appl. Opt. 35(13), 2304–2314 (1996). [CrossRef]

19. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef] [PubMed]

20. V. Turzhitsky, J. D. Rogers, N. N. Mutyal, H. K. Roy, and V. Backman, “Characterization of Light Transport in Scattering Media at Subdiffusion Length Scales with Low-Coherence Enhanced Backscattering,” IEEE J. Sel. Top. Quantum Electron. 16(3), 619–626 (2010). [CrossRef]

21. E. Akkermans, P. E. Wolf, and R. Maynard, “Coherent backscattering of light by disordered media: Analysis of the peak line shape,” Phys. Rev. Lett. 56(14), 1471–1474 (1986). [CrossRef] [PubMed]

22. G. M. Palmer and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part I: Theory and validation on synthetic phantoms,” Appl. Opt. 45(5), 1062–1071 (2006). [CrossRef] [PubMed]

23. L. V. Wang, “Rapid modeling of diffuse reflectance of light in turbid slabs,” J. Opt. Soc. Am. A 15(4), 936–944 (1998). [CrossRef] [PubMed]

Constant	C₁	C₂
a	−0.29773	−0.064755
b	0.11547	−0.46123
c	0.068216	0.15039
d	1.3086	1.5305
e	−0.0045095	−0.0049453
f	−1.1378	−1.0031
g	7.4495·10⁻⁵	3.9844·10⁻⁵
h	0.25636	−0.37841
i	0.047610	0.13420
j	−0.0071922	0.0082885

A predictive model of backscattering at subdiffusion length scales

Abstract

1. Introduction

1. Monte Carlo Simulation of Reflectance

2. Model of Reflectance: Henyey-Greenstein Phase Function

3. Model of Reflectance: Whittle-Matérn Phase Function

4. Experimental Measurement of P(r) from Tissue Phantom

5. Conclusions

Acknowledgements

References and links

Cited By

Figures (8)

Tables (1)

Equations (10)

Biomedical Optics Express