|Home | About | Journals | Submit | Contact Us | Français|
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 tissue1–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 generations9–12. Propagation of electrons in semiconductors13, anisotropic light diffusion in liquid crystals14, and coherent backscattering15 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,
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 time16–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 isotropic20. It is not sensitive to the particular form of the phase function, reducing phase function information to a single parameter, , called the average cosine of the phase function or the anisotropy factor, which in turn is incorporated into the reduced scattering coefficient, . Hence, 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 functions21–23. One such approach, called the δ-P1 approximation, is to solve the equation for the delta-Eddington phase function, which results in some improvements for characterizing forward scattering media24. 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 P3 approximation uses only the 0th through 3rd terms in the Legendre series25. 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 used26,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 ) 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
where pDI(s · s′)=[1− g +2g ·†(1−s·s′)]/4π is the delta-isotropic phase function28. This decomposition of the phase function enables us to separate the radiance into the three analogous parts described above.
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, Ld(r, s) a phase function dependent part, (Lp(r, s) and an on-axis part due to the narrow collimated beam, Lc(r, s).
After substitution for L(r, s) in equation (3), we seek an expression for Lc(r, s) which satisfies the boundary condition for a narrow collimated beam and a portion of equation (3), . Note that this portion is described by rather than μs and specifically accounts for low-angle scattered photons29 (see Fig. 1 for illustration). The boundary condition corresponding to a narrow collimated beam entering a semi-infinite medium is for inward directions (), where p0 is the power of the incoming beam, ρ and z are the cylindrical coordinates and is a unit vector along the z axis. This yields30 and the boundary conditions and for .
The transport equation (3) can then be written as a system of two equations for the remaining parts of the radiance Ld(r, s) and Lp(r, s). It is natural then to distribute the two remaining Lc(r, s) dependent terms so that the isotropic, phase function independent term is combined is combined with the equation for the predominately isotropic Ld(r, s) and the anisotropic, phase function dependent term is combined with the equation for the anisotropic Lp(r, s).
The first equation,
leads to the standard diffusion approximation and the second equation,
retains all of the phase function dependent terms. Here we use and neglect the integral term (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 Rd(ρ) presented in ref. 31,
where is the effective attenuation coefficient, is the coordinate of the extrapolated boundary condition, , 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, Rp(ρ), for the radially dependent reflectance, obtained by projecting Lp (r,s) on the surface of the medium and integrating over the emerging solid angle
The physical nature of the phase function correction (7) can be qualitatively understood if we notice thatz + r1 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/r1) and finally traveling distance r1 to the detector (see Fig. 1). Moreover, the phase function in the integral (−z/r1) 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 . The reduced scattering coefficient characterizes the attenuation of light due to diffuse scattering, i.e., only diffusely scattered photons are removed from the z + r1 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 . The z + r1 path is not the actual path of real photons, but rather the virtual path of quasi-ballistic photons which are described by equation (7). The contribution of the term describing this behavior can be positive or negative since it depends on the difference of the contributions of the actual phase function and the delta-isotropic phase function. Physically, the PFC term accounts for the contribution of these photons, which have undergone multiple low-angle scattering events plus a single large angle scattering event along the z + r1 path. That single large angle turn is described by the phase function, thus providing the phase function correction. An illustration of this scattering behavior is shown in Figure 1.
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 code32. 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 and within 2% for the full range. We include results for two phase functions commonly used in turbid media: the Henyey-Greenstein phase function34 and the “Mie” phase function, which can be described using Mie theory35. 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 , 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 δ-P1 and P3 approximations. For the case of the Mie phase function PFC shows similar results, though the comparison of δ-P1 and P3 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 and high absorption are shown in Figure 3a. The DA gives significant errors for , 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 factorg. 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 Rp(ρ) and Rd(ρ) and their asymptotic rates of decay for the g = 0.9 case. The contribution of the correction to the reflectance at becomes negligible, resulting in an accuracy that is determined by the particular diffusion theory model being used. Note that although we have chosen a well-accepted diffusion theory model to describe isotropic scattering, there may be other solutions to this component of the reflectance.
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 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 Lp(r,s) over the entire solid angle
where . The comparison with Monte Carlo simulations shows good agreement within the medium including the vicinity of the POE (see Fig. 4) at depths of . 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 theory35 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 mm−1 and the anisotropy factor to be g = 0.89. The absorption coefficient of the scattering sample is approximately μa = 0.001 mm−1 (see Supplementary Methods for details). The reflectance signal measured from the polystyrene microsphere sample is presented in Figure 5b. We note that the measurements were performed tens of micrometers from the POE, where the signal strength has significant variation. This distance is comparable to cellular sizes in biological tissue. Thus, measurements of the reflectance from tissue on this scale could be sensitive to alterations in cell architecture associated with normal function or early, pre-cancerous changes.
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 techniques36,37 that have been used with the standard diffusion approximation to extract 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 δ-P1 approximation and the P3 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 we use the standard diffusion approximation expression for Ld (r, s)
where ϕd (r) is the fluence rate and is the current is the current density39. By substituting in the integral, we get
From here, we see directly that due to the unity normalization of both phase functions.
Let's now simplify the integral term and re-write it in a spherical coordinate system (Fig. 7a). The integral can be expressed as μs|jd(r)|∫∫ p (cosθ)(sinθj · sin sinθ + cosθ · cosθj)sinθdθd, where is the unit vector in the direction of vector jd(r), θ = c o−1ss ·(s′), and The part proportional to the s i nθj s i n s · i n θ is zero since The part proportional to c o sθ c· o s θj gives |jd (r)| cosθj ∫∫ p(cosθ)cosθsinθdθd. Here |jd(r)|cosθ0 = s · jd (r) and, by definition ∫∫p(cosθ)cosθ sin θdθdg. Thus we see that ∫∫p(s, s′)(s′ · jd (r))ds′ = g (s · jd (r)). By substituting this into the remaining portion of the integral we get
since both phase functions p(s · s′) and pDI (s ·s′) have the same anisotropy factor g. Thus indeed μs ∫∫ [p(s · s′) − pDI (s · s′)]Ld (r, s′)ds′ = 0.
To evaluate the contribution of the second integral term μs ∫∫[p(s · s′)−(2π)−1 g · δ (1−s · s′)]Lp (r, s′)ds′, we calculated the phase function correction to the radially dependent reflectance including this integral. Then we subtracted Rp (ρ) obtained by neglecting this integral to calculate ΔRp (ρ). We then plot Rp (ρ) and can see it is just several percent of Rp (ρ for small (see Fig. 7b). We also include Rd (ρ) to show that at greater both corrections are neglectable. Thus neglecting μs ∫∫[p (s · s′)−(2π)−1 g · δ(1−s · s′)]Lp (r, s′)ds′ is justified.
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 glycerin38. 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 water40 at the wavelength 458 nm is 2.7·10−5 mm−1. The absorption coefficient of glycerol41 is approximately 0.003 mm−1 and the absorption coefficient of polystyrene42 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.