Home | About | Journals | Submit | Contact Us | Français |

**|**Biomed Opt Express**|**v.1(3); 2010 October 1**|**PMC3018048

Formats

Article sections

- 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
- References and links

Authors

Related links

Biomed Opt Express. 2010 October 1; 1(3): 1034–1046.

Published online 2010 September 30. doi: 10.1364/BOE.1.001034

PMCID: PMC3018048

Vladimir Turzhitsky,^{1,}^{*} Andrew Radosevich,^{1} Jeremy D. Rogers,^{1} Allen Taflove,^{2} and Vadim Backman^{1}

Received 2010 July 12; Revised 2010 August 27; Accepted 2010 September 26.

Copyright ©2010 Optical Society
of America

This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported License, which permits download and redistribution, provided that the original work is properly cited. This license restricts the article from being modified or used commercially.

This article has been cited by other articles in PMC.

We provide a methodology for accurately predicting elastic backscattering radial distributions from random media with two simple empirical models. We apply these models to predict the backscattering based on two classes of scattering phase functions: the Henyey-Greenstein phase function and a generalized two parameter phase function that is derived from the Whittle-Matérn correlation function. We demonstrate that the model has excellent agreement over all length scales and has less than 1% error for backscattering at subdiffusion length scales for tissue-relevant optical properties. The presented model is the first available approach for accurately predicting backscattering at length scales significantly smaller than the transport mean free path.

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.

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 \Delta {n}^{2}{\left(r/{l}_{c}\right)}^{m-3/2}{K}_{m-3/2}\left(r/{l}_{c}\right)$, where *Δn ^{2}* is the variance of the refractive index fluctuations,

$$\begin{array}{l}F(\mathrm{cos}\theta )=\frac{2\widehat{g}(m-1)\cdot {\left[1-2\widehat{g}\mathrm{cos}\theta +{\widehat{g}}^{2}\right]}^{-m}}{{(1-\widehat{g})}^{2-2m}-{(1+\widehat{g})}^{2-2m}}\\ \widehat{g}=1-\frac{\sqrt{1+4{\left(k{l}_{c}\right)}^{2}}-1}{2{(k{l}_{c})}^{2}}\\ k{l}_{c}=\frac{\sqrt{\widehat{g}}}{1-\widehat{g}},\end{array}$$

(1)

where *k = 2π/λ* and the normalization is such that $\int F\left(\mathrm{cos}\theta \right)}\cdot d\mathrm{cos}\theta =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(\theta )=\frac{2{(k{l}_{c})}^{2}(m-1)\cdot {\left[1+{\left(2k{l}_{c}\mathrm{sin}\left(\theta /2\right)\right)}^{2}\right]}^{-m}}{1-{(1+{(2k{l}_{c})}^{2})}^{1-m}}.$$

(2)

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=\langle \mathrm{cos}\theta \rangle =\{\begin{array}{l}\frac{{(1-\widehat{g})}^{2-2m}\cdot (1+{\widehat{g}}^{2}-2\widehat{g}(m-1))-{(1+\widehat{g})}^{2-2m}\cdot (1+{\widehat{g}}^{2}+2\widehat{g}(m-1))}{2\widehat{g}(m-2)\left[{(1+\widehat{g})}^{2-2m}-{(1-\widehat{g})}^{2-2m}\right]}\text{\hspace{1em}}m\ne 2\\ \frac{1+{\widehat{g}}^{2}}{2\widehat{g}}+\frac{\mathrm{ln}(1+\widehat{g})-\mathrm{ln}(1-\widehat{g})}{g\left[{(1-\widehat{g})}^{-2}-{(1+\widehat{g})}^{-2}\right]}\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}}m=2\end{array}$$

(3)

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:

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).

