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

**|**HHS Author Manuscripts**|**PMC3370954

Formats

Article sections

Authors

Related links

Nat Commun. Author manuscript; available in PMC 2012 December 13.

Published in final edited form as:

PMCID: PMC3370954

NIHMSID: NIHMS352343

Edward Vitkin,^{‡} Vladimir Turzhitsky,^{‡} Le Qiu,^{‡} Lianyu Guo, Irving Itzkan, Eugene B. Hanlon, and Lev T. Perelman^{*}

Center for Advanced Biomedical Imaging and Photonics Department of Obstetrics, Gynecology and Reproductive Biology and Department of Medicine Beth Israel Deaconess Medical Center Biological and Biomedical Sciences Program Harvard University Boston, Massachusetts 02215 USA

The publisher's final edited version of this article is available at Nat Commun

See other articles in PMC that cite the published article.

From astronomy to cell biology, the manner in which light propagates in turbid media has been of central importance for many decades. However, light propagation near the point-of-entry (POE) in turbid media has never been analytically described, until now. Here we report a straightforward and accurate method that overcomes this longstanding, unsolved problem in radiative transport. Our theory properly treats anisotropic photon scattering events and takes the specific form of the phase function into account. As a result, our method correctly predicts the spatially dependent diffuse reflectance of light near the POE for any arbitrary phase function. We demonstrate that the theory is in excellent agreement with both experimental results and Monte Carlo simulations for several commonly used phase functions.

Recent advances in optical spectroscopy of human epithelium, high resolution imaging and early cancer detection demand an accurate description of light propagation near the point where light enters the tissue^{1–8}. However, an analytical method capable of describing the propagation of photons near their point of entry (POE) in turbid media has eluded researchers in diverse fields for generations^{9–12}. Propagation of electrons in semiconductors^{13}, anisotropic light diffusion in liquid crystals^{14}, and coherent backscattering^{15} are only a few examples that can benefit from such a method. The most accurate treatment of such problems, when interference can be neglected, uses the radiative transport equation,

$$\mathbf{s}\cdot \nabla L\left(r,s\right)=-({\mu}_{a}+{\mu}_{s})L\left(r,s\right)+{\mu}_{s}\int \int p\left(s,{s}^{\prime}\right)L\left(r,{s}^{\prime}\right)d{s}^{\prime}+S\left(r,{s}^{\prime}\right)$$

(1)

where *L*(**r**, **s**) is the radiance of light at position **r** traveling in the direction of unit vector **s**, *S*(**r**, **s**) describes the distribution of internal sources, and μ_{a} and μ_{s} are the absorption and scattering coefficients of the medium, respectively. The unity normalized function *p*(**s, s**′) is the phase function of the medium, which describes the angular distribution for a single scattering event (here, **s**′ is the incoming and **s** the outgoing direction of scatter)^{9}.

Unfortunately, no analytical solution to this equation has been discovered for practical scenarios. An often used numerical method is Monte Carlo simulation. Monte Carlo simulation gives accurate results but does not provide any physical insight and often requires excessive computational time. Recent GPU based Monte Carlo methods significantly reduce computational time^{16–19} thus extending the utility of Monte Carlo simulations to practical applications, but inverse problems with a large number of parameters will greatly benefit from the analytical approach provided here. In the absence of suitable numerical methods, the most widely used approach is to obtain an approximate solution using the diffusion approximation (DA) to the radiative transport equation. The DA works best when scattering is isotropic or near isotropic^{20}. It is not sensitive to the particular form of the phase function, reducing phase function information to a single parameter, $g=\int p(\mathbf{s}\cdot {\mathbf{s}}^{\prime})(\mathbf{s}\cdot {\mathbf{s}}^{\prime})d{\mathbf{s}}^{\prime}$, called the average cosine of the phase function or the anisotropy factor, which in turn is incorporated into the reduced scattering coefficient, ${\mu}_{s}^{\prime}={\mu}_{s}(1-g)$. Hence, ${\mu}_{s}^{\prime}$ becomes the only parameter in the DA to describe scattering. While the diffusion approximation can provide analytical solutions in some cases, it does not provide accurate results for the critically important region near the POE.

This limitation of the DA has been a major impediment, particularly for light scattering in biological media. Many investigators have attempted to amend the DA by solving the radiative transport equation using more complicated phase functions^{21–23}. One such approach, called the δ-*P*_{1} approximation, is to solve the equation for the delta-Eddington phase function, which results in some improvements for characterizing forward scattering media^{24}. Another approach is to decompose the phase function into a series of Legendre polynomials and to construct a system of radiative transport equations corresponding to each Legendre term. Because it becomes very cumbersome to solve the system of partial differential equations as the number of Legendre terms increases, the *P*_{3} approximation uses only the 0^{th} through 3^{rd} terms in the Legendre series^{25}. Both approximations are substantial improvements over the DA because they take the phase function into account and broaden the parameter space with which the approximation can be used^{26,27}. However, until now, no proposed analytical method has addressed the important regime of light scattering near the POE (e.g., at source-detector separations less than $1\u2215{\mu}_{s}^{\prime}$) well enough to be widely used as common practice.

Here, we introduce a novel approach that results in the first successful attempt to analytically describe spatially dependent diffuse reflectance near the POE in turbid media with arbitrary phase function and absorption. The phase function corrected diffusion approximation (PFC) takes advantage of the delta-isotropic phase function, described below, to derive a two-part expression for the diffuse reflectance. The first part is identical to the DA for a pencil beam illumination, while the second part approximates the changes in reflectance due to the specific form of the phase function. The resulting predictions are shown to be in excellent agreement with experiment and also with Monte Carlo simulations at all distances, including the critical area near POE.

We consider the problem of calculating the reflectance of light from a semi-infinite turbid medium with no internal sources (*S*(**r, s** = 0) illuminated with a narrow collimated beam (see Fig. 1), an important, general problem often treated using the DA.

The radiance in this problem can be divided into three main parts: an on-axis part due to the narrow collimated beam, a predominately isotropic, phase function independent part and an anisotropic, phase function dependent part. Similarly, any arbitrary phase function can be expressed as a sum of a delta function, an isotropic part and an anisotropic part. We achieve this by expressing an arbitrary phase function, *p*(**s**· **s**′), as the sum of the delta-isotropic phase function with the same anisotropy factor and a deviation term

$$p\left(s\cdot {s}^{\prime}\right)={p}_{\mathit{DI}}\left(s\cdot {s}^{\prime}\right)+[p\left(s\cdot {s}^{\prime}\right)-{p}_{\mathit{DI}}\left(s\cdot {s}^{\prime}\right)].$$

(2)

where *p _{DI}*(

By substituting equation (2) into equation (1) and using the reduced scattering coefficient, we can express equation (1) as

$$S\cdot \nabla L\left(r,s\right)=-({\mu}_{a}+{\mu}_{s}^{\prime})L\left(r,s\right)+\frac{{\mu}_{s}^{\prime}}{4\pi}\int \int L\left(r,{s}^{\prime}\right)d{s}^{\prime}+{\mu}_{s}\int \int [p\left(s\cdot {s}^{\prime}\right)-{P}_{DI}\left(s\cdot {s}^{\prime}\right)]L\left(r,{s}^{\prime}\right)d{s}^{\prime}.$$

(3)

Without the second integral on the right-hand side, equation (3) would be the equation of radiative transport for the case of an isotropically scattering medium, which, as mentioned above, can be described using the DA. Thus, indeed, it is logical to present the radiance *L*(**r**, **s**) as the sum of a predominantly isotropic diffuse part, *L _{d}*(

After substitution for *L*(**r**, **s**) in equation (3), we seek an expression for *L _{c}*(

The transport equation (3) can then be written as a system of two equations for the remaining parts of the radiance *L _{d}*(

The first equation,

$$s\cdot \nabla {L}_{d}\left(r,s\right)=-({\mu}_{a}+{\mu}_{s}^{\prime}){L}_{d}\left(r,s\right)+\frac{{\mu}_{s}^{\prime}}{4\pi}\int \int {L}_{d}\left(r,{s}^{\prime}\right)d{s}^{\prime}+{P}_{0}\frac{{\mu}_{s}^{\prime}}{8{\pi}^{2}}\frac{\delta \left(\rho \right)}{\rho}\text{exp}[-({\mu}_{a}+{\mu}_{s}^{\prime})z]$$

(4)

leads to the standard diffusion approximation and the second equation,

$$s\cdot \nabla {L}_{p}\left(r,s\right)=-({\mu}_{s}+{\mu}_{s}^{\prime}){L}_{p}\left(r,s\right)+{P}_{0}\frac{{\mu}_{s}}{2\pi}\frac{\delta \left(\rho \right)}{\rho}\left[p\left(s\cdot \widehat{z}\right)-\frac{1-g}{4\pi}\right]\text{exp}[-({\mu}_{s}+{\mu}_{s}^{\prime})z],$$

(5)

retains all of the phase function dependent terms. Here we use ${\mu}_{s}\int \int \left[p\left({s\cdot s}^{\prime}\right)-{p}_{DI}\left({s\cdot s}^{\prime}\right)\right]{L}_{d}(\mathbf{r},{\mathbf{s}}^{\prime})d{\mathbf{s}}^{\prime}=0$ and neglect the integral term ${\mu}_{s}\int \int [p\left({s\cdot s}^{\prime}\right)-{\left(2\pi \right)}^{-1}g\cdot \delta (1-{s\cdot s}^{\prime})]{L}_{p}\left({r,s}^{\prime}\right)d{\mathbf{s}}^{\prime}$ (The justification for the treatment of these two integrals appears in the Methods section).

Equation (4) exactly matches the equation of radiative transport for a turbid medium with an isotropic phase function, illuminated with a narrow collimated beam, and has been solved by a standard diffusion approximation procedure. This results in the radially dependent diffuse reflectance *R _{d}*(ρ) presented in

$${R}_{d}\left(\rho \right)=\frac{{\mu}_{s}^{\prime}}{4\pi}{\int}_{0}^{\infty}\left[z{r}_{1}^{-3}(1+{\mu}_{\mathit{eff}}{r}_{1})\text{exp}(-{\mu}_{\mathit{eff}}{r}_{1})+(z+2{z}_{b}){r}_{2}^{-3}(1+{\mu}_{\mathit{eff}}{r}_{2})\text{exp}(-{\mu}_{\mathit{eff}}{r}_{2})\right]\text{exp}[-({\mu}_{a}+{\mu}_{s}^{\prime})z]dz,$$

(6)

where ${\mu}_{\mathit{eff}}={\left[3{\mu}_{a}({\mu}_{a}+{\mu}_{s}^{\prime})\right]}^{1\u22152}$ is the effective attenuation coefficient, ${z}_{b}=2\u2215\left[3({\mu}_{a}+{\mu}_{s}^{\prime})\right]$ is the coordinate of the extrapolated boundary condition, ${r}_{1}=\sqrt{{\rho}^{2}+{z}^{2}},{r}_{2}=\sqrt{{\rho}^{2}+{(z+2{z}_{b})}^{2}}$, and ρ is the exit radius measured from the POE in the surface plane of the semi-infinite medium (see Fig. 1).

Equation (5) describes the phase function correction and is also integrable. The solution yields the phase function correction term, *R _{p}*(ρ), for the radially dependent reflectance, obtained by projecting

$${R}_{p}\left(\rho \right)={\mu}_{s}{\int}_{0}^{\infty}\left[p\left(\frac{-z}{{r}_{1}}\right)-\frac{1-g}{4\pi}\right]\text{exp}\left[-({\mu}_{a}+{\mu}_{s}^{\prime})(z+{r}_{1})\right]\frac{z}{{r}_{1}^{3}}dz.$$

(7)

.

The physical nature of the phase function correction (7) can be qualitatively understood if we notice that*z* + *r*_{1} is the distance a photon would travel through the turbid medium to the detector if it undergoes just a single scattering event, first traveling in the collimated beam to the depth *z*, then scattering at the angle cos^{−1} (−*z/r*_{1}) and finally traveling distance *r*_{1} to the detector (see Fig. 1). Moreover, the phase function in the integral (−*z/r*_{1}) represents the probability of scattering at an angle cos^{−1}. However, note that the attenuation (the exponential in equation (7)) has a characteristic inverse scattering length of ${\mu}_{s}^{\prime}$. The reduced scattering coefficient ${\mu}_{s}^{\prime}$ characterizes the attenuation of light due to diffuse scattering, i.e., only diffusely scattered photons are removed from the *z* + *r*_{1} path; unscattered and low-angle scattered photons remain. If we were considering only unscattered photons, the attenuation factor would have the scattering coefficient μ* _{s}* in the exponent, instead of the reduced scattering coefficient ${\mu}_{s}^{\prime}$. The

To check the accuracy of PFC, we compared the reflectance calculated by Monte Carlo simulation, which is considered to be the gold standard, to the reflectance calculated using PFC, i.e., the sum of the standard DA term, equation (6) and the phase function correction term, equation (7). The Monte Carlo program used was validated against a well accepted code^{32}. Following the procedure described in ^{ref. 33}, we obtained the 95% confidence interval of all the Monte Carlo simulations presented. The margin of error is within 1% of the mean for $\rho {\mu}_{s}^{\prime}<0.1$ and within 2% for the full range. We include results for two phase functions commonly used in turbid media: the Henyey-Greenstein phase function^{34}
${p}_{HG}\left({s\cdot s}^{\prime}\right)=[(1-{g}^{2})\u22154\pi ]\u2215{[1+{g}^{2}-2g\left({s\cdot s}^{2}\right)]}^{3\u22152}$ and the “Mie” phase function, which can be described using Mie theory^{35}. The Mie phase function was calculated for a Gaussian size distribution with a standard deviation of 10% and a mean size chosen such that the resulting phase function had the desired anisotropy factor, *g*. The results are presented in Figures 2 and and33.

Figure 2a shows the sensitivity of PFC to the particular choice of phase function. While the DA strongly disagrees with the results of the Monte Carlo simulations, especially in the vicinity of the source $(\rho {\mu}_{s}^{\prime}<0.1)$, PFC shows excellent agreement for both the Henyey-Greenstein and Mie phase functions. As an example, a comparison is presented for *g* 0.95, although for other values of the anisotropy factor *g* the agreement is equally good. Importantly, both the Henyey-Greenstein and Mie phase functions used in these simulations have identical anisotropy factors yet show marked differences at small source-detector separations. As demonstrated in Figure 2a, these phase function-dependent changes in reflectance are well characterized by PFC. A comparison between PFC and recent attempts to improve the DA is shown in Figure 2b, together with a plot of the percentage error between the approximations and Monte Carlo simulation. The PFC diffusion theory shows dramatically improved agreement with Monte Carlo simulations over both the δ-*P*_{1} and *P*_{3} approximations. For the case of the Mie phase function PFC shows similar results, though the comparison of δ-*P*_{1} and *P*_{3} approximations with Monte Carlo simulations demonstrate improvement compared to Henyey-Greenstein phase function close to the POE (Supplementary Fig. S1).

Figure 3 presents several scenarios where PFC corrects the major deficiencies of the DA. Scenarios of low $({\mu}_{a}\u2215{\mu}_{s}^{\prime}=0.01)$ and high $({\mu}_{a}\u2215{\mu}_{s}^{\prime}=1)$ absorption are shown in Figure 3a. The DA gives significant errors for $\rho {\mu}_{s}^{\prime}<0.1$, i.e. small distances from the source, however PFC does not suffer from this problem and agrees with the Monte Carlo simulations for all distances. For large source-detector separations and high absorptions, the amount of light emerging from the scattering medium is usually so small that detection is impractical. Figure 3b shows the effects of varying the anisotropy factor*g*. PFC agrees with Monte Carlo simulations for both *g* = 0 and *g* = 0.9 while the DA is only accurate for anisotropy factor *g* = 0, i.e., totally isotropic scattering. Percentage errors between theory and Monte Carlo for Fig. 3 are shown in Supplementary Fig S2. The inset in Figure 3b depicts the separate contributions of *R _{p}*(ρ) and

The correction is intended to work best near the POE where the specific form of the phase function plays a role, while the diffusion approximation is known to work well for large radii where the specific form of the phase function does not matter. The intermediate regime is described less accurately due to the neglected integral term. The average deviation in this range of $\rho {\mu}_{s}^{\prime}\u25fb0.21$ is less than 10%, which is significantly better than the other methods (see Fig. 2b).

Though the primary intention of this work is to present a solution to the problem of diffuse reflectance near the POE, we note that the correction also works within the medium. The approach presented above employs the standard diffusion approximation to provide the isotropic solution, which is then corrected with the PFC. However, the standard diffusion approximation is not particularly accurate in predicting the distribution of light near the POE within the medium, even for *g* = 0. Thus, to check the accuracy of the correction within the medium, we employed Monte Carlo simulation for both isotropic and anisotropic cases using the Henyey-Greenstein phase function and calculated the phase function related change in the fluence rate at several depths within the medium. The PFC fluence rate is obtained by integrating *L _{p}*(

$${\varphi}_{p}(\rho ,z)={\mu}_{s}{\int}_{0}^{\infty}\left[p\left(\frac{z-{z}^{\prime}}{{\rho}_{1}}\right)-\frac{1-g}{4\pi}\right]\text{exp}[-({\mu}_{a}+{\mu}_{s}^{\prime})({z}^{\prime}+{\rho}_{1})]\frac{d{z}^{\prime}}{{\rho}_{1}^{2}}$$

(8)

where ${\rho}_{1}=\sqrt{{\rho}^{2}+{(z-{z}^{\prime})}^{2}}$. The comparison with Monte Carlo simulations shows good agreement within the medium including the vicinity of the POE (see Fig. 4) at depths of $z{\mu}_{s}^{\prime}\le 0.5$. Accurate treatment of larger depths would require including transverse spreading of the collimated beam.

We also compared the reflectance calculated using PFC with experimental measurements performed using the setup shown in Fig. 5a (see Methods for details). A high numerical aperture (N.A.) objective was used in this setup because equations (6) and (7) describe a 2π collection solid angle and Monte Carlo simulations show the illumination N.A. has virtually no effect on the shape of the reflectance near the POE (Supplementary Fig. S3). The scattering medium consisted of a mixture of an aqueous glycerin solution and polystyrene microspheres with a diameter of 0.457 μm (Polysciences, Inc.). Using Mie theory^{35} and assuming a Gaussian size distribution of the microspheres with a standard deviation of 0.011 μm, as listed by the manufacturer, we calculated the reduced scattering coefficient of the scattering sample to be ${\mu}_{s}^{\prime =5.6}$ mm^{−1} and the anisotropy factor to be *g* = 0.89. The absorption coefficient of the scattering sample is approximately μ* _{a}* = 0.001 mm

Figure 6 shows the experimental reflectance spectra and a comparison with the DA and PFC theories. While the DA exhibits significant deviation from the experiment for distances smaller than 100 μm from the POE, PFC shows excellent agreement.

One of the major applications of this theory will be solving the inverse problem to obtain optical properties of unknown turbid samples especially on a significantly smaller scale than it can be done presently. Fortunately, the same inversion techniques^{36,37} that have been used with the standard diffusion approximation to extract ${\mu}_{s}^{\prime}$ and μ_{a} can be applied with the phase function correction, provided that adequate assumptions are made regarding the phase function. Aside from higher accuracy near the source, PFC also allows for the inversion of parameters of the phase function. For example, the anisotropy factor can be measured by assuming either Mie or Henyey-Greenstein scattering and using *g* as a free parameter in the fitting to experimental data. Moreover, the entire backscattering portion of the phase function can potentially be reconstructed.

A requirement of the technique is the availability of an accurate isotropic solution (*g* = 0). Though this solution is available for diffuse reflectance, an accurate solution for the fluence rate near the POE is not commonly available. One approach for obtaining the fluence rate within the medium would be to use Monte Carlo simulation for isotropic scattering and combine it with PFC to obtain the fluence for arbitrary phase functions. This is the approach we used for obtaining the results in Figure 6. Thus the result of a single Monte Carlo simulation can be used to obtain solutions for a large array of phase functions and anisotropy parameters. It is also likely that PFC could be used to describe angular dependences in a similar fashion.

In conclusion, PFC, the phase function corrected diffusion theory for light transport in turbid media, addresses the deficiencies of the standard diffusion theory near the POE by accounting for the specific form of the phase function. The technique accurately predicts the correct scattering behavior for two frequently used but very different types of phase functions, the Henyey-Greenstein and Mie theory phase functions (Fig. 2a). The PFC approach results in a substantial improvement over the predictions of the δ-*P*_{1} approximation and the *P*_{3} approximation (Fig. 2b), each of which seeks to improve the performance of the DA. PFC also demonstrates excellent agreement with experimental results (Fig. 6).

We anticipate that this work will be utilized extensively in characterizing photon scattering in turbid media such as human tissue and in a variety of other applications.

To evaluate the first integral term ${\mu}_{s}\int \int [p\left({s\cdot s}^{\prime}\right)-{p}_{DI}\left({s\cdot s}^{\prime}\right)]{L}_{d}\left({r,s}^{\prime}\right)d{\mathbf{s}}^{\prime}$ we use the standard diffusion approximation expression for *L _{d}* (

$${L}_{d}\left(r,s\right)=\frac{1}{4\pi}[{\varphi}_{d}\left(\mathbf{r}\right)-3{s\cdot j}_{d}\left(\mathbf{r}\right)]$$

(9)

where ϕ_{d} (**r**) is the fluence rate and ${\mathbf{j}}_{d}\left(\mathbf{r}\right)=-\frac{1}{3({\mu}_{a}+{\mu}^{\prime}s)}\nabla {\varphi}_{d}\left(\mathbf{r}\right)$ is the current is the current density^{39}. By substituting $\frac{1}{4\pi}[{\varphi}_{d}\left(\mathbf{r}\right)-3\mathbf{s}\cdot {\mathbf{j}}_{d}\left(\mathbf{r}\right)]$ in the integral, we get

$${\mu}_{s}\int \int [p\left(s\cdot {s}^{\prime}\right)-{p}_{DI}\left(s\cdot {s}^{\prime}\right)]{L}_{d}\left(r,{s}^{\prime}\right)d{s}^{\prime}=\phantom{\rule{0ex}{0ex}}=\frac{{\mu}_{s}}{4\pi}\int \int [p\left(s\cdot {s}^{\prime}\right)-{p}_{DI}\left(s\cdot {s}^{\prime}\right)]{\varphi}_{d}\left(\mathbf{r}\right)d{s}^{\prime}-\frac{3{\mu}_{s}}{4\pi}\int \int [p\left(s\cdot {s}^{\prime}\right)-{p}_{DI}\left(s\cdot {s}^{\prime}\right)]\left({{s}^{\prime}j}_{d}\left(r\right)\right)d{s}^{\prime}$$

(10)

From here, we see directly that $\frac{{\mu}_{s}}{4\pi}\int \int [p\left({s\cdot s}^{\prime}\right)-{p}_{DI}\left({s\cdot s}^{\prime}\right)]{\varphi}_{d}\left(\mathbf{r}\right)d{\mathbf{s}}^{\prime}=\frac{{\mu}_{s}{\varphi}_{d}\left(\mathbf{r}\right)}{4\pi}\int \int [p\left({s\cdot s}^{\prime}\right)-{p}_{DI}\left({s\cdot s}^{\prime}\right)]d{\mathbf{s}}^{\prime}=0$ due to the unity normalization of both phase functions.

Let's now simplify the integral term $\int \int p\left({s,s}^{\prime}\right)\left({s\cdot j}_{d}\left(\mathbf{r}\right)\right)d{\mathbf{s}}^{\prime}$ and re-write it in a spherical coordinate system (Fig. 7a). The integral can be expressed as μ_{s}|**j**_{d}(**r**)|∫∫ *p* (cosθ)(sinθ_{j} · sin sinθ + cosθ · cosθ_{j})sinθ*d*θ_{d}, where $\widehat{\mathbf{j}}$ is the unit vector in the direction of vector **j**_{d}(**r**), θ = c o^{−1}ss ·(**s**′), and ${\theta}_{j}=c\phantom{\rule{1em}{0ex}}{o}^{-1}s{\widehat{\mathbf{j}}}_{d}(\cdot \mathbf{s})$ The part proportional to the s i nθ_{j} s i n s · i n θ is zero since ${\int}_{0}^{2\pi}s\phantom{\rule{1em}{0ex}}\stackrel{.}{\phi}d\phi =0$ The part proportional to c o sθ c· o s θ_{j} gives |**j**_{d} (**r**)| cosθ_{j} ∫∫ *p*(cosθ)cosθsinθ*d*θ*d*. Here |**j**_{d}(**r**)|cosθ_{0} = s · **j**_{d} (**r**) and, by definition ∫∫*p*(cosθ)cosθ sin θ*d*θ*d**g*. Thus we see that ∫∫*p*(s, s′)(s′ · **j**_{d} (**r**))*d*s′ = *g* (s · **j**_{d} (**r**)). By substituting this into the remaining portion of the integral $\frac{3{\mu}_{s}}{4\pi}\int \int [p\left({s\cdot s}^{\prime}\right)-{p}_{DI}\left({s\cdot s}^{\prime}\right)]\left({{s}^{\prime}\cdot j}_{d}\left(\mathbf{r}\right)\right)d{\mathbf{s}}^{\prime}$ we get

$$\frac{3{\mu}_{s}}{4\pi}\int \int [p\left(s\cdot {s}^{\prime}\right)-{p}_{DI}\left(s\cdot {s}^{\prime}\right)]\left({{s}^{\prime}\cdot j}_{d}\left(\mathbf{r}\right)\right)d{s}^{\prime}=\frac{3{\mu}_{s}}{4\pi}[g\left({s\cdot j}_{d}\left(\mathbf{r}\right)\right)-g\left({s\cdot j}_{d}\left(\mathbf{r}\right)\right)]=0,$$

(11)

since both phase functions *p*(s · s′) and *p*_{DI} (s ·s′) have the same anisotropy factor *g*. Thus indeed μ* _{s}* ∫∫ [

To evaluate the contribution of the second integral term μ_{s} ∫∫[*p*(s · s′)−(2π)^{−1 }*g* · δ (1−s · s′)]*L*_{p} (**r**, s′)*d*s′, we calculated the phase function correction to the radially dependent reflectance including this integral. Then we subtracted *R _{p}* (ρ) obtained by neglecting this integral to calculate Δ

A schematic of the experimental set-up is shown in Figure 5a. The set-up employs a confocal microscope (FluoView 1000, Olympus) with several modifications. The incident 458 nm laser beam is expanded to fill the back aperture of the objective and passes through a 20/80 beamsplitter (Chroma). The beam is delivered to the scattering sample through a planar apochromatic oil immersion objective (N.A.=1.32, 63× magnification, Olympus). Light scattered from the sample is collected by the same objective and is reflected by a beamsplitter toward the imaging lens and the imaging CCD (DU-434-FI, Andor Technology). The spatial scale for imaging the scattering sample is derived by calibrating the CCD with a counting chamber (Millipore), with a single CCD pixel imaged to 0.327 μm on the sample.

The sample is placed in an imaging chamber assembled from a glass slide, four adhesive silicone isolators (JTR20R-A2-1.0, Grace Bio-Labs Inc.), and a No. 1 cover glass. The resulting chamber had a 20 mm internal diameter and a 4 mm depth. We placed optical gel (LS-3238, Fiber Optic Center) on the inside surface of the cover glass. After curing, the optical gel formed a layer of approximately 0.1 mm thickness with a refractive index of 1.388. The same refractive index of 1.388 at the 458 nm wavelength is achieved in the scattering sample by mixing a 66% volume fraction of the polystyrene microspheres aqueous solution and a 34% volume fraction of glycerin^{38}. By matching indices we minimize the effects of specular reflection and obtain a matched boundary condition at the interface.

The absorption coefficient of the sample can be estimated as a weighted average of the absorption coefficients of the glycerin, water, and polystyrene in the sample. The absorption coefficient of water^{40} at the wavelength 458 nm is 2.7·10^{−5} mm^{−1}. The absorption coefficient of glycerol^{41} is approximately 0.003 mm^{−1} and the absorption coefficient of polystyrene^{42} is less then 1·10^{−5} mm^{−1}. Thus the absorption coefficient of the scattering sample is approximately 0.001 mm^{−1}.

This work was supported by US National Science Foundation grants CBET-0922876 and CBET-0943180, US National Institutes of Health grants R01 EB003472, R01 EB006462 and R33 RR017361, and in part by the Department of Veterans Affairs, Office of Research and Development.

**Author contributions** E.V. and L.T.P. conceived and developed the theory, L.Q. and L.G. performed the experiments, V.T. and L.Q. performed the data analysis, I.I., E.B.H. and L.T.P. evaluated the method, E.V., V.T., L.Q., I.I., and E.B.H. contributed to the writing of the manuscript, and L.T.P. initiated and supervised the project and wrote the manuscript.

**Competing financial interests:** The authors declare no competing financial interests.

1. Qiu L, et al. Multispectral scanning during endoscopy guides biopsy of dysplasia in Barrett's esophagus. Nature Med. 2010;16:603–606. [PMC free article] [PubMed]

2. Backman V, et al. Detection of preinvasive cancer cells. Nature. 2000;406:35–36. [PubMed]

3. Gurjar RS, et al. Imaging human epithelial properties with polarized light-scattering spectroscopy. Nature Med. 2001;7:1245–1248. [PubMed]

4. Perelman LT, et al. Observation of periodic fine structure in reflectance from biological tissue: A new technique for measuring nuclear size distribution. Phys. Rev. Lett. 1998;80:627–630.

5. Reif R, A'Amar O, Bigio IJ. Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media. Appl. Opt. 2007;46:7317–7328. [PubMed]

6. Amelink A, Sterenborg H. Measurement of the local optical properties of turbid media by differential path-length spectroscopy. Appl. Opt. 2004;43:3048–3054. [PubMed]

7. Boas DA, Oleary MA, Chance B, Yodh AG. Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media -analytic solution and applications. Proc. Natl Acad. Sci. USA. 1994;91:4887–4891. [PubMed]

8. Itzkan I, et al. Confocal light absorption and scattering spectroscopic microscopy monitors organelles in live cells with no exogenous labels. Proc. Natl Acad. Sci. USA. 2007;104:17255–17260. [PubMed]

9. Chandrasekhar S. Radiative Transfer. Dover; New York: 1960.

10. Yoo KM, Liu F, Alfano RR. When does the diffusion approximation fail to describe photon transport in random media. Phys. Rev. Lett. 1990;64:2647–2650. [PubMed]

11. Blanco S, Fournier R. Short-path statistics and the diffusion approximation. Phys. Rev. Lett. 2006;97:230604-1–230604-4. [PubMed]

12. Chen NG, Bai J. Monte Carlo approach to modeling of boundary conditions for the diffusion equation. Phys. Rev. Lett. 1998;80:5321–5324.

13. Mal'shukov AG, Wang LY, Chu CS, Chao KA. Spin Hall effect on edge magnetization and electric conductance of a 2D semiconductor strip. Phys. Rev. Lett. 2005;95:146601-1–146601-4. [PubMed]

14. vanTiggelen BA, Maynard R, Heiderich A. Anisotropic light diffusion in oriented nematic liquid crystals. Phys. Rev. Lett. 1996;77:639–642. [PubMed]

15. Akkermans E, Wolf PE, Maynard R. Coherent backscattering of light by disordered media - analysis of the peak line-shape. Phys. Rev. Lett. 1986;156:1471–1474. [PubMed]

16. Alerstam E, Svensson T, Andersson-Engels S. Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration. J. Biomed. Opt. 2008;13:060504-1–060504-3. [PubMed]

17. Fang Q, Boas DA. Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units. Opt. Express. 2009;17:20178–20190. [PMC free article] [PubMed]

18. Ren N, et al. GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues. Optics Express. 2010;18:6811–6823. [PubMed]

19. Alerstram E, et al. Next-generation acceleration and code optimization for light transport in turbid media using GPUs. Biomed. Opt. Express. 2010;1:658–675. [PMC free article] [PubMed]

20. Yoon G, Prahl SA, Welch AJ. Accuracies of the diffusion approximation and its similarity relations for laser irradiated biological media. Appl. Opt. 1989;28:2250–2255. [PubMed]

21. Venugopalan V, You JS, Tromberg BJ. Radiative transport in the diffusion approximation: an extension for highly absorbing media and small source-detector separations. Phys. Rev. 1998;E 58:2395–2407.

22. Durian DJ, Rudnick J. Photon migration at short times and distances and in cases of strong absorption. J. Opt. Soc. Am. 1997;A 14:235–245.

23. Kienle A, Patterson MS. Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium. J. Opt. Soc. Am. 1997;A 14:246–254. [PubMed]

24. Seo I, Hayakawa CK, Venugopalan V. Radiative transport in the delta-P-1 approximation for semi-infinite turbid media. Med. Phys. 2008;35:681–693. [PMC free article] [PubMed]

25. Hull EL, Foster TH. Steady-state reflectance spectroscopy in the P-3 approximation. J. Opt. Soc. Am. 2001;A 18:584–599.

26. Menon S, Su Q, Grobe R. Determination of g and mu using multiply scattered light in turbid media. Phys. Rev. Lett. 2005;94:153904-1–153904-4. [PubMed]

27. Chin LCL, Lloyd B, Whelan WM, Vitkin IA. Interstitial point radiance spectroscopy of turbid media. J. Appl. Phys. 2009;105:102025-1–102025-11.

28. Groenhuis RAJ, Ferwerda HA, Tenbosch JJ. Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory. Appl. Opt. 1983;22:2456–2462. [PubMed]

29. Yoon G, Prahl SA, Welch AJ. Accuracies of the diffusion approximation and its similarity relations for laser irradiated biological media. Appl. Opt. 1989;28:2250–2255. [PubMed]

30. Zege EP, Ivanov AP, Katsev IL. Image Transfer Through a Scattering Medium. Springer-Verlag; New York: 1991.

31. Farrell TJ, Patterson MS, 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. 1992;19:879–888. [PubMed]

32. Wang LH, Jacques SL, Zheng LQ. MCML - Monte-Carlo modeling of light transport in multilayered tissues. Comput. Methods Programs Biomed. 1995;47:131–146. [PubMed]

33. Liu Q, Ramanujam N. Scaling method for fast Monte Carlo simulation of diffuse reflectance spectra from multilayered turbid media. J. Opt. Soc. Am. 2007;A 24:1011–1025. [PubMed]

34. Henyey LG, Greenstein JL. Diffuse radiation in the galaxy. Astrophys. J. 1941;93:70–83.

35. Toublanc D. Henyey-Greenstein and Mie phase functions in Monte Carlo radiative transfer computations. Appl. Opt. 1996;35:3270–3274. [PubMed]

36. Zonios G, et al. Diffuse reflectance spectroscopy of human adenomatous colon polyps in vivo. Appl. Opt. 1999;38:6628–6637. [PubMed]

37. Sun J, et al. Influence of fiber optic probe geometry on the applicability of inverse models of tissue reflectance spectroscopy: computational models and experimental measurements. Appl. Opt. 2006;45:8152–8162. [PubMed]

38. Fang H, et al. Confocal light absorption and scattering spectroscopic microscopy. Appl. Opt. 2007;46:1760–1769. [PubMed]

39. Wang LV, Wu HI. Biomedical Optics: Principles and Imaging. Wiley; New Jersey: 2007.

40. Hale GM, Querry MR. Optical constants of water in the 200-nm to 200-microm wavelength region. Appl. Opt. 1973;12:555–563. [PubMed]

41. Reich O, Lohmannsroben H-G, Schael F. Optical sensing with photon density waves: Investigation of model media. Phys. Chem. Chem. Phys. 2003;5:5182–5187.

42. Firbank M, Arridge SR, Schweiger M, Delpy DT. An investigation of light transport through scattering bodies with non-scattering regions. Phys. Med. Biol. 1996;41:767–783. [PubMed]

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