## 1.

## Introduction

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 method^{3, 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 approximation^{14, 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) function^{18} and the Delta-Eddington function^{19} 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
$g$
, 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.

## 2.

## Phase Approximation Model

The biological tissue scatters light strongly in the forward direction,^{22} and the scattering phase function can be modeled by a generalized Delta-Eddington function:^{19}

## Eq. 1

$$p(v\cdot {v}^{\prime})=\frac{1}{4\pi}[(1-f)(1+3hv\cdot {v}^{\prime})+2f\delta (1-v\cdot {v}^{\prime})],$$^{2, 19}These parameters are related to a single free parameter $g$ , where $g$ 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 RTE

## Eq. 2

$$v\cdot \nabla L(\mathbf{r},v)+({\mu}_{a}+{\mu}_{s})L(\mathbf{r},v)={\mu}_{s}{\oint}_{{S}^{2}}L(\mathbf{r},{v}^{\prime})p(v,{v}^{\prime})\phantom{\rule{0.2em}{0ex}}\mathrm{d}{v}^{\prime}+\frac{1}{4\pi}Q\left(\mathbf{r}\right),\phantom{\rule{1em}{0ex}}\mathbf{r}\u220a\Omega ,$$## Eq. 3

$$v\cdot \nabla L(\mathbf{r},v)+{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}L(\mathbf{r},v)=\frac{{\stackrel{\u0303}{\mu}}_{s}}{4\pi}[\Phi \left(\mathbf{r}\right)+3hv\cdot \mathbf{J}\left(\mathbf{r}\right)]+\frac{1}{4\pi}Q\left(\mathbf{r}\right),\phantom{\rule{1em}{0ex}}\mathbf{r}\u220a\Omega ,$$## Eq. 4

$$\Phi \left(\mathbf{r}\right)={\oint}_{{S}^{2}}L(\mathbf{r},v)\phantom{\rule{0.2em}{0ex}}\mathrm{d}v\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\mathbf{J}\left(\mathbf{r}\right)={\oint}_{{S}^{2}}vL(\mathbf{r},v)\phantom{\rule{0.2em}{0ex}}\mathrm{d}v.$$^{23}

## Eq. 5

$$L(\mathbf{r},v)=\frac{1}{4\pi}{\int}_{0}^{R}\{{\stackrel{\u0303}{\mu}}_{s}[\Phi (\mathbf{r}-\rho v)-3hv\cdot \mathbf{J}(\mathbf{r}-\rho \mathbf{v})]+Q(\mathbf{r}-\rho v)\}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\phantom{\rule{0.2em}{0ex}}[-{\int}_{0}^{\rho}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(\mathbf{r}-tv)\phantom{\rule{0.2em}{0ex}}\mathrm{d}t]\phantom{\rule{0.2em}{0ex}}\mathrm{d}\rho +L(\mathbf{r}-Rv,v)\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\phantom{\rule{0.2em}{0ex}}[-{\int}_{0}^{R}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(\mathbf{r}-tv)\phantom{\rule{0.2em}{0ex}}\mathrm{d}t],$$## Eq. 6

$$\Phi \left(\mathbf{r}\right)=\frac{1}{4\pi}{\oint}_{{S}^{2}}{\int}_{0}^{R}\{{\stackrel{\u0303}{\mu}}_{s}[\Phi (\mathbf{r}-\rho v)-3hv\cdot \mathbf{J}(\mathbf{r}-\rho v)]+Q(\mathbf{r}-\rho v)\}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\phantom{\rule{0.2em}{0ex}}[-{\int}_{0}^{\rho}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(\mathbf{r}-tv)dt]\mathrm{d}\rho \phantom{\rule{0.2em}{0ex}}\mathrm{d}v+{\oint}_{{S}^{2}}L(\mathbf{r}-Rv,v)\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\phantom{\rule{0.2em}{0ex}}[-{\int}_{0}^{R}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(\mathbf{r}-tv)\phantom{\rule{0.2em}{0ex}}\mathrm{d}t]\phantom{\rule{0.2em}{0ex}}\mathrm{d}v.$$## Eq. 7

$$J\left(r\right)=\frac{1}{4\pi}{\oint}_{{S}^{2}}{\int}_{0}^{R}\{{\stackrel{\u0303}{\mu}}_{s}[\Phi (\mathbf{r}-\rho v)-3hv\cdot \mathbf{J}(r-\rho v)]+Q(\mathbf{r}-\rho v)\}\mathrm{exp}(-{\int}_{0}^{\rho}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(r-tv)\phantom{\rule{0.2em}{0ex}}\mathrm{d}t)v\phantom{\rule{0.2em}{0ex}}\mathrm{d}\rho \phantom{\rule{0.2em}{0ex}}\mathrm{d}v+{\oint}_{{S}^{2}}L(r-Rv,v)\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\phantom{\rule{0.2em}{0ex}}(-{\int}_{0}^{R}{\stackrel{\u0303}{\mu}}_{\mathrm{tr}}(r-tv)\phantom{\rule{0.2em}{0ex}}\mathrm{d}t)v\phantom{\rule{0.2em}{0ex}}\mathrm{d}v.$$## Eq. 8

$$L(\mathbf{r}-Rv,v)={r}_{d}L(\mathbf{r}-Rv,{v}^{\prime}),\phantom{\rule{1em}{0ex}}{v}^{\prime}=v-2(v\cdot \mathbf{n})\mathbf{n},v\cdot \mathbf{n}<0,$$^{14, 16}by ${r}_{d}=-1.4399{\eta}^{-2}+0.7099{\eta}^{-1}+0.6681+0.0636\eta $ , $\eta ={\eta}_{\mathrm{issue}}\u2215{\eta}_{\mathrm{air}}$ . From the first-order approximation of the radiance $L(\mathbf{r}-Rv,{v}^{\prime})$ , the reflected radiance $L(\mathbf{r}-Rv,v)$ on the boundary can be approximated by

## Eq. 9

$$L(\mathbf{r}-Rv,v)\approx \frac{{r}_{d}}{4\pi}\{\Phi (\mathbf{r}-Rv)+3[v-2(v\cdot \mathbf{n})\mathbf{n}]\cdot \mathbf{J}(\mathbf{r}-Rv)\},\phantom{\rule{1em}{0ex}}\mathbf{r}-Rv\u220a\partial \Omega .$$## Eq. 10

$$\{\begin{array}{ccc}\Phi \left(\mathbf{r}\right)\hfill & =\hfill & \frac{1}{4\pi}{\int}_{\Omega}\{{\stackrel{\u0303}{\mu}}_{s}[\Phi \left({\mathbf{r}}^{\prime}\right)-3hv\cdot \mathbf{J}\left({\mathbf{r}}^{\prime}\right)]+Q\left({\mathbf{r}}^{\prime}\right)\}G(\mathbf{r},{\mathbf{r}}^{\prime})\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\mathbf{r}}^{\prime}\hfill \\ \hfill & \hfill & +\frac{{r}_{d}}{4\pi}{\oint}_{\partial \Omega}[\Phi \left({\mathbf{r}}^{\prime}\right)+3{v}^{\prime}\cdot \mathbf{J}\left({\mathbf{r}}^{\prime}\right)]G(\mathrm{r},{\mathbf{r}}^{\prime})(v\cdot \mathbf{n})\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\mathbf{r}}^{\prime}\hfill \\ \mathbf{J}\left(\mathbf{r}\right)\hfill & =\hfill & \frac{1}{4\pi}{\int}_{\Omega}\{{\stackrel{\u0303}{\mu}}_{s}[\Phi \left({\mathbf{r}}^{\prime}\right)-3hv\cdot \mathbf{J}\left({\mathbf{r}}^{\prime}\right)]+Q\left({\mathbf{r}}^{\prime}\right)\}\mathbf{G}(\mathrm{r},{\mathbf{r}}^{\prime})v\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\mathbf{r}}^{\prime}\hfill \\ \hfill & \hfill & +\frac{{r}_{d}}{4\pi}{\oint}_{\partial \Omega}[\Phi \left({\mathbf{r}}^{\prime}\right)+3{v}^{\prime}\cdot \mathbf{J}\left({\mathbf{r}}^{\prime}\right)]G(\mathbf{r},{\mathbf{r}}^{\prime})(v\cdot \mathbf{n})v\phantom{\rule{0.2em}{0ex}}\mathrm{d}{\mathbf{r}}^{\prime},\hfill \end{array}$$^{2}

