Radiative transfer equation and diffusion theory for photon transport in biological tissue

In the RTE, six different independent variables define the radiance at any spatial and temporal point (${\displaystyle x}$, ${\displaystyle y}$, and ${\displaystyle z}$ from ${\displaystyle {\vec {r}}}$, polar angle ${\displaystyle \theta }$ and azimuthal angle ${\displaystyle \phi }$ from ${\displaystyle {\hat {s}}}$, and ${\displaystyle t}$). By making appropriate assumptions about the behavior of photons in a scattering medium, the number of independent variables can be reduced. These assumptions lead to the diffusion theory (and diffusion equation) for photon transport. Two assumptions permit the application of diffusion theory to the RTE:

• Relative to scattering events, there are very few absorption events. Likewise, after numerous scattering events, few absorption events will occur and the radiance will become nearly isotropic. This assumption is sometimes called directional broadening.
• In a primarily scattering medium, the time for substantial current density change is much longer than the time to traverse one transport mean free path. Thus, over one transport mean free path, the fractional change in current density is much less than unity. This property is sometimes called temporal broadening.

Both of these assumptions require a high-albedo (predominantly scattering) medium.[1]

Radiance can be expanded on a basis set of spherical harmonics ${\displaystyle Y}$_{n, m}. In diffusion theory, radiance is taken to be largely isotropic, so only the isotropic and first-order anisotropic terms are used: ${\displaystyle L({\vec {r}},{\hat {s}},t)\approx \ \sum _{n=0}^{1}\sum _{m=-n}^{n}L_{n,m}({\vec {r}},t)Y_{n,m}({\hat {s}})}$ where ${\displaystyle L}$_{n, m} are the expansion coefficients. Radiance is expressed with 4 terms; one for n = 0 (the isotropic term) and 3 terms for n = 1 (the anisotropic terms). Using properties of spherical harmonics and the definitions of fluence rate ${\displaystyle \Phi ({\vec {r}},t)}$ and current density ${\displaystyle {\vec {J}}({\vec {r}},t)}$, the isotropic and anisotropic terms can respectively be expressed as follows:

• ${\displaystyle L_{0,0}({\vec {r}},t)Y_{0,0}({\hat {s}})={\frac {\Phi ({\vec {r}},t)}{4\pi }}}$
• ${\displaystyle \sum _{m=-1}^{1}L_{1,m}({\vec {r}},t)Y_{1,m}({\hat {s}})={\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

Hence we can approximate radiance as[1]

${\displaystyle L({\vec {r}},{\hat {s}},t)={\frac {1}{4\pi }}\Phi ({\vec {r}},t)+{\frac {3}{4\pi }}{\vec {J}}({\vec {r}},t)\cdot {\hat {s}}}$

Substituting the above expression for radiance, the RTE can be respectively rewritten in scalar and vector forms as follows (The scattering term of the RTE is integrated over the complete ${\displaystyle 4\pi }$ solid angle. For the vector form, the RTE is multiplied by direction ${\displaystyle {\hat {s}}}$ before evaluation.):[1]

${\displaystyle {\frac {\partial \Phi ({\vec {r}},t)}{c\partial t}}+\mu _{a}\Phi ({\vec {r}},t)+\nabla \cdot {\vec {J}}({\vec {r}},t)=S({\vec {r}},t)}$

${\displaystyle {\frac {\partial {\vec {J}}({\vec {r}},t)}{c\partial t}}+(\mu _{a}+\mu _{s}'){\vec {J}}({\vec {r}},t)+{\frac {1}{3}}\nabla \Phi ({\vec {r}},t)=0}$

The diffusion approximation is limited to systems where reduced scattering coefficients are much larger than their absorption coefficients and having a minimum layer thickness of the order of a few transport mean free path.

Using the second assumption of diffusion theory, we note that the fractional change in current density ${\displaystyle {\vec {J}}({\vec {r}},t)}$ over one transport mean free path is negligible. The vector representation of the diffusion theory RTE reduces to Fick's law ${\displaystyle {\vec {J}}({\vec {r}},t)={\frac {-\nabla \Phi ({\vec {r}},t)}{3(\mu _{a}+\mu _{s}')}}}$, which defines current density in terms of the gradient of fluence rate. Substituting Fick's law into the scalar representation of the RTE gives the diffusion equation:[1]

${\displaystyle {\frac {1}{c}}{\frac {\partial \Phi ({\vec {r}},t)}{\partial t}}+\mu _{a}\Phi ({\vec {r}},t)-\nabla \cdot [D\nabla \Phi ({\vec {r}},t)]=S({\vec {r}},t)}$

${\displaystyle D={\frac {1}{3(\mu _{a}+\mu _{s}')}}}$ is the diffusion coefficient and μ'_s${\displaystyle =(1-g)}$μ_s is the reduced scattering coefficient.
Notably, there is no explicit dependence on the scattering coefficient in the diffusion equation. Instead, only the reduced scattering coefficient appears in the expression for ${\displaystyle D}$. This leads to an important relationship; diffusion is unaffected if the anisotropy of the scattering medium is changed while the reduced scattering coefficient stays constant.[1]

A solution to the diffusion equation for the simple case of a short-pulsed point source in an infinite homogeneous medium is presented in this section. The source term in the diffusion equation becomes ${\displaystyle S({\vec {r}},t,{\vec {r}}',t')=\delta ({\vec {r}}-{\vec {r}}')\delta (t-t')}$, where ${\displaystyle {\vec {r}}}$ is the position at which fluence rate is measured and ${\displaystyle {\vec {r}}'}$ is the position of the source. The pulse peaks at time ${\displaystyle t'}$. The diffusion equation is solved for fluence rate to yield the Green function for the diffusion equation:

${\displaystyle \Phi ({\vec {r}},t;{\vec {r}}',t)={\frac {c}{[4\pi Dc(t-t')]^{3/2}}}\exp \left[-{\frac {\mid {\vec {r}}-{\vec {r}}'\mid ^{2}}{4Dc(t-t')}}\right]\exp[-\mu _{a}c(t-t')]}$

The term ${\displaystyle \exp \left[-\mu _{a}c(t-t')\right]}$ represents the exponential decay in fluence rate due to absorption in accordance with Beer's law. The other terms represent broadening due to scattering. Given the above solution, an arbitrary source can be characterized as a superposition of short-pulsed point sources. Taking time variation out of the diffusion equation gives the following for a time-independent point source ${\displaystyle S({\vec {r}})=\delta ({\vec {r}})}$:

${\displaystyle \Phi ({\vec {r}})={\frac {1}{4\pi Dr}}\exp(-\mu _{\mathrm {eff} }r)}$

${\displaystyle \mu _{\mathrm {eff} }={\sqrt {\frac {\mu _{a}}{D}}}}$ is the effective attenuation coefficient and indicates the rate of spatial decay in fluence.[1]

Consideration of boundary conditions permits use of the diffusion equation to characterize light propagation in media of limited size (where interfaces between the medium and the ambient environment must be considered). To begin to address a boundary, one can consider what happens when photons in the medium reach a boundary (i.e. a surface). The direction-integrated radiance at the boundary and directed into the medium is equal to the direction-integrated radiance at the boundary and directed out of the medium multiplied by reflectance ${\displaystyle R_{F}}$:

${\displaystyle \int _{{\hat {s}}\cdot {\hat {n}}<0}L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega =\int _{{\hat {s}}\cdot {\hat {n}}>0}R_{F}({\hat {s}}\cdot {\hat {n}})L({\vec {r}},{\hat {s}},t){\hat {s}}\cdot {\hat {n}}d\Omega }$

where ${\displaystyle {\hat {n}}}$ is normal to and pointing away from the boundary. The diffusion approximation gives an expression for radiance ${\displaystyle L}$ in terms of fluence rate ${\displaystyle \Phi }$ and current density ${\displaystyle {\vec {J}}}$. Evaluating the above integrals after substitution gives:[3]

${\displaystyle {\frac {\Phi ({\vec {r}},t)}{4}}+{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}=R_{\Phi }{\frac {\Phi ({\vec {r}},t)}{4}}-R_{J}{\vec {J}}({\vec {r}},t)\cdot {\frac {\hat {n}}{2}}}$

• ${\displaystyle R_{\Phi }=\int _{0}^{\pi /2}2\sin \theta \cos \theta R_{F}(\cos \theta )d\theta }$
• ${\displaystyle R_{J}=\int _{0}^{\pi /2}3\sin \theta (\cos \theta )^{2}R_{F}(\cos \theta )d\theta }$

Substituting Fick's law (${\displaystyle {\vec {J}}({\vec {r}},t)=-D\nabla \Phi ({\vec {r}},t)}$) gives, at a distance from the boundary z=0,[3]

${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$

• ${\displaystyle A_{z}=2D{\frac {1+R_{\mathrm {eff} }}{1-R_{\mathrm {eff} }}}}$
• ${\displaystyle R_{\mathrm {eff} }={\frac {R_{\Phi }+R_{J}}{2-R_{\Phi }+R_{J}}}}$

It is desirable to identify a zero-fluence boundary. However, the fluence rate ${\displaystyle \Phi (z=0,t)}$ at a physical boundary is, in general, not zero. An extrapolated boundary, at ${\displaystyle z}$_b for which fluence rate is zero, can be determined to establish image sources. Using a first order Taylor series approximation,

${\displaystyle \left.\Phi (z=-A_{z},t)\approx \Phi (z=0,t)-A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}\right|_{z=0}}$

which evaluates to zero since ${\displaystyle \Phi ({\vec {r}},t)=A_{z}{\frac {\partial \Phi ({\vec {r}},t)}{\partial z}}}$. Thus, by definition, ${\displaystyle z}$_b must be ${\displaystyle -A}$_z as defined above. Notably, when the index of refraction is the same on both sides of the boundary, ${\displaystyle R}$_F is zero and the extrapolated boundary is at ${\displaystyle z}$_b${\displaystyle =-2D}$.[3]

Using boundary conditions, one may approximately characterize diffuse reflectance for a pencil beam normally incident on a semi-infinite medium. The beam will be represented as two point sources in an infinite medium as follows (Figure 2):[1][4]

1. Set scattering anisotropy ${\displaystyle g}$₂${\displaystyle =0}$ for the scattering medium and set the new scattering coefficient μ_s2 to the original μ_s1 multiplied by ${\displaystyle (1-g}$₁${\displaystyle )}$, where ${\displaystyle g}$₁ is the original scattering anisotropy.
2. Convert the pencil beam into an isotropic point source at a depth of one transport mean free path ${\displaystyle l}$' below the surface and power = ${\displaystyle a}$'.
3. Implement the extrapolated boundary condition by adding an image source of opposite sign above the surface at ${\displaystyle l}$'${\displaystyle +2z}$_b.

The two point sources can be characterized as point sources in an infinite medium via

${\displaystyle \Phi _{\infty }(r,\theta ,z;r',\theta ',z')={\frac {1}{4\pi D\rho }}\exp(-\mu _{\mathrm {eff} }\rho )}$

${\displaystyle \rho }$ is the distance from observation point ${\displaystyle (r,\theta ,z)}$ to source location ${\displaystyle (r',\theta ',z')}$ in cylindrical coordinates. The linear combination of the fluence rate contributions from the two image sources is

${\displaystyle \Phi (r,\theta ,z;r',\theta ',z')=a'\Phi _{\infty }(r,\theta ,z;r',\theta ',z')-a'\Phi _{\infty }(r,\theta ,z;r',\theta ',-z'-2z_{b})}$

This can be used to get diffuse reflectance ${\displaystyle R}$_d${\displaystyle (r)}$ via Fick's law:

${\displaystyle \left.R_{d}(r)=D{\frac {\partial \Phi }{\partial z}}\right|_{z=0}={\frac {a'z'(1+\mu _{\mathrm {eff} }\rho _{1})\exp(-\mu _{\mathrm {eff} }\rho _{1})}{4\pi \rho _{1}^{3}}}+{\frac {a'(z'+4D)(1+\mu _{\mathrm {eff} }\rho _{2})\exp(-\mu _{\mathrm {eff} }\rho _{2})}{4\pi \rho _{2}^{3}}}}$

${\displaystyle \rho _{1}}$ is the distance from the observation point ${\displaystyle (r,0,0)}$ to the source at ${\displaystyle (0,0,z')}$ and ${\displaystyle \rho _{2}}$ is the distance from the observation point to the image source at ${\displaystyle (0,0,-z'-2z}$_b${\displaystyle )}$.[1][4]

Let ${\displaystyle \Phi ({\vec {r}},t)}$ be the Green function solution to the diffusion equation for a homogeneous medium of optical properties ${\displaystyle \mu _{a}}$, ${\displaystyle \mu _{s}'}$, then the Green function solution for a homogeneous medium which differs from the former only by optical properties ${\displaystyle {\bar {\mu }}_{a}}$, ${\displaystyle {\bar {\mu }}_{s}'}$, such that ${\displaystyle {\bar {\mu }}_{a}/\mu _{a}={\bar {\mu }}_{s}'/\mu _{s}'}$, can be obtained with the following rescaling:[5]

${\displaystyle {\bar {\Phi }}({\bar {\vec {r}}},{\bar {t}})=\left({\frac {{\bar {\mu }}_{s}'}{\mu _{s}'}}\right)^{3}\Phi ({\vec {r}},t)}$

where ${\displaystyle {\bar {\vec {r}}}={\vec {r}}{\frac {\mu _{s}'}{{\bar {\mu }}_{s}'}}}$ and ${\displaystyle {\bar {t}}=t{\frac {\mu _{s}'}{{\bar {\mu }}_{s}'}}}$.

Such property can also be extended to the radiance in the more general general framework of the RTE, by substituting the transport coefficients ${\displaystyle {\bar {\mu }}_{s}'}$, ${\displaystyle \mu _{s}'}$ with the extinction coefficients ${\displaystyle {\bar {\mu }}_{t}}$, ${\displaystyle \mu _{t}}$.

The usefulness of the property resides in taking the results obtained for a given geometry and set of optical properties, typical of a lab scale setting, rescaling them and extending them to contexts in which it would be complicated to perform measurements due to the sheer extension or inaccessibility.[6]

Let ${\displaystyle \Phi ({\vec {r}},t,\mu _{a}=0)}$ be the Green function solution to the diffusion equation for a non-absorbing homogeneous medium. Then, the Green function solution for the medium when its absorption coefficient is ${\displaystyle \mu _{a}}$ can be obtained as: [5]

${\displaystyle \Phi ({\vec {r}},t,\mu _{a})=\Phi ({\vec {r}},t,\mu _{a}=0)\exp(-\mu _{a}vt)}$

Again, the same property also holds for radiance within the RTE.