The propagation of light through the biological tissue is a complicated process involving both absorption and scattering. The photon propagation model provides insight into the interaction between light and tissues and is essential for tomographic imaging with visible and near-IR light.1, 2 The radiative transfer equation (RTE) is considered the golden standard for biomedical applications.1, 2 Analytical solutions for the RTE are available for few simple geometries, and numerical solutions, such as the discrete ordinate method3, 4 and the spherical harmonic method,5 often lead to enormous computational cost, especially to solve inverse problems such as optical tomography,6 bioluminescence tomography, 7, 8, 9, 10 and fluorescence tomography.11, 12 Monte Carlo (MC) is a statistical simulation method in which the paths of photons are traced as they are scattered and absorbed within the medium.13 The MC method is well established to produce accurate estimates for light propagation in tissues. However, due to its statistical nature, the MC method has the disadvantage of requiring a long computation time; therefore, it is usually used as a reference method for other approaches. The popular diffusion approximation14, 15 (DA) to RTE is widely used in the field of biophotonics because of its high computational efficiency. Nevertheless, DA assumes weakly anisotropic scattering and works well only in a highly scattering and weakly absorbing medium.16 It is also not suitable for modeling light propagation in the visible spectrum in which biological tissue present significant photon absorption.17
In the RTE, the phase functions describe the scattering characteristics of the medium and significantly influence the precision of the solution and the efficiency of the computation. Because the exact form of the phase function is currently unknown, the popular Henyey-Greenstein (HG) function18 and the Delta-Eddington function19 are usually used to approximate the true phase functions in biomedical applications. These two functions can be written in closed forms with a single free parameter , called the anisotropic factor, which is often considered to be independent of the tissue scattering and absorption. Based on the Delta-Eddington phase function, a generalized diffusion model was presented to simulate photon propagation in highly absorbing medium and smaller source-detector separations.20 The inverse problem of optical parameters was presented based on the RTE with the Delta-Eddington function.21 The HG function was proven to be the most accurate in terms of the angular dependence of single scattering events in biological tissues.2, 13 However, the RTE with the HG phase function is difficult to simplify further. In this paper, we present a generalized Delta-Eddington function to approximate the real phase function. Based on this new definition of Delta-Eddington phase function, the RTE can be reduced to a system of integral equations with respect to both the photon fluence rate and the flux vector. The solution of the system of integral equations enables a highly accurate prediction for the photon propagation in biological tissue over a wide range of optical parameters of biomedical interest.
Phase Approximation Modelis the weight factor measuring the anisotropy of the photon scattering, which we call the anisotropy weight; and is a asymmetry factor of the phase function to modulate the weakly anisotropic scattering. The phase function is a linear combination of the weakly anisotropic scattering and the strongly peaked forward scattering. The original Delta-Eddington phase function rigidly defines the parameters and as and , respectively.2, 19 These parameters are related to a single free parameter , where is the anisotropic factor defined as the mean of the cosine of the scattering angles. Here, the proposed generalized Delta-Eddington function defines that the anisotropy weight and asymmetry factor are related to the photon absorbing and scattering coefficients in the medium. The optical properties of the medium are characterized by the anisotropy weight and asymmetry factor along with the absorption and scattering coefficients.
Based on the generalized Delta-Eddington phase function [Eq. 1], the RTEis the region of interest, is the photon radiance (in ), is the isotropic source (in ), is the scattering coefficient (in ), and is the absorption coefficient (in ), , and . The photon fluence rate and photon flux vector are, respectively, defined by3 is a linear, first-order differential equation for the photon propagation in a heterogeneous medium, and the radiance can be formulated as23 is a scalar so that , representing the distance from point to the boundary along the direction . Integrating Eq. 5 over all the solid angles, we have5 by the unit vector , and integrating Eq. 5 over all the solid angles, we obtain an equation for the flux vector: , in Eqs. 6, 7 would vanish for the matched refractive indices on the boundary. However, the refractive index in the biological tissues is higher than the refractive index of the surrounding air, and photons will be internally reflected at the boundary. In this case, describes the reflected radiance from a fraction of photons reflected on the boundary , and can be expressed as is the outward unit normal at on , and the internal reflection coefficient can be calculated14, 16 by , . From the first-order approximation of the radiance , the reflected radiance on the boundary can be approximated by9 into Eqs. 6, 7 and performing a variable transformation from the polar coordinates to the Cartesian coordinates , we obtain a system of integral equations with respect to the photon fluence rate and flux vector: 10 the PA equation because it is derived from an approximate phase function. The PA equation is a well-posed system of integral equations of the second kind, and it enables an accurate prediction to the photon fluence rate and flux over a wide range of optical coefficients. The photon fluence rate presents the distribution of photon density in the spatial domain, while the photon flux describes the axial photon energy transfer and is directly related to experimental measurement data on the boundary. The exiting photon flux at on the surface of object can be linked with the photon flux,2 , , and (or, equivalently, , , , and ) in the PA equation represent the optical parameters of the medium. When the absorption coefficient and scattering coefficient in the medium are unknown, we can directly determine the optical parameters , , and using optical tomography techniques based on the PA equation.6 If the absorption coefficient , the scattering coefficient , and conventional anisotropic factor in the tissue are known, based on a simple homogenous numerical phantom (such as a spherical or cylindrical object) with a known light source setting, the exiting photon flux on the boundary of the phantom can be generated using MC simulation with a appropriate phase function, such as the HG phase function. Then, an optimization procedure can be performed to determine the anisotropy weight and asymmetry factor by matching the exiting photon flux on the phantom boundary obtained from the PA equation to the MC-simulated counterpart. These exiting photon fluxes distribute on the boundary of the phantom and reflect the variation with different distances from the light source. Hence, the anisotropy weight and asymmetry factor are independent of the phantom geometry.
To perform the numerical computation based on Eq. 10, the region of interest can be discretized into finite elements with vertex nodes, and the photon fluence rate and flux are approximated in terms of nodal-based basis functions24 ,12 into Eq. 10, we obtain a system of linear equations with respect to the photon fluence rate and flux: , , and represent the corresponding discrete integral kernels and source vector in Eq. 10.
Extensive numerical experiments were conducted to validate the accuracy of the PA equation by MC simulation. The experiments are based on a cylindrical phantom of diameter and height. We set up a spherical source of radius with a power of at position of in the phantom. The phantom was then discretized into 25,335 tetrahedral elements and 4833 nodes, as shown in Fig. 1 . A total of 1226 virtual detectors were allocated over the surface of the phantom to record the exiting photons. A wide range of optical parameters were respectively assigned to the phantom to mimic both high-absorption and low-absorption media. A reduced scattering albedo (RSA) defined by was used to measure the ratio of photon scatting to absorption. The relative refractive index of the boundary was set at 1.37. An HG phase function with an anisotropic coefficient of 0.9 was assumed for the MC simulation. The MC simulation and the computational scheme of linear Eq. 13 were respectively performed to compute the photon fluence rate and exiting photon flux at detectors for three sets of optical parameters in Table 1 . These results showed that the solutions of the PA equation were in excellent agreement with the results obtained from the MC simulation with a relative error below 3.8%. Here the relative error is defined as . Figures 2a to 4a present a comparison between the PA and MC data for RSAs 3.57, 7.25, and 100, respectively. In contrast, we also performed the numerical experiments with the same phantom setting to compute the photon fluence rates based on the DA model. As a result, the relative errors between the photon fluence rate obtained from the DA model and from the MC simulation were as high as 31.2 and 14.6% for an RSAs of 3.57 and 7.25, respectively, as shown in Figs. 2b and 3b . As we had expected, the DA model had a satisfactory performance with relative errors of 3.8% for an RSA of 100, as shown in Fig. 4b . The numerical experiments demonstrated that the solution obtained from the PA equation is more accurate than the results from the DA model as compared to the MC data over a broad range of optical parameters of biomedical interest. The computational time of the PA equation was about for both fluence rate and flux in our numerical experiments, while that of the MC simulation with the photon number of was about , and that of the DA model in the same setting about . The running time is measured under a 2.8-GHz Intel Xeon CPU with bytes of memory.
Optical parameters used in the numerical experiments.
Discussion and Conclusions
The phase function in the RTE describes the scattering characteristics of the medium and is strictly related to the scattering and absorption coefficients of the medium. The generalized Delta-Eddington phase function is characterized by anisotropy weight and asymmetry factor and shows that the anisotropy weight and asymmetry factor are related to the photon absorbing and scattering in the medium. The proposed concept improves the popular interpretations of the Delta-Eddington and HG phase functions, which are independent of absorbing and scattering in a medium. Based on the generalized Delta-Eddington phase function, the RTE can be significantly reduced to a system of integral equations, which enables the simultaneous computation of both the photon fluence rate and flux vector. MC simulation experiments demonstrated that the solution of the system of integral equations is highly accurate to model photon propagation over a wide range of optical parameters. The system of integral equations is suitable to model the photon propagation in heterogeneous media and can be applied to optical molecular imaging for a small animal with complex anatomical structures.
This work was partially supported by the National Institutes of Health under Grants EB001685, EB006036, EB008476, and CA127189.