$$F(\mathrm{cos}\theta )=\{\begin{array}{l}\frac{\widehat{g}\cdot {\left[1-2\widehat{g}\mathrm{cos}\theta +{\widehat{g}}^{2}\right]}^{-1}}{\mathrm{ln}(1+\widehat{g})-\mathrm{ln}(1-\widehat{g})}\text{\hspace{1em}\hspace{1em}}m=1\\ \text{\hspace{1em}\hspace{1em}\hspace{1em}}1/2\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}}\widehat{g}=0\end{array}$$

(4)

In these two cases, *g* becomes:

$$g=\{\begin{array}{l}\frac{1+\widehat{g}}{2\widehat{g}\left[\mathrm{ln}(1+\widehat{g})-\mathrm{ln}(1-\widehat{g})\right]}\text{\hspace{1em}\hspace{1em}}m=1,\text{\hspace{0.17em}}\widehat{g}\ne 0\\ \text{\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{1em}\hspace{1em}}0\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}\hspace{0.17em}}\widehat{g}=0\end{array}.$$

(5)

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 *ξ*:

$$\mathrm{cos}\theta =\frac{1+{\widehat{g}}^{2}-{\left[\xi \left({(1-g)}^{2-2m}-{(1+g)}^{2-2m}\right)+{(1+g)}^{2-2m}\right]}^{\frac{1}{1-m}}}{2g}$$

(6)

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].

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)ds=1$, where s = r/ls* and $P(r)={\displaystyle \underset{0}{\overset{2\pi}{\int}}p(r,\theta )rd\theta =2\pi \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)dr=1$.

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 **...**

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)\left[{P}_{0.9}-{P}_{0}\right]\\ c(g)=a{g}^{b},\end{array}$$

(7)

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

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 **...**

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)PC1+{c}_{2}(g)PC2+{c}_{3}(g)PC3,$$

(8)

where *c _{1}*,

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:

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 **...**

$$\begin{array}{l}{P}_{g,m}={P}_{g=0}+{c}_{1}(g,m){P}_{1}\Delta +{c}_{2}(g,m){P}_{2}\Delta \\ {P}_{1}\Delta ={P}_{g=0.9,m=1.5}-{P}_{g=0}\\ {P}_{2}\Delta ={P}_{g=0.9,m=1.01}-{P}_{g=0}\end{array}$$

(9)

The coefficients *c _{1}* and

$${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}xy+{i}_{i}{x}^{2}y+{j}_{i}x{y}^{2},$$

(10)

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.

Alternatively, a principle component model can also be adopted similar to the one described by Eq. (8), except that the coefficients *c _{1}*,

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 **...**

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

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.

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 [

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.

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

1. Farrell T. J., Patterson M. S., Wilson B., “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).10.1118/1.596777 [PubMed] [Cross Ref]

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

3. Seo I., Hayakawa C. K., Venugopalan V., “Radiative transport in the delta-P1 approximation for semi-infinite turbid media,” Med. Phys. 35(2), 681–693 (2008).10.1118/1.2828184 [PMC free article] [PubMed] [Cross Ref]

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

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

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

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

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

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

10. Wang L. H., Jacques S. L., Zheng L. Q., “Mcml - Monte-Carlo Modeling of Light Transport in Multilayered Tissues,” Comput Meth. Prog Biol. 47(2), 131–146 (1995).10.1016/0169-2607(95)01640-F [PubMed] [Cross Ref]

11. Bevilacqua F., Depeursinge C., “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).10.1364/JOSAA.16.002935 [Cross Ref]

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

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

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

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

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

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

18. Kienle A., Lilge L., Patterson M. S., Hibst R., Steiner R., Wilson B. C., “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).10.1364/AO.35.002304 [PubMed] [Cross Ref]

19. Farrell T. J., Patterson M. S., Wilson B., “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).10.1118/1.596777 [PubMed] [Cross Ref]

20. Turzhitsky V., Rogers J. D., Mutyal N. N., Roy H. K., Backman V., “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).10.1109/JSTQE.2009.2032666 [PMC free article] [PubMed] [Cross Ref]

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

Articles from Biomedical Optics Express are provided here courtesy of **Optical Society of America**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |