In recent years, the simplified spherical harmonics (${\mathrm{SP}}_{\mathrm{N}}$) equations were applied for modeling light propagation in biological tissues.^1–5 They were considered to effectively remedy the applicability of the diffusion approximation (DA) and reduce the complexity of solving the full spherical harmonics (${\mathrm{P}}_{\mathrm{N}}$) or the discrete ordinates (${\mathrm{S}}_{\mathrm{N}}$) equations. In a previous paper, we derived the analytical solutions to the ${\mathrm{SP}}_{\mathrm{N}}$ equations in turbid medium.³ It will thus be useful to proceed with its extension for studying the response of fluorophores, since fluorescence spectroscopy/imaging has some important applications in the field of biomedical optics.^4–6

In this study, we derive the analytical solutions to the fluorescence simplified spherical harmonics (${\mathrm{FSP}}_{\mathrm{N}}$) equations based on the Green’s function method in the homogeneous turbid medium. The eigen-decomposition strategy developed in our previous paper and the analytical methods to fluorescence diffusion approximation (FDA) proposed by Yalavarthy and Ayyalasomayajula⁶ are used. The analytical method validation is presented against the Monte Carlo (MC) simulations and FDA data.

The principle of the ${\mathrm{SP}}_{\mathrm{N}}$ equations can be extended to the ${\mathrm{FSP}}_{\mathrm{N}}$ equations as a two-stage system of excitation and emission. Here, the ${\mathrm{FSP}}_3$ equations in the steady-state domain are given by the following system of coupled partial differential equations:

where ${\mu }_{ai(v)}={\mu }_{a0(v)}+{\mu }_{s(v)}(1-g^i)$ is the absorption coefficient of order $i$, $g$ is the anisotropy factor, ${\mu }_{s(v)}$ is the scattering coefficient, and subscript $v$ represents excitation $x$ or emission $m$. ${\varphi }_{i(v)}(\mathbf{r})$ is the light fluence rate, and ${\mathrm{S}}_{(v)}(\mathbf{r})$ is the light source. For the excitation process, ${\mathrm{S}}_{(x)}(\mathbf{r})$ denotes the excitation light source. For the emission process, ${\mathrm{S}}_{(m)}(\mathbf{r})={\varphi }_{(x)}(\mathbf{r})N_r(\mathbf{r})$, where ${\varphi }_{(x)}={\varphi }_{1(x)}-(2/3){\varphi }_{2(x)}$ is the total excitation light fluence rate, and $N_r(\mathbf{r})$ is the fluorescence yield distribution, which is the product of the quantum efficiency $\eta$ and the absorption coefficient of the fluorophore ${\mu }_{\mathrm{af}}$. For ${\mathrm{FSP}}_{\mathrm{N}}$ equations, the partially reflective boundary conditions can be used.¹

In the homogeneous medium with a spatially uniform distribution of fluorophores, Eq. (1) can be rewritten as

where ${\mathrm{\Phi }}_{(v)}(\mathbf{r})={[{\varphi }_{1(v)}(\mathbf{r}){\varphi }_{2(v)}(\mathbf{r})]}^T$, ${\epsilon }_{(v)}={[{\epsilon }_{1(v)}{\epsilon }_{2(v)}]}^T$, with the coefficient matrix and the vector ${\epsilon }_{(v)}={[3{\mu }_{a1(v)}-14/3{\mu }_{a3(v)}]}^T$.

By our proposed transformation matrix ${\mathrm{\Phi }}_{(v)}^{\prime }{=\mathbf{B}}_{(v)}{\mathrm{\Phi }}_{(v)}$ and ${\epsilon }_{(v)}^{\prime }={\mathbf{B}}_{(v)}^{-1}{\epsilon }_{(v)}$, where $\mathbf{B}$ is the eigenvector and ${\lambda }_{1(v)}^2,{\lambda }_{2(v)}^2$ are the two real eigenvalues of coefficient matrix $A_{(v)}$,³ the excitation and emission equations can be reduced to the form of a system of independent differential equations that are mathematically consistent with the FDA equation