## Eq. 11

$$\Gamma \left(\mathbf{r}\right)=\mathbf{J}\left(\mathbf{r}\right)\cdot \mathbf{n}\u2215(1-{r}_{d}),\phantom{\rule{1em}{0ex}}\mathrm{for}\phantom{\rule{0.3em}{0ex}}0<{r}_{d}<1,\phantom{\rule{1em}{0ex}}\mathbf{r}\u220a\partial \Omega ,$$^{6}If the absorption coefficient ${\mu}_{a}$ , the scattering coefficient ${\mu}_{s}$ , and conventional anisotropic factor $g$ 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 $f$ and asymmetry factor $h$ 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 $f$ and asymmetry factor $h$ are independent of the phantom geometry.

## 3.

## Numerical Experiments

To perform the numerical computation based on Eq. 10, the region of interest
$\Omega $
can be discretized into finite elements with
$N$
vertex nodes, and the photon fluence rate
$\Phi \left(\mathbf{r}\right)$
and flux
$\mathbf{J}\left(\mathbf{r}\right)$
are approximated in terms of nodal-based basis functions^{24}
${\phi}_{j}\left(\mathbf{r}\right)(j=1,2,\dots ,N)$
,

## Eq. 12

$$\{\begin{array}{c}\Phi \left(\mathbf{r}\right)=\sum _{j=1}^{N}\Phi \left({\mathbf{r}}_{j}\right){\phi}_{j}\left(\mathbf{r}\right)\hfill \\ \mathbf{J}\left(\mathbf{r}\right)=\sum _{j=1}^{N}\mathbf{J}\left({\mathbf{r}}_{j}\right){\phi}_{j}\left(\mathbf{r}\right)\hfill \end{array}.$$## Eq. 13

$$\{\begin{array}{c}\left\{\Phi \right\}={\mathbf{M}}_{1}\left\{\Phi \right\}+{\mathbf{H}}_{1}\left\{\mathbf{J}\right\}+\left\{{\mathbf{Q}}_{1}\right\}\hfill \\ \left\{\mathbf{J}\right\}={\mathbf{M}}_{2}\left\{\Phi \right\}+{\mathbf{H}}_{2}\left\{\mathbf{J}\right\}+\left\{{\mathbf{Q}}_{2}\right\}\hfill \end{array},$$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 $20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ diameter and $20\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ height. We set up a spherical source of $1.0\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ radius with a power of $10\phantom{\rule{0.3em}{0ex}}\mathrm{nw}$ at position of $(-4.0,0.0,10.0)$ 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 ${\mu}_{s}(1-g)\u2215{\mu}_{a}$ was used to measure the ratio of photon scatting to absorption. The relative refractive index $\eta $ 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 ${\sum}_{\mathrm{detectors}}\mid {\Gamma}_{\mathrm{MC}}-{\Gamma}_{\mathrm{PA}}\mid \u2215{\sum}_{\mathrm{detectors}}{\Gamma}_{\mathrm{MC}}$ . 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 $45\phantom{\rule{0.2em}{0ex}}\mathrm{min}$ for both fluence rate and flux in our numerical experiments, while that of the MC simulation with the photon number of ${10}^{7}$ was about $80\phantom{\rule{0.2em}{0ex}}\mathrm{min}$ , and that of the DA model in the same setting about $3\phantom{\rule{0.2em}{0ex}}\mathrm{min}$ . The running time is measured under a 2.8-GHz Intel Xeon CPU with $4\phantom{\rule{0.2em}{0ex}}\mathrm{G}$ bytes of memory.

## Table 1

Optical parameters used in the numerical experiments.

μa(mm−1) | μs(mm−1) | μs(1−g)∕μa | g | η | f | h |
---|---|---|---|---|---|---|

0.35 | 12.50 | 3.57 | 0.9 | 1.37 | 0.938 | 0.90 |

0.20 | 14.50 | 7.25 | 0.9 | 1.37 | 0.950 | 0.90 |

0.008 | 8.0 | 100 | 0.9 | 1.37 | 0.973 | 0.98 |

## 4.

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

## Acknowledgments

This work was partially supported by the National Institutes of Health under Grants EB001685, EB006036, EB008476, and CA127189.