In Ref. 6, the Green’s functions for the time-domain and frequency-domain FDA have been derived in detail. Here, we can further write the universal formula of the steady-state Green’s function for regular geometries

where $g_{\mathrm{geo}}^{{\psi }_{\mathrm{fl}}}(r-r^{\prime })$ represents the Green’s function of the emission light fluence rate, with the isotropic source $\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$ placed at ${\mathbf{r}}^{\prime }$, and the subscript “geo” represents the relevant geometry under consideration. ${\gamma }_{x,m}^2=c/[3({\mu }_{\mathrm{ax},\mathrm{am}}+{\mu }_{\mathrm{sx},\mathrm{sm}}^{\prime })]$, ${\mu }_{\mathrm{ax},\mathrm{am}}$ and ${\mu }_{\mathrm{sx},\mathrm{sm}}^{\prime }$ are the absorption and reduced scattering coefficients, respectively, and $c$ is the speed of light in the medium. ${\zeta }^2=({\gamma }_m^2-{\gamma }_x^2)/({\gamma }_m^2{\mu }_{\mathrm{ax}}c-{\gamma }_x^2{\mu }_{\mathrm{am}}c)$, $g_{\mathrm{geo}}^{{\psi }_{x,m}}(r-r^{\prime })$ are the Green’s functions evaluated by substituting ${\mu }_{\mathrm{ax}}$ and ${\mu }_{\mathrm{am}}$ into excitation and emission DA equations, respectively.

Based on Eqs. (3) and (4), we are encouraged to further derive the Green’s functions for the ${\mathrm{FSP}}_3$ equations. For the excitation equation, by applying the transformation relationship of ${\mathrm{\Phi }}_{(x)}={\mathbf{B}}_{(x)}{\mathrm{\Phi }}_{(x)}^{\prime }$, the Green’s function of the total excitation light fluence rate is given by

For the emission equation, using the result of Eq. (4) and the principle of linear superposition theory, we can write the Green’s function of ${\varphi }_{i(m)}^{\prime }$ as

where $B_{1j(x,m)}$ and $B_{2j(x,m)}$ denote the entry of ${\mathbf{B}}_{(x,m)}$ at the first and second row and the $j$’th column, respectively, $g_{j(x)}^{\prime }$ and $g_{j(m)}^{\prime }$ are the Green’s functions evaluated by substituting ${\lambda }_{j(x)}$, ${\epsilon }_{j(x)}^{\prime }$, and ${\lambda }_{i(m)}$, ${\epsilon }_{i(m)}^{\prime }$ into Eq. (3).

The Green function of the total emission light fluence rate for the ${\mathrm{FSP}}_3$ equations can be obtained with

The exiting flux at the boundary for the ${\mathrm{FSP}}_3$ equations can be written as¹

where the constants $J_i(i=0,1,2,3)$ are explicitly given in Ref. 1.

From the above derivation, we note that if $g_{i(x,m)}^{\prime }$ are given, the Green’s functions for the ${\mathrm{FSP}}_3$ equations can be readily obtained. As with our previous paper, the analytical solutions for infinite medium and two-dimensional (2-D) circle are given in the following section.

In the case of an infinite homogeneous medium excited by an isotropic point source, the Green’s functions $g_{i(x,m)}^{\prime }$ can be expressed as^3,7

The Green’s functions $g_{i(x,m)}^{\prime }$ of 2-D circular domain which represents an infinite cylindrical medium excited by anaxially parallel and infinite isotropic line source at $r^{\prime }$ can be expressed as^3,7

where $I_n$ and $K_n$ are the modified Bessel functions of the first and second kinds of $n$’th respectively, $r$ is the radius of circle, $a_{n(x,m)}^j$ are unknown coefficients which can be determined by boundary conditions.¹

To demonstrate the effectiveness of the proposed method, the analytical solutions to the ${\mathrm{FSP}}_3$ equations are compared with FDA data and MC simulations. The turbid medium under consideration is homogeneous and the entire medium can act as a source for fluorescent light. In the following comparative investigations, we mainly consider the effect of the different absorption coefficients in biological tissue, and keep the others to be constant: the scattering coefficient ${\mu }_s=10{\mathrm{mm}}^{-1}$, the fluorophore absorption coefficient ${\mu }_{af}=0.003$, the quantum efficiency of the fluorophore $\eta =0.23$, the anisotropy factor $g=0.9$, and the refractive index $n=1.4$.

First, the analytical solution to the ${\mathrm{FSP}}_3$ equations in an infinite turbid medium with relatively weak absorption ${\mu }_{\mathrm{ax}}=0.03{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{am}}=0.02{\mathrm{mm}}^{-1}$ was compared with the FDA calculations and MC simulations. The fluorescence fluence rate versus source–detector distance $r$ (mm) and a plot of the percentage error [$({\varphi }_{\mathrm{fl}-}{\varphi }_{\mathrm{MC}})/{\varphi }_{\mathrm{MC}}\times 100\%$] are shown in Figs. 1(a1) and 1(a2), respectively. As expected, when the source–detector separation is large, both the ${\mathrm{FSP}}_3$ and FDA solutions show an excellent agreement with the MC data. However, when the detector is placed near the isotropic source, the ${\mathrm{FSP}}_3$ shows smaller model error than FDA. Next, a similar comparison was performed for a highly absorbing medium. We increased the absorption coefficient to ${\mu }_{\mathrm{ax}}=2{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{am}}=1{\mathrm{mm}}^{-1}$. We see that the error of the FDA solution is significantly increased; however, the ${\mathrm{FSP}}_3$ solutions are still in a relatively good agreement with the MC in Figs. 1(b1) and 1(b2).

In 2-D circular domain, the radius of circle is $r=5\mathrm{mm}$, the isotropic point light source is placed at ($r-1/{\mu }_{a1}$, 0 deg) and 15 detectors are distributed at equal spacing from ($r$, 22.5 deg) to ($r$, 337.5 deg). First, we compared the exiting fluorescence flux from the ${\mathrm{FSP}}_3$, FDA, and MC calculations with the small absorption coefficients ${\mu }_{\mathrm{ax}}=0.03{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{am}}=0.02{\mathrm{mm}}^{-1}$. Figures 2(a1)and 2(a2) describe the exiting fluorescence flux from the boundary versus the detection angle and a plot of the percentage error ($({\mathrm{\Gamma }}_{\mathrm{fl}-}{\mathrm{\Gamma }}_{\mathrm{MC}})/{\mathrm{\Gamma }}_{\mathrm{MC}}\times 100\%$), respectively. Since the absorption coefficient is small, it is natural to observe that the results obtained from ${\mathrm{FSP}}_3$ and FDA are close to MC data, and the maximum model error is within 6.7%. Then we increased the absorption coefficient to ${\mu }_{\mathrm{ax}}=0.03{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{am}}=0.02{\mathrm{mm}}^{-1}$. As can been seen in Figs. 2(b1) and 2(b2), the error of FDA is more obvious than ${\mathrm{FSP}}_3$, and for the detector at 180 deg the ${\mathrm{FSP}}_3$ and FDA errors reach a maximum of about $-29.5\%$ and $-16.8\%$, respectively.

In conclusion, we have derived the analytical Green’s function for the ${\mathrm{FSP}}_3$ equations for the homogeneous turbid medium and shown the effectiveness of our proposed method. As expected, compared with the FDA, the ${\mathrm{FSP}}_3$ equations show the advantages in the case of the small source–detector separation and high absorption. To the authors’ knowledge, this is the first time that the analytical solutions are presented. Although the analytical Green’s functions derived in this paper are limited to the specific geometries (infinite medium and 2-D circle), the proposed derivation process is principly concise and universally extendable to other regular geometries. Moreover, the method is applicable for an arbitrary order of ${\mathrm{FSP}}_{\mathrm{N}}$ equations.