Evaluation of a radiative transfer equation and diffusion approximation hybrid forward solver for fluorescence molecular imaging

Dimitris S. Gorpas; Stefan Andersson-Engels

doi:10.1117/1.JBO.17.12.126010

3 December 2012 Evaluation of a radiative transfer equation and diffusion approximation hybrid forward solver for fluorescence molecular imaging

Dimitris S. Gorpas, Stefan Andersson-Engels

Author Affiliations +

Journal of Biomedical Optics, Vol. 17, Issue 12, 126010 (December 2012). https://doi.org/10.1117/1.JBO.17.12.126010

Abstract

The solution of the forward problem in fluorescence molecular imaging strongly influences the successful convergence of the fluorophore reconstruction. The most common approach to meeting this problem has been to apply the diffusion approximation. However, this model is a first-order angular approximation of the radiative transfer equation, and thus is subject to some well-known limitations. This manuscript proposes a methodology that confronts these limitations by applying the radiative transfer equation in spatial regions in which the diffusion approximation gives decreased accuracy. The explicit integro differential equations that formulate this model were solved by applying the Galerkin finite element approximation. The required spatial discretization of the investigated domain was implemented through the Delaunay triangulation, while the azimuthal discretization scheme was used for the angular space. This model has been evaluated on two simulation geometries and the results were compared with results from an independent Monte Carlo method and the radiative transfer equation by calculating the absolute values of the relative errors between these models. The results show that the proposed forward solver can approximate the radiative transfer equation and the Monte Carlo method with better than 95% accuracy, while the accuracy of the diffusion approximation is approximately 10% lower.

1. Introduction

Fluorescence molecular imaging has been extensively studied in recent years as a possible tomographic tool for the geometrical and molecular assessment of regions of diseased tissue, such as malignant tumors embedded into tissues.¹^,² Possible applications include the monitoring and evaluation of therapeutic and/or diagnostic actions of pharmaceutical substances²^,³ and breast cancer diagnosis.¹ Other medical imaging applications that have been proposed include lung inflammation, as well as applications in immunology and cardiology.⁴^–⁶

The above potential applications, together with the interesting properties of fluorescence molecular imaging, make it an attractive technique for integration into specific preclinical, clinical and pharmaceutical applications. Using fluorescence molecular imaging avoids the exposure to ionizing radiation and allows biocompatible fluorescence probes to be used. Further, data with a good signal-to-noise ratio can be obtained.⁷^,⁸ Recent advantages in both computational units and fluorescence probe technology have been important for the development of the field. Modern workstations possess multicore processors and at least 4 GB RAM, permitting the rapid processing of large data volumes. Other advances that have boosted the development of molecular imaging based on fluorescence techniques include the progress in fluorescent probes, where the recent development of the up-converting nanoparticles⁹^–¹¹ is highly promising, and the progress in optical imaging modalities, with extremely sensitive fluorescence acquisition modules and high-quality optics.

However, despite the scientific attention that fluorescence molecular imaging has received in recent years and the aforementioned technological advances, some fundamental challenges have not yet been overcome. One of these challenges is the accurate and robust description of light propagation through tissue, which constitutes the core of what is known as the “forward problem.” The most commonly applied solution to the forward problem currently available is the diffusion approximation.¹²^,¹³ This model, being a first-order angular approximation of the radiative transfer equation,¹⁴^,¹⁵ although very fast in convergence, it also suffers from some important limitations. The output from this model has poor accuracy close to the boundaries and to the light sources, even in the case of scattering dominated media.¹⁴ An alternative method that is often used, the Monte Carlo approximation, requires intense computational power to converge to a solution and thus most of the approaches published are basically two-dimensional solutions of the forward problem.¹⁶^–¹⁸

The radiative transfer equation is the model considered by many researchers as the most accurate solution to the forward problem. The numerical solution methods are the gold standard to solve this model, and the methods most commonly used for spatial discretization are the finite difference, the discrete ordinates, the finite volume, and the finite element methods.¹⁹^–²² The last two methods are generally considered to be more flexible than the others in dealing with multiple boundary conditions or complex geometries, and thus they are most commonly used to solve problems related to light propagation through scattering-dominated media, such as tissue.²² The angular discretization is numerically assessed using discrete ordinates, finite elements, and spherical harmonics.²¹^–²³

This manuscript proposes a methodology for solving the forward problem in fluorescence molecular imaging that successfully confronts the limitations of the diffusion approximation, while requiring less computational power than the Monte Carlo method. It is based on a hybrid model initially proposed as a two-dimensional forward solver in diffuse optical tomography,²³^–²⁵ known as coupled radiative transfer equation and diffusion approximation model. The methodology described in this manuscript is an extension of that model in three dimensions and for the fluorescence molecular imaging application. Since fluorescence molecular imaging requires the model to be applied twice, for the excitation and emission fields, the term “dual coupled” is introduced to describe the proposed model. Furthermore, the method presented in this manuscript uses a simplified phase function,²⁶ instead of the Henyey-Greenstein phase function that is used in most radiative transfer equation solutions presented in literature.¹³^,²¹^,²³^–²⁵

The first three-dimensional solution of the coupled model was presented in a fluorescence molecular imaging application,²⁷ where the forward solver was confronted with the finite element method, based on an unstructured spatial mesh. This nonadaptive discretization scheme introduced significant computational load, especially in the case of radiative transfer equation. Its solution with a single computational unit required coarse spatial and angular discretization schemes that produced oscillating results. In order to solve the forward problem on a single computational unit, either an adaptive discretization scheme should be adopted or a structured mesh should be constructed. The later technique is more preferable in applications involving evaluation of algorithms and thus it was the one adopted in this work.

In this work, the proposed model was evaluated on two nonhomogeneous simulation geometries, and we have investigated in depth its applicability as a forward solver in fluorescence molecular imaging. The first geometry was a cube enclosing ellipsoids that corresponded to different tissues and was used to evaluate the proposed model against an independent Monte Carlo forward solver.²⁸ This work is among the very few that compare the results of a three-dimensional radiative transfer equation forward solver with the results derived from Monte Carlo method²² and it is among the first reports in literature where the three-dimensional coupled radiative transfer equation and diffusion approximation model is verified against an independent Monte Carlo method.

The second geometry was the MOBY digital mouse, developed at the Duke Center for In Vivo Microscopy.²⁹ This simulation geometry has been used in many medical imaging studies, especially computed tomography and single-photon emission computed tomography.³⁰^–³³ A few studies have used this digital mouse in fluorescence molecular imaging.²¹^,³⁴^,³⁵ Using this simulation geometry, various scenarios were assessed, so to compare the proposed model against the widely used diffusion approximation and the more accurate model of radiative transfer equation.

The outline of this paper is as follows. Section 2 briefly presents the dual-coupled radiative transfer equation and diffusion approximation model, and outlines the most important approaches to its solution. The digital mouse is presented, and the optical properties of its various organs are summarized. Section 3 presents the scenarios developed to evaluate the model, while Sec. 4 presents and discusses results obtained with the method. Finally, Sec. 5 discusses future work and the prospects for this forward solver.

2. Methodology

The proposed methodology for solving the forward problem in fluorescence molecular imaging is based on the coupled radiative transfer equation and diffusion approximation model,²⁵^,³⁶ expanded to form a dual-coupled model that can deal with both the excitation field and the emission field. Figure 1 is a block diagram that summarizes the solution of the dual-coupled radiative transfer equation and diffusion approximation model. The use of this model for developing a solution in the frequency domain is briefly described below.

Fig. 1

A block diagram of the dual-coupled radiative transfer equation and diffusion approximation model. The inputs ( $In_0 i$ , $i = 1, \dots, 5$ ) and the outputs ( $Out_0 j$ , $j = 1, 2$ ) are shown.

2.1.

Dual-Coupled Radiative Transfer Equation and Diffusion Approximation

The coupled radiative transfer equation and diffusion approximation model was initially proposed for diffuse optical tomography measurements, and the results presented were mainly two-dimensional.²³^–²⁵ The model was extended to fluorescence molecular imaging in 2010, and three-dimensional simulations were presented,²⁷ while an investigation into applying the model to solve the inverse problem in diffuse optical tomography was carried out in 2011.³⁶ We present here the results of an in-depth study into applying the dual-coupled radiative transfer equation and diffusion approximation model in fluorescence molecular imaging, and present results obtained from a tissue-like simulation phantom.

The analytical expression of this model is as follows:²³

Eq. (1a)

\frac{i \cdot ω}{c} \cdot I_{x / m} (r, \hat{s}) + \hat{s} \cdot \nabla I_{x / m} (r, \hat{s}) + [μ_{α, x / m} (r) + μ_{s, x / m} (r)] \cdot I_{x / m} (r, \hat{s}) - μ_{s, x / m} (r) \cdot \int_{4 π} p_{x / m} (\hat{s}, {\hat{s}}^{'}) \cdot I_{x / m} (r, {\hat{s}}^{'}) \cdot d {\hat{s}}^{'} = Λ_{x / m} (r, \hat{s}) |_{r \in V_{RTE}}

Eq. (1b)

I_{x / m} (r, \hat{s}) = {\begin{cases} 0 & r \in S_{RTE, out} ∖ d S_{src}, \hat{s} \cdot \hat{n} < 0 \\ I_{src} (r, \hat{s}) & r \in d S_{src}, \hat{s} \cdot \hat{n} < 0 \end{cases},

Eq. (1c)

I_{x / m} (r, \hat{s}) = U_{x / m} (r) - \frac{3}{4 \cdot π} \cdot [D_{x / m} (r) \cdot \nabla U_{x / m} (r)] \cdot \hat{s} |_{r \in S_{interface}},

Eq. (1d)

\begin{array}{r} \frac{i \cdot ω}{c} \cdot U_{x / m} (r) - \nabla [D_{x / m} (r) \cdot \nabla U_{x / m} (r)] + μ_{α, x / m} (r) \cdot U_{x / m} (r) \\ = Λ_{0, x / m} (r) \end{array} |_{r \in V_{DA}},

Eq. (1e)

U_{x / m} (r) = - 2 \cdot A \cdot D_{x / m} (r) \cdot \frac{\partial}{\partial \hat{n}} U_{x / m} (r) r \in S_{DA, out}, \hat{s} \cdot \hat{n} < 0,

Eq. (1f)

U_{x / m} (r) = \frac{1}{4 \cdot π} \cdot \int_{4 π} I_{x / m} (r, \hat{s}) \cdot d \hat{s} r \in S_{interface},

where

ω

is the angular modulation frequency of the excitation source in

rad \cdot s^{- 1}

,

c

is the average speed of light in the medium,

μ_{s} (r)

and

μ_{a} (r)

are the scattering and absorption coefficients of the medium, both expressed in

{cm}^{- 1}

,

Λ (r, \hat{s})

is the power radiation from sources inside the medium, and

I (r, \hat{s})

is the radiance at location

r

in direction

\hat{s}

within the medium, measured in

W \cdot {cm}^{- 2} \cdot {sr}^{- 1}

. Furthermore,

\hat{n}

is the surface normal unit vector, pointing outwards,

I_{src} (r, \hat{s})

is the radiance of the excitation source,

U (r)

is the fluence rate or photon density in

W \cdot {cm}^{- 2}

,

D (r)

is the photon diffusion coefficient,³⁷ and

A

is the boundary reflection coefficient, which depends on the mismatch of the refractive indices between the tissue and the air.³⁸ Finally,

p (\hat{s}, {\hat{s}}^{'})

is the scattering phase function.

The well-known Henyey-Greenstein formula³⁹ is the most commonly used phase function in biomedical optics, but we used a phase function that fits better with Mie theory. This phase function is:²⁶^,⁴⁰

Eq. (2)

p (\hat{s}, {\hat{s}}^{'}) = K \cdot {[1 + ε \cdot \cos (\hat{s}, {\hat{s}}^{'})]}^{n_{p}},

where

ε

and

n_{p}

are two independent parameters and

K

is a normalization factor, which for

ε = 1

can be expressed as

Eq. (3)

K_{SAM} = \frac{n_{p} + 1}{2^{n_{p} + 2} \cdot π},

and is known as the “simplified approximate MIE (SAM)” phase function.⁴⁰

Equation (1a) corresponds to the radiative transfer equation, which can be applied in the corresponding subregion, $V_{RTE}$ , of the inspected volume, while Eq. (1b) is the vacuum boundary condition.¹⁴ Similarly, Eq. (1d) corresponds to the diffusion approximation, solved within the related subregion, $V_{DA}$ , and Eq. (1e) is its Robin-type boundary condition.⁴¹^,⁴² Those two light transport models were coupled at the interface $S_{interface}$ through Eqs. (1c) and (1f). The location of the interface boundary was determined by the scattering coefficient of the medium, and it was chosen far enough from the surface to ensure that the diffusion is a valid approximation. Furthermore, it was considered that this interface was located within the diffusion approximation subdomain.²⁵

The excitation case is described by the index $x$ in the equation system of Eq. (1), and the only source that existed in this case was the incident excitation light source $I_{src} (r, \hat{s})$ . All internal light sources had a radiated power of zero, $Λ_{x} (r, \hat{s}) = 0$ . The fluorescence emission is described by the index $m$ , and here, in contrast, no external sources were present, but only the internal fluorescence sources. Thus, $I_{src} (r, \hat{s}) = 0$ and

Eq. (4)

Λ_{m} (r, \hat{s}) = Λ_{0, m} (r) = \frac{η \cdot μ_{α, x}^{fluo} (r)}{1 + i \cdot ω \cdot τ (r)} \cdot U_{x} (r),

where

η

and

τ (r)

are the quantum yield and lifetime of the fluorophore, respectively. The index fluo describes the optical properties of the fluorophore. The model described by Eq. (1) was based on the assumption that the presence of the fluorophore did not alter the scattering coefficient of the medium. The absorption coefficient varied with the fluorophore concentration and was modeled as the sum of the absorption coefficients of the tissue and of the fluorophore.

2.1.1.

Numerical solution of the forward problem

We have used the finite element method to solve the equations of the proposed dual-coupled radiative transfer equation and diffusion approximation model. The integro-differential equation system presented in Eq. (1) was transformed to a weak variational formalism by applying the Galerkin finite element method⁴³ and using the streamline diffusion modification⁴⁴ in the radiative transfer equation subdomain, in order to minimize the convergence oscillating results. The final linear algebraic system, in its matrix formalism, is as follows:²³

Eq. (5)

[\begin{matrix} A_{RTE, x / m} & B_{RTE, x / m} \\ B_{DA, x / m} & A_{DA, x / m} \end{matrix}] \cdot [\begin{matrix} α_{x / m} \\ a_{x / m} \end{matrix}] = [\begin{matrix} C_{RTE, x / m} \\ C_{DA, x / m} \end{matrix}],

where

α_{x / m} = [α_{i l, x / m}]

is the radiance at the nodal points

i

and in directions

l

of the radiation transfer equation subdomain, and

a_{x / m} = [a_{k, x / m}]

is the fluence rate at the nodal points

k

of the diffusion approximation subdomain. Moreover, the matrix

A_{RTE, x / m}

contains the finite element approximation of the radiative transfer equation, and is of the form

Eq. (6)

A_{RTE, x / m} (h_{1}, h_{2}) = A_{0}^{x / m} (h_{1}, h_{2}) + A_{1}^{x / m} (h_{1}, h_{2}) + A_{2}^{x / m} (h_{1}, h_{2}) + A_{3}^{x / m} (h_{1}, h_{2}) + A_{4}^{x / m} (h_{1}, h_{2}),

where

h_{1} = N_{a} \cdot (j - 1) + q

,

h_{2} = N_{a} \cdot (i - 1) + l

,

i = 1, ..., N_{n}

,

l = 1, ..., N_{a}

,

j = 1, ..., N_{n}

and

q = 1, ..., N_{a}

, where

N_{n}

is the number of spatial nodes and

N_{a}

is the number of angular directions in the radiative transfer equation subdomain. The matrices on the right side of Eq. (1) are given by

Eq. (7a)

A_{0}^{x / m} (h_{1}, h_{2}) = \frac{i \cdot ω}{c} \cdot \int_{V_{RTE}} ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r \cdot \int_{4 π} ψ_{l} (\hat{s}) \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} + \frac{i \cdot ω}{c} \cdot \int_{V_{RTE}} \int_{4 π} δ_{x / m} (r) \cdot \hat{s} \cdot \nabla ψ_{j} (r) \cdot ψ_{q} (\hat{s}) \cdot ψ_{l} (\hat{s}) \cdot d \hat{s} \cdot ψ_{i} (r) \cdot d r,

Eq. (7b)

A_{1}^{x / m} (h_{1}, h_{2}) = - \int_{V_{RTE}} \int_{4 π} \hat{s} \cdot \nabla ψ_{j} (r) \cdot ψ_{q} (\hat{s}) \cdot ψ_{l} (\hat{s}) \cdot d \hat{s} \cdot ψ_{i} (r) \cdot d r + \int_{V_{RTE}} \int_{4 π} δ_{x / m} (r) \cdot \hat{s} \cdot \nabla ψ_{j} (r) \cdot ψ_{q} (\hat{s}) \cdot \hat{s} \cdot \nabla ψ_{i} (r) \cdot ψ_{l} (\hat{s}) \cdot d \hat{s} \cdot d r,

Eq. (7c)

A_{2}^{x / m} (h_{1}, h_{2}) = \int_{S_{RTE}} ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r_{S_{RTE}} \cdot \int_{4 π} {(\hat{s} \cdot \hat{n})}_{+} \cdot ψ_{l} (\hat{s}) \cdot ψ_{q} (\hat{s}) \cdot d \hat{s}

Eq. (7d)

A_{3}^{x / m} (h_{1}, h_{2}) = \int_{V_{RTE}} [μ_{α, x / m} (r) + μ_{s, x / m} (r)] \cdot ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r \cdot \int_{4 π} ψ_{l} (\hat{s}) \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} + \int_{V_{RTE}} \int_{4 π} {\begin{cases} δ_{x / m} (r) \cdot [μ_{a, x / m} (r) + μ_{s, x / m} (r)] \\ \cdot \hat{s} \cdot \nabla ψ_{j} (r) \cdot ψ_{q} (\hat{s}) \cdot ψ_{l} (\hat{s}) \cdot d \hat{s} \cdot ψ_{i} (r) \cdot d r \end{cases}},

Eq. (7e)

A_{4}^{x / m} (h_{1}, h_{2}) = - \int_{V_{RTE}} μ_{s, x / m} (r) \cdot ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r \cdot \int_{4 π} \int_{4 π} p_{x / m} (\hat{s}, {\hat{s}}^{'}) \cdot ψ_{l} ({\hat{s}}^{'}) \cdot d {\hat{s}}^{'} \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} - \int_{V_{RTE}} \int_{4 π} {\begin{cases} δ_{x / m} (r) \cdot [μ_{a, x / m} (r) + μ_{s, x / m} (r)] \cdot \hat{s} \cdot \nabla ψ_{j} (r) \cdot ψ_{q} (\hat{s}) \\ \cdot \int_{4 π} p_{x / m} (\hat{s}, {\hat{s}}^{'}) \cdot ψ_{l} ({\hat{s}}^{'}) \cdot d {\hat{s}}^{'} \cdot d \hat{s} \cdot ψ_{i} (r) \cdot d r \end{cases}} .

In Eq. (7),

ψ_{i} (r)

and

ψ_{l} (\hat{s})

are the nodal basis functions of the spatial and angular finite element discretization of the radiative transfer equation subdomain, while

ψ_{j} (r)

and

ψ_{q} (\hat{s})

are the spatial and angular basis functions of the radiative transfer equation piecewise linear functions. The Galerkin methodology is prone to the formation of unstable systems and oscillating results. In order to avoid this, the radiative transfer equation test function,

ψ (r, \hat{s})

, was expanded using the following formalism:

Eq. (8)

ψ (r, \hat{s}) : = ψ (r, \hat{s}) + δ (r) \cdot \hat{s} \cdot ψ (r, \hat{s}),

which expresses the streamline diffusion modification.⁴⁴ The smoothing parameter,

δ (r)

, is a spatial nodal-wise varying constant that depends on the local absorption and scattering coefficients. It can be estimated from:

Eq. (9)

δ (r) = \min {\frac{1}{μ_{α, x / m} (r) + μ_{s, x / m} (r)}, \frac{25}{(μ_{α, x / m} (r) + μ_{s, x / m} (r)) \cdot r_{j}}},

where

r_{j}

is the distance between spatial node

j

and the excitation source location.²⁷

The matrix of diffusion approximation finite elements is $A_{DA, x / m}$ in the system of Eq. (1). This matrix is given by:

Eq. (10)

A_{DA, x / m} (p, k) = {Ke}_{x / m} (p, k) + {Me}_{x / m} (p, k) + {Pe}_{x / m} (p, k),

where

p = 1, ..., N

and

k = 1, ..., N

, with

N

being the number of nodal points in the diffusion approximation subdomain. Furthermore, the mass matrix,⁴³

K e_{x / m} (p, k)

, the stiffness matrix,⁴³

M e_{x / m} (p, k)

, and the boundary condition matrix,

P e_{x / m} (p, k)

, of the finite element approximation are given by

Eq. (11a)

{Ke}_{x / m} (p, k) = \frac{i \cdot ω}{c} \cdot \int_{V_{DA}} y_{k} (r) \cdot y_{p} (r) \cdot d r + \int_{V_{DA}} μ_{α, x / m} (r) \cdot y_{k} (r) \cdot y_{p} (r) \cdot d r,

Eq. (11b)

{Me}_{x / m} (p, k) = \int_{V_{DA}} D_{x / m} (r) \cdot \nabla y_{k} (r) \cdot \nabla y_{p} (r) \cdot d r,

Eq. (11c)

{Pe}_{x / m} (p, k) = \frac{1}{2 \cdot A} \cdot \int_{S_{DA}} y_{k} (r) \cdot y_{p} (r) \cdot d r_{S_{DA}},

where

y_{k} (r)

are the nodal basis functions of the finite element discretization of the diffusion approximation subdomain, and

y_{p} (r)

are the spatial basis functions of the diffusion approximation piecewise linear functions.

The matrices $B_{RTE, x / m}$ and $B_{DA, x / m}$ of Eq. (1) contain the finite elements of the interface boundary conditions and are given by

Eq. (12a)

B_{RTE, x / m} (h_{1}, k) = - \int_{S_{interface}} y_{k} (r) \cdot ψ_{j} (r) \cdot d r_{S_{interface}} \cdot \int_{4 π} {(\hat{s} \cdot n)}_{-} \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} + \frac{3}{4 \cdot π} \cdot \int_{S_{interface}} D_{x / m} (r) \cdot \int_{4 π} {(\hat{s} \cdot \hat{n})}_{-} \cdot [\nabla y_{k} (r) \cdot \hat{s}] \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} \cdot ψ_{j} (r) \cdot d r_{S_{interface}},

Eq. (12b)

B_{DA, x / m} (p, h_{2}) = - \frac{1}{4 \cdot π} \cdot \int_{S_{interface}} D_{x / m} (r) \cdot [\hat{n} \cdot \nabla ψ_{i} (r)] \cdot y_{p} (r) \cdot d r_{S_{interface}} \cdot \int_{4 π} ψ_{l} (\hat{s}) \cdot d \hat{s} .

The matrix on the right of Eq. (5) is the finite element matrix of the sources, and it can be rewritten in the form

Eq. (13)

[\begin{matrix} C_{RTE, x / m} \\ C_{DA, x / m} \end{matrix}] = [\begin{matrix} c_{RTE, x / m} \cdot α \\ c_{DA, x / m} \cdot a \end{matrix}] .

The only light source in the excitation field was the incident source,

I_{src} (r, \hat{s})

, which was also located only in the radiative transfer equation subdomain. Thus, the diffusion approximation part of Eq. (13) equals zero, while the finite element matrix

c_{RTE, x}

is

Eq. (14)

c_{RTE, x} (h_{1}, h_{2}) = \int_{S_{RTE}} ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r_{S_{RTE}} \cdot \int_{4 π} {(\hat{s} \cdot \hat{n})}_{-} \cdot ψ_{l} (\hat{s}) \cdot ψ_{q} (\hat{s}) \cdot d \hat{s} .

Further,

α = α_{src} = [α_{i l, src}]

is the intensity of the excitation source at the surface nodal points

i

in the directions

l

.

The finite element matrices $c_{RTE, m}$ and $c_{DA, m}$ that describe the emission field are

Eq. (15a)

c_{RTE, m} (h_{1}, h_{2}) = \frac{1}{4 \cdot π} \cdot \int_{V_{RTE}} \frac{η \cdot μ_{α, x}^{fluo} (r)}{1 + i \cdot ω \cdot τ (r)} \cdot ψ_{i} (r) \cdot ψ_{j} (r) \cdot d r \cdot \int_{4 π} ψ_{l} (\hat{s}) d \hat{s} \cdot \int_{4 π} ψ_{q} (\hat{s}) \cdot d \hat{s},

Eq. (15b)

c_{DA, m} (p, k) = \int_{V_{DA}} \frac{η \cdot μ_{α, x}^{fluo} (r)}{1 + i \cdot ω \cdot τ (r)} \cdot y_{k} (r) \cdot y_{p} (r) \cdot d r,

while

α = α_{x}

and

a = a_{x}

are the radiance and fluence rate calculated from the dual-coupled radiative transfer equation and diffusion approximation model for the excitation case. The importance of applying the same finite element approximation and region discretization between the excitation and emission cases becomes apparent at this point. If these were not applied, the outcome of the excitation field should be translated to the emission approximation, which would increase the complexity of the method and decrease its time efficiency.

The block diagram of Fig. 1 shows that the assembly of the diffusion approximation finite element matrix was implemented in one step, while the radiative transfer equation was divided into spatial and angular finite element approximations. Thus, after the solution of all the spatial and angular nodal integrals, the radiative transfer equation finite element matrix was assembled by applying the Kronecker product.³⁵

2.2.

Digital Tissue Geometry and Region Discretization

Most forward problem solvers used for fluorescence molecular imaging and for diffuse optical tomography are evaluated using simulation geometrics with homogeneous background optical properties.¹³^,²⁷^,³⁶ Furthermore, only a few forward solvers are based on the radiative transfer equation, and only a few reports have presented three-dimensional results.²²^,²⁷^,⁴⁵ We have studied a forward solver based on the radiative transfer equation, and examined its application to imaging a three-dimensional tissue-like simulation geometry. We have compared results from the proposed model with those obtained from the three-dimensional radiative transfer equation solution from the same digital phantom.

The virtual phantom we have used was based on the MOBY mouse.²⁹ The various organs of this digital mouse phantom have been implemented using nonuniform rational b-spline surfaces. High-resolution three-dimensional magnetic resonance microscopy data, obtained from the Duke Center for In Vivo Microscopy, was used as the basis for the formulation of the surfaces.²⁹

This virtual animal model was discretized using the open-source Amide medical image data analysis tool with a step size of $d = 0.05 cm$ on the $x$ , $y$ and $z$ axes. This animal model is presented in Fig. 2, where both volume visualization and the discretized version are shown, with the superficial tissue excluded to reveal the various organs.

Fig. 2

The MOBY phantom as a volume image (a) and under the discretized scheme (b) with the body tissue excluded, allowing the inner organs to be seen clearly.

Each organ of the digital mouse was labeled with a unique number, which corresponds to a color label in Fig. 2. These numbers were assigned during the discretization process and were used to assign optical properties to the organs of the MOBY phantom. Table 1 summarizes the optical properties in the near infrared (NIR) spectral region, in the wavelength range 600 to 800 nm, obtained from Refs. 13, 21, 46, and 47.

Table 1

Optical properties of the MOBY phantom.

Organ	Absorption coefficient $μ_{a}$ ( ${cm}^{- 1}$ )	Scattering coefficient $μ_{s}$ ( ${cm}^{- 1}$ )	Mean cosine $g$
Body (muscle)	0.027	118	0.97
Skeleton	0.08	69	0.80
Heart (tissue)	0.35	167	0.98
Liver	6.50	144	0.95
Lung	8.40	36	0.95
Stomach	1.20	200	0.90
Pancreas	1.20	200	0.90
Kidney	0.01	73	0.85
Spleen	2.80	59	0.78
Small intestine	1.20	200	0.90
Large intestine	1.20	200	0.90
Bladder	1.25	51	0.95
Testicles	0.35	394	0.69
Brain (white matter)	1.58	51	0.96
Brain (gray matter)	2.63	60	0.88
Heart (blood)	1.30	1246	0.99

All optical properties were considered to be independent of wavelength in order to simplify the simulation. Furthermore, the mean cosine was considered to be independent of the tissue type and equal to $g ≅ 0.9$ .

The digital phantom was initially simulated to be in a chamber filled with an Intralipid-ink solution that matched the bulk optical properties of the body tissue. It was possible in this way to study the excitation and emission, avoiding complications arising from partial refraction and reflection at the animal boundary. This is a widely adopted technique even for in vivo small animal fluorescence molecular imaging and for diffuse optical tomography.⁴^,²¹

Although this simulation geometry can be adequately addressed by diffusion approximation, the scope of these simulation experiments was to validate that the proposed model outperforms diffusion approximation even in situations where the later model is expected to provide relatively accurate results. However, the advantage gained by the application of the radiative transfer equation was demonstrated by replacing the Intralipid-Ink solution with a low scattering index matched interface. The simulated absorption coefficient of that interface was equal to $μ_{a} = 0.02 {cm}^{- 1}$ and the reduced scattering coefficient equal to $μ_{s}^{'} = 0.05 {cm}^{- 1}$ and they were wavelength independent.

The simulation experiments were conducted in a cubic region of dimension $2 \times 2 \times 2 {cm}^{3}$ , discretized with a step size of $d = 0.05 cm$ . This produced a three-dimensional grid of 68,921 nodes. The Delaunay triangulation technique was used to formulate structured tetrahedral elements, defined by

Eq. (16)

[\begin{matrix} ζ_{1} \\ ζ_{2} \\ ζ_{3} \\ ζ_{4} \end{matrix}] = \frac{1}{6 \cdot V} \cdot [\begin{matrix} a_{1} & b_{1} & c_{1} & d_{1} \\ a_{2} & b_{2} & c_{2} & d_{2} \\ a_{3} & b_{3} & c_{3} & d_{3} \\ a_{4} & b_{4} & c_{4} & d_{4} \end{matrix}] \cdot [\begin{matrix} 1 \\ x \\ y \\ z \end{matrix}],

where

ζ_{i}

are the isoparametric or simplex coordinates,

(x, y, z)

are the global coordinates of the nodal points,

V

is the volume of the element, and the 16 elements of the matrix

[\begin{matrix} a_{i} & b_{i} & c_{i} & d_{i} \end{matrix}]

depend only on the vertex locations of each tetrahedron.⁴³ The number of tetrahedral elements produced using this spatial discretization was 384,000.

The spatial finite element matrices of both the radiative transfer equation and the diffusion approximation subdomains were estimated from Eq. (16) by transforming all the spatial integrals into simplex coordinates. The integrals thus obtained were solved using⁴³

Eq. (17)

\int_{V_{e l}} ζ_{1}^{i} \cdot ζ_{2}^{j} \cdot ζ_{3}^{k} \cdot ζ_{4}^{l} \cdot d r = \frac{i! \cdot j! \cdot k! \cdot l!}{(i + j + k + l + 3)!} \cdot 6 \cdot V .

All the spatial integrals in Eq. (5) were assessed using the same method as that used for these basic transformations.

The azimuthal technique was used for the angular discretization of the radiative transfer equation subregion, for an angular discretization of $N_{θ} \times N_{φ} = 6 \times 4$ angles, where $N_{θ}$ is the number of control angles over the polar angle $θ$ and $N_{φ}$ is the number of control angles over the azimuthal angle $φ$ . The number of angular nodes considered for each spatial node was 28, using this angular discretization process. Moreover, each octant of the angular space $4 π$ , at any spatial nodal point, was discretized using constant latitude and longitude into $N_{θ} \times N_{φ} = 24$ elementary solid angles $Δ s^{m n}$ as described by Grissa et al.:⁴⁸

Eq. (18)

Δ {\hat{s}}^{m n} = \int_{φ^{-}}^{φ^{+}} \int_{θ^{-}}^{θ^{+}} \sin θ \cdot d θ \cdot d φ .

The basis functions were defined piecewise, element by element, and were taken as a bilinear function over every control angle element. The functions that describe this restriction are known as “shape functions” and they are⁴⁹

Eq. (19a)

ψ_{1}^{m n} (θ, φ) = \frac{θ - θ^{-}}{θ^{+} - θ^{-}} \cdot \frac{φ - φ^{-}}{φ^{+} - φ^{-}},

Eq. (19b)

ψ_{2}^{m n} (θ, φ) = \frac{θ - θ^{+}}{θ^{-} - θ^{+}} \cdot \frac{φ - φ^{-}}{φ^{+} - φ^{-}},

Eq. (19c)

ψ_{3}^{m n} (θ, φ) = \frac{θ - θ^{+}}{θ^{-} - θ^{+}} \cdot \frac{φ - φ^{+}}{φ^{-} - φ^{+}}

Eq. (19d)

ψ_{4}^{m n} (θ, φ) = \frac{θ - θ^{-}}{θ^{+} - θ^{-}} \cdot \frac{φ - φ^{+}}{φ^{-} - φ^{+}} .

These shape functions were used during the solution of the angular integrals of the linear system of Eq. (5). The solution of the phase function defined in Eq. (3) was obtained by first taking the binomial expansion

Eq. (20)

p (\hat{s}, {\hat{s}}^{'}) = K \cdot {[1 + ε \cdot \cos (\hat{s}, {\hat{s}}^{'})]}^{n_{p}} = K \cdot {[1 + ε \cdot \sin θ \cdot \sin θ^{'} \cdot \cos (φ - φ^{'}) + ε \cdot \cos θ \cdot \cos θ^{'}]}^{n_{p}} = K \cdot \sum_{m = 0}^{n_{p}} {\begin{cases} \frac{n_{p}!}{m! \cdot (n_{p} - m)!} \cdot \sum_{j = 0}^{n_{p} - m} [\frac{(n_{p} - m)!}{j! \cdot (n_{p} - m - j)!} \cdot {(ε \cdot \cos θ \cdot \cos θ^{'})}^{j}] \\ \cdot {[ε \cdot \sin θ \cdot \sin θ^{'} \cdot \cos (φ - φ^{'})]}^{m} \end{cases}} .

The solution was subsequently based on the shape functions of Eq. (19). All angular integrals were solved using the same method as that used to solve the spatial integrals, while analytical solutions were obtained whenever required.⁴³^,⁵⁰

2.3.

Solution to the Forward Problem

The discretization procedure described above was applied to the MOBY digital mouse, and the dual-coupled radiative transfer equation and diffusion approximation model was subsequently solved. Figure 3 shows the geometry defined for the digital mouse. The interface between the two subregions of the coupled model was located 0.3 cm below the top surface of the cubic region. The interface divided the cubic region into the required two subdomains. Of these, the radiative transfer equation subdomain consisted of 11,767 nodal points (57,600 tetrahedral elements) and 28 angular nodal points (24 angular elements); while the diffusion approximation subdomain consisted of 58,835 nodal points (326,400 tetrahedral elements).

Fig. 3

The geometry of the digital mouse model. The tetrahedral element and the angular nodes are shown in the upper right corner, along with the counter-clockwise labeling of the elementary nodes.

The finite element matrices of Eqs. (6) and (15) were assembled under the premises described above, and the resulting linear matrix system of Eq. (5) was solved using the stabilized biconjugate gradient (BiCGStab) method.⁵¹ This method provides smooth and fast convergence of nonsymmetric large-scale problems, and is admirably suited to solving the system described by Eq. (5). The method is often used in situations described by radiative transfer equations, since it does not require intense computational resources. This allows denser discretization schemes to be used than those used with other iterative methods. The method works reasonably well for radiative transfer equation problems⁵² and we have therefore not considered it necessary to test other solvers. The tolerance of the relative residual, which defines the convergence of the method, was set to $10^{- 10}$ .

Finally, all the results presented in the following sections were estimated using the equations:

Eq. (21a)

M_{A} (r) = | U (r ∶ ω) | = \sqrt{Re {[U (r ∶ ω)]}^{2} + Im {[U (r ∶ ω)]}^{2}},

Eq. (21b)

M_{θ} (r) = \arg [U (r ∶ ω)] = a \tan 2 {Im [U (r ∶ ω)], Re [U (r ∶ ω)]},

where

M_{A}

is the fluence rate modulation magnitude and

M_{θ}

is its phase shift. The exitance from radiative transfer was estimated after using Eq. (1f) to estimate the fluence rate from the radiance values.

All simulations presented in this manuscript were performed on a 2.67 GHz Intel Core i7 computational unit with 20 GB physical memory.

3. Evaluation Process and Simulation Scenarios

The evaluation of the dual-coupled radiative transfer equation and diffusion approximation model was implemented through two stages. Initially, this model was compared with an independent Monte Carlo forward solver. The diffusion approximation and the full radiative transfer equation models were also compared with the Monte Carlo forward solver. The MOSE software²⁸ was used for that purpose. The simulation geometry was based on standard building blocks, which approximate the shapes of biological tissues, and was identical for all four forward solvers.

More specifically, it was a cubic region with bulk tissue optical properties, inside of which there were five subregions; two ellipsoids with optical properties corresponding to kidneys, an ellipsoid with optical properties corresponding to lung, an ellipsoid with optical properties corresponding to stomach and an ellipsoid with optical properties corresponding to skeleton. The optical properties of these regions were assigned according to Table 1.

A sphere located at depth $d z = - 0.65 cm$ from the top surface of this geometry and with radius $r = 0.15 cm$ simulated the fluorophores inclusion. The optical properties of the simulated fluorophore matched those of Alexa Fluor 680. Thus, the absorption coefficient was $μ_{a, x}^{fluo} = 0.42 {cm}^{- 1}$ at the excitation wavelength, $λ_{x} = 680 nm$ , and $μ_{a, m}^{fluo} = 0.25 {cm}^{- 1}$ at the emission wavelength, $λ_{m} = 700 nm$ . Furthermore, the quantum efficiency was set to $η = 0.36$ , while the lifetime of the fluorophore was set to $τ (r) = 1.2 \cdot 10^{- 9} s$ . The lifetime of the fluorophore was independent of its position in the simulation geometry. Figure 4 presents this simulation geometry.

Fig. 4

The geometry used for the evaluation of the radiative transfer equation, the diffusion approximation and the hybrid models against the Monte Carlo forward solver.

Next we evaluated the dual-coupled radiative transfer equation and diffusion approximation model for solving the forward problem in fluorescence molecular imaging by comparing the results with those obtained from the widely applied diffusion approximation. We also compared the results with those from the radiative transfer equation. More specifically, the forward problem for the MOBY phantom was solved using the three models: the diffusion approximation, the radiative transfer equation, and the dual coupled. The spatial discretization used was the same for all models, as was also the angular discretization, where relevant. Table 2 presents the qualitative quantities of the finite element approximations for the three models.

Table 2.

Spatial and angular discretization schemes for the diffusion approximation (DA) the radiative transfer equation (RTE) and the proposed dual-coupled model (RTE-DA), and the size of the corresponding assembly matrices.

	DA	RTE	RTE-DA
Tetrahedral elements	384,000	384,000	57,600–RTE 326,400–DA
Nodal points	68,921	68,921	11,767–RTE 58,835–DA
Angular elements	—	24	24–RTE
Angular nodes	—	28	28 – RTE
Assembly matrix size	[ $68,921 \times 68,921$ ]	[ $1,929,788 \times 1,929,788$ ]	[ $388,311 \times 388,311$ ]

The dual-coupled radiative transfer equation and diffusion approximation and the diffusion approximation models were compared with the radiative transfer equation by calculating

Eq. (22)

{ARE}_{A / θ, DA / RTE - DA} = | \frac{M_{A / θ, RTE} - M_{A / θ, DA / RTE - DA}}{M_{A / θ, RTE}} | \times 100 %,

where ARE is the absolute values of the relative error between the magnitude of the fluence rate or the phase shift derived from the radiative transfer equation solution,

M_{A / θ, RTE}

, and the magnitude of the fluence rate or the phase shift obtained from the diffusion approximation,

M_{A / θ, DA}

, or from the dual-coupled model,

M_{A / θ, RTE - DA}

.

Results obtained by considering fluorophore inclusions within the digital mouse were compared. The optical properties of the simulated fluorophore matched those of Alexa Fluor 680. The simulation scenarios used to evaluate the proposed model considered the fluorophore inclusions to be located at four positions within the MOBY phantom. More specifically, a sphere with radius $R = 0.15 cm$ was initially placed at depth $d z = - 0.2 cm$ from the simulation geometry top surface, Fig. 5(a). The second scenario was that an identical sphere was located deeper within the simulation geometry, at depth $d z = - 0.5 cm$ , Fig. 5(b). Scenario 3 considered the mouse kidneys to be labeled with the Alexa Fluor 680 probe, Fig. 5(c), while Scenario 4 considered the bladder to be labeled, Fig. 5(d). Finally, the last two scenarios were identical with Scenarios 1 and 3 correspondingly, with the difference that the Intralipid solution was replaced with the low scattering index matched interface. Since diffusion approximation was not valid for these two last scenarios, only the dual coupled and the radiative transfer equation models addressed them. Moreover, the coupling interface of the proposed model, for these two scenarios, was translated 1 cm deeper inside the animal model, to ensure that the diffusion approximation of the hybrid model was valid at all nodes of the corresponding domain.

Fig. 5

The locations of the fluorophores for the simulation scenarios. A single inclusion at depth $d z = - 0.2 cm$ (a) or $d z = - 0.5 cm$ (b) from the top surface of the MOBY digital mouse. The third scenario (c) was that the mouse kidneys were labeled, while the fourth (d) was that the bladder was labeled.

All simulation scenarios included an epi-illumination acquisition system, and thus the top surface of the cubic region was the area inspected. The excitation source was incident on the same plane of the region, and it was a point source with a Dirac intensity profile and an angular modulation frequency $ω = 100 MHz$ , directed inwards along the z-axis. For the diffusion approximation model, this source was simulated as an isotropic point source located one transport length ( $1 / μ_{s}^{'} ≅ 0.05 cm$ ) below the region top surface.

4. Results and Discussion

We now present and discuss solutions of the forward problem in fluorescence molecular imaging for the six simulation scenarios described above, as well as for the comparison with the independent Monte Carlo forward solver.

4.1.

Evaluation Against Monte Carlo

Initially we evaluated the dual-coupled radiative transfer equation and diffusion approximation model against an independent Monte Carlo forward solver. The logarithms of the fluence rate modulation amplitudes (normalized with respect to their maximum values) of the exitance at the region top surface for the diffusion approximation, the dual-coupled model, the radiative transfer equation model, and the Monte Carlo model are shown in Fig. 6. In this figure are also shown the corresponding phase shifts for the excitation and emission fields.

Fig. 6

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), the radiative transfer equation (c), and the Monte Carlo (d) models. The logarithm of the normalized fluence rate modulation amplitude, on the simulation geometry top surface, of the excitation source is presented in the plots indexed (1.1), and for the emission field in those indexed (2.1), while the corresponding phase shift, in degrees, is presented in those indexed (1.2) and (2.2).

The simulation geometry optical properties, for all four models, were assigned according to Fig. 4 and Table 1. For the Monte Carlo-based simulation, the number of the launched photons was $10^{9}$ , which is considered a sufficient number for accurate three-dimensional results.¹³^,²⁸^,⁵³ Higher photon numbers would increase the accuracy of the Monte Carlo method, especially for the emission field and for the phase shift in particular. However, we compared Monte Carlo results for $10^{7}$ , $5 \cdot 10^{7}$ , $10^{8}$ , $5 \cdot 10^{8}$ , and $10^{9}$ photons and the main improvement was the noise reduction, especially in the emission field. With $10^{9}$ photons the Monte Carlo results were adequately smooth for the scope of this work and for that reason no further increase to the launch photons number was tested.

A strong similarity between the results from all four models is obvious in Fig. 6. This similarity is also theoretically expected, as this simulation experiment can be adequately addressed even from the diffusion approximation model. A more detailed assessment of the four models is presented in Fig. 7, where the logarithm of the normalized fluence rate and the phase shift on the top surface of the region is shown.

Fig. 7

Logarithm of the normalized fluence rate modulation amplitude (top row) and the phase shift (bottom row) of the excitation (a and c) and the emission (b and d) fields, at the region top surface, along the $y$ -axis (a and b) and along the $x$ -axis (c and d).

The results presented in Fig. 7 show good agreement between all four models, and especially between the dualcoupled, radiative transfer equation, and the Monte Carlo models. Considering the independent Monte Carlo to be the most accurate model, the relative accuracy of the other three was quantified by Eq. (22). In Fig. 8 the errors of the three models versus the Monte Carlo method are presented.

Fig. 8

Absolute values of fluence rate modulation amplitude (top row) and phase shift (bottom row) relative errors of the diffusion approximation, the dual-coupled, and the radiative transfer equation models versus the Monte Carlo forward solver, as a function of the distance from the excitation source along the x-axis (a and b) and the y-axis (c and d) for the excitation field (a and c) and the emission field (b and d).

Quantifying the metrics of Fig. 8 concludes that the radiative transfer equation and the dual-coupled models present more than 92% relative accuracy when compared with Monte Carlo, while this quantity is approximately 88% for the diffusion approximation model. These accuracy levels indicate that the discretization schemes of the three models that are based on the finite element method are sufficiently dense for investigating their applicability as forward solvers in three-dimensional fluorescence molecular imaging. Figure 9 further verifies this assumption, where the logarithms of the normalized fluence rate amplitudes on the plane $y = 1 cm$ are presented for all models. According to the simulation geometry of Fig. 4, this plane is just in front of the ellipsoid simulating the liver. The absorption coefficient of that ellipsoid was responsible for the light extinction observed in Fig. 9 as the black area close to the center of the plots.

Fig. 9

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), the radiative transfer equation (c), and the Monte Carlo (d) models. The logarithm of the normalized fluence rate modulation amplitude, on the $y = 1 cm$ surface of the simulation geometry, of the excitation source is presented in the top row, and for the emission field in bottom row.

4.2.

Simulation Results from the First Two Scenarios

Figure 10 presents the results from the first scenario, showing the logarithm of the fluence rate modulation amplitudes for the diffusion approximation, the dual-coupled radiative transfer equation and diffusion approximation, and the radiative transfer equation models, as slices at the $y = 0$ plane for the excitation and emission fields. The figure shows also the corresponding phase shifts given by Eq. (21b).

Fig. 10

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), and the radiative transfer equation (c) models from the first scenario. The logarithm of the fluence rate modulation amplitude of the excitation source is presented in the plots indexed (1.1), and for the emission field in those indexed (1.2), while the corresponding phase shift is presented in those indexed (2.1) and (2.2).

The fluence rate and phase shift distributions of Fig. 10 can be interpreted by comparing them with the tissue-like digital phantom. Figure 11 shows the excitation and emission fluence rates of the dual-coupled radiative transfer equation and diffusion approximation model, along with the optical properties distribution of the MOBY phantom over the $y = 0$ plane. Light extinction deep inside the discretized volume coincides with the location of the liver.

Fig. 11

The logarithm of the fluence rate modulated amplitude for the excitation (a) and the emission (b) fields, calculated from the dual-coupled radiative transfer equation and diffusion approximation model, and the digital mouse organs (c) in the plane $y = 0$ .

Similar results were obtained from the second simulation scenario, in which the virtual fluorophores were located deeper inside the digital mouse phantom. Figure 12 presents the results obtained from the three forward solvers. All three forward solvers provide similar photon distributions within the inspected tissue.

Fig. 12

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), and the radiative transfer equation (c) models for the second scenario. The logarithm of the fluence rate modulation amplitude is presented in the plots indexed (1.1) for the excitation field, and for the emission in those indexed (1.2), while the corresponding phase shifts are presented in those indexed (2.1) and (2.2).

The fluorescence that is emitted deeply inside the inspected region is almost completely absorbed by the liver, as Fig. 13 shows. The normalized emission fluence rate modulation amplitude is plotted along the $z$ -axis. The peak of fluorescence emission is located at the area of the fluorophore inclusion for all forward solvers, while this signal decreases in the liver region.

Fig. 13

Logarithm of the normalized fluence rate modulation amplitude, and the corresponding phase shift of the emitted fluorescence for the first [(a) fluence rate, (b) phase shift] and the second [(c) fluence rate, (d) phase shift] simulation scenarios along the $z$ -axis.

Results from the dual-coupled radiative transfer equation and diffusion approximation model were more similar to those from the radiative transfer equation for both fluence rate and phase shift, than to those from the diffusion approximation. This is very obvious close to the top surface of the inspected area. The absolute values of the relative errors calculated from Eq. (22) between the models for the second simulation scenario confirm this (Fig. 14).

Fig. 14

Absolute values of fluence rate modulation amplitude [(a) excitation, (b) emission] and phase shift [(c) excitation, (d) emission] relative errors of the diffusion approximation and the dual-coupled models versus the radiative transfer equation results, as a function of the distance from the excitation source along the $z$ -axis, for the second simulation scenario.

The relative errors of both fluence rate and phase shift for the diffusion approximation model versus the radiative transfer equation model were relatively large, especially close to the top surface of the region. Results from the diffusion approximation model and the radiative transfer equation model become more similar deeper inside the region, as expected. In contrast, results from the dual-coupled radiative transfer equation and the diffusion approximation model were quite similar to those from the radiative transfer equation model. More specifically, the accuracy of the dual-coupled model relative to the radiative transfer equation was approximately 97%, while the corresponding accuracy of the diffusion approximation was approximately 90%.

4.3.

Simulation Results from the Third and Fourth Scenarios

The following two simulation scenarios represent two pathological conditions, where the animal kidneys and bladder have been labeled with the fluorophore, Fig. 5(c) and 5(d). The kidneys of the animal model were labeled in the third scenario. The most interesting result is the exitance from the top surface of the model. Figure 15 presents this signal for all three forward solvers.

Fig. 15

Logarithm of the fluence rate modulation amplitude (top row) and phase shift (bottom row) of the diffusion approximation [(a) and (d)], the dual-coupled [(b) and (e)], and the radiative transfer equation [(c) and (f)] models, as calculated in the third simulation scenario.

Figure 15 shows that none of the forward solvers presents two peaks of emitted fluorescence fluence rate, even though the fluorophore is present in both kidneys. The reason for this is that one of the kidneys is located deeper inside the mouse, close to the pancreas, which absorbs a large fraction of the excitation light and the emitted fluorescence. This can be seen in Fig. 16, which presents the dual-coupled radiative transfer equation and diffusion approximation model fluence rate on the $y = 0$ plane, together with the corresponding slice of the animal model. The fluorescence distribution on the top surface shows that the labeled inclusion is not a simple shape. A scanning excitation source is therefore required to reveal both inclusions, Fig. 17.

Fig. 16

The logarithm of the fluence rate modulated amplitude for the excitation (a) and the emission (b) fields, calculated by the dual-coupled radiative transfer equation and diffusion approximation model, and to the digital mouse organs (c), on the plane $y = 0$ .

Fig. 17

Logarithm of the fluence rate modulation amplitude (top row) and phase shift (bottom row) on the top surface of the region, calculated by the dual-coupled model in the third simulation scenario, with the excitation source located at positions $(- 0.6, 0, 1)$ [(a) and (d)], $(0, 0, 1)$ [(b) and (e)], and $(0.6, 0, 1)$ [(c) and (f)], directed inwards. The scanning process revealed both inclusions.

The fourth scenario involved labeling the mouse bladder with fluorophore. Figure 18 presents results from the three forward solvers as slices through the three-dimensional coordinate system. The geometry presented in Fig. 3 was modified for this scenario such that the bladder was located within the discretized region and the z-axis was reversed. Figure 19 shows the logarithm of the fluence rate modulation amplitude for excitation and emission, and a slice of the animal model on the $y = 0$ plane, in the modified geometry.

Fig. 18

Results obtained by solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), and the radiative transfer equation (c) models for the fourth scenario. The logarithm of the fluence rate modulation amplitude is presented in the plots indexed (1.1) for the excitation field, and for the emission in those indexed (1.2), while the corresponding phase shifts are presented in those indexed (2.1) and (2.2). For each index the coloscale remains the same.

Fig. 19

The logarithm of the fluence rate modulated amplitude for the excitation (a) and the emission (b) fields, calculated by the dual-coupled radiative transfer equation and diffusion approximation model, and the digital mouse organs (c) on the plane $y = 0$ .

The diffusion approximation and the dual-coupled models for these two scenarios were compared using a similar method as that used for the previous two scenarios. More specifically, the absolute values of the errors relative to the results from the radiative transfer equation model were estimated using Eq. (22). These errors, however, were estimated in the third and fourth scenarios over the entire discretized region and over the top surface of the region, in order to evaluate the differences between exitances calculated by the three models. The mean relative accuracy was approximately 90% for the entire region, and 98% for the top surface, after applying the dual-coupled model in both scenarios. The corresponding values for the diffusion approximation were 89% and 85%.

The relative accuracies of both models were almost equal when considering the entire region. They were also close to the relative accuracies of the diffusion approximation in the previous two scenarios. The discretization procedure adopted for the dual-coupled model (Fig. 3) resulted in approximately 11,000 spatial nodes for the radiative transfer equation and almost 59,000 nodes for the diffusion approximation, and thus we expected that the relative accuracy when considering the entire region would approach the accuracy obtained from the diffusion approximation model. The accuracy of the exitance, however, calculated by the dual-coupled model was very similar to the accuracy obtained by the radiative transfer equation forward solver, in contrast to the diffusion approximation, for which the accuracy level was approximately 15% smaller than those obtained from the dual-coupled and the radiative transfer equation models.

4.4.

Simulation Results from the Fifth and Sixth Scenarios

Only the dual-coupled radiative transfer equation and diffusion approximation and the radiative transfer equation models addressed the last two simulation scenarios. The diffusion approximation was not used for these two scenarios due to the low scattering index matched interface. Figure 20 shows the logarithm of the normalized fluence rate modulation amplitudes resulted from the fifth simulation scenario, as slices at the $y = 0$ plane for the excitation and emission fields. In this figure are also presented the corresponding phase shifts. The equivalent results from the sixth scenario are presented in Fig. 21.

Fig. 20

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), and the radiative transfer equation (c) models for the fifth scenario. The logarithm of the fluence rate modulation amplitude is presented in the plots indexed (1.1) for the excitation field, and for the emission in those indexed (1.2), while the corresponding phase shifts are presented in those indexed (2.1) and (2.2).

Fig. 21

Results from solving the diffusion approximation (a), the dual-coupled radiative transfer equation and diffusion approximation (b), and the radiative transfer equation (c) models for the sixth scenario. The logarithm of the fluence rate modulation amplitude is presented in the plots indexed (1.1) for the excitation field, and for the emission in those indexed (1.2), while the corresponding phase shifts are presented in those indexed (2.1) and (2.2).

The two forward solvers present almost identical results for both these simulation scenarios. Moreover, as it can be observed in Figs. 20 and 21, both of them preserved the propagation direction of the excitation source inside the low scattering interface, while the light became diffused once entered the mouse tissue at $z = 0.8 cm$ . This is a great advantage of the hybrid model, since it can be used in situations where the diffusion approximation is not a valid model and yet present results very close to the radiative transfer equation and the Monte Carlo method under a much more computational efficient modus.

Nevertheless, the forward solvers studied in this work were developed without considering time efficiency, and all of them are based on similar algorithmic routines. The dual-coupled model required approximately 2 h to complete a simulation, which is significantly faster than the radiative transfer equation model, which required approximately 8 h. The diffusion approximation model gave the solution in the shortest time, approximately 100 s. However, time comparisons at this stage of the research are misleading, since most of the time consumed by the dual-coupled and the radiative transfer equation models is required to solve the phase function integral and to assemble the final finite element matrix. The diffusion approximation model will not become more efficient if more physical memory or a faster HDD is used, but such measures will significantly reduce the time required by the other two models. The radiative transfer equation model, for example, has been additionally solved on a computational unit with 8 GB RAM, which is 12 GB less than that used for the results presented above. This 12 GB reduction in physical memory increased nine times the time required to solve the forward problem, since larger amounts of data needed to be written to virtual memory. For these reasons, no time comparisons were also considered between the three forward solvers and the Monte Carlo method, which required approximately 14 h to solve the forward problem for $10^{9}$ photons.

5. Conclusions

We have applied the coupled radiative transfer equation and diffusion approximation²⁵^,³⁶ model, modified such that it constitutes a dual-coupled fluorescence molecular imaging forward solver, to a digital mouse for the first time. Furthermore, this paper is among the very few that present three-dimensional results from solving the radiative transfer equation.²²^,²⁷^,⁴⁵ An independent Monte Carlo method was used to evaluate the proposed model. Our results were in good agreement with those obtained by the Monte Carlo, indicating accurate modeling of photon migration through scattering media. The robustness of our results is further increased by the use of the MOBY phantom.²⁹ Results from the radiative transfer equation were considered to be the ground truth against which the proposed model and the diffusion approximation were compared.

We have shown that the synthetic data obtained from the frequency domain of the dual-coupled model are well correlated with radiative theory, while the diffusion approximation fails to provide accurate results close to the top boundary of the region. The excitation source was located at this top boundary and emission was detected there in the work presented here. The proposed model gives results that are twice as accurate in photon propagation modeling. More specifically, the relative accuracy of the suggested model is greater than 98%, versus the radiative transfer equation model, close to the top surface of the region. The corresponding accuracy of the diffusion approximation model is almost 85%, as resulted from the first four simulation scenarios presented in this paper. The last two simulation scenarios considered situations where the diffusion approximation fails to accurately simulate light transport. The proposed model, however, achieved equivalent results to the radiative transfer equation, but in a much more computational feasible modus.

Future work will investigate the time efficiency of the proposed model. Thus, one major project will optimize the algorithms with respect to time efficiency and memory management. The time efficiency of the three models can then be compared. Moreover, adaptive unstructured discretization schemes, both spatial and angular, may be used, and the usefulness of the radiative transfer equation model as a forward solver will be investigated. Other work will investigate the possibility of using GPU (graphics processing unit) programming for the three-dimensional solution of the proposed model, and its comparison with the GPU Monte Carlo simulations.

Finally, this forward solver will be used to solve the inverse problem in fluorescence molecular imaging, and results obtained with it will be compared with results obtained using the diffusion approximation and Monte Carlo forward solvers. The long-term scope of this research is to develop a forward solver that combines accuracy and time efficiency, providing an alternative to the widely adopted diffusion approximation and Monte Carlo models.

Acknowledgments

The authors would like to thank Prof. W. P. Segars from Duke University for the license to use the MOBY phantom, which made it possible to implement this research and investigate the appropriateness of the dual-coupled radiative transfer and diffusion approximation model as a forward solver in more realistic scenarios. Stefan Andersson-Engels also acknowledges financial support from the Swedish Research Council and through the Linnaeus grant to Lund Laser Centre, LLC.

References

1.

E. M. Sevick-Muraca, “Translation of near-infrared fluorescence imaging technologies: emerging clinical applications,” Annu. Rev. Med., 63 (1), 217 –231 (2012). http://dx.doi.org/10.1146/annurev-med-070910-083323 ARMCAH 0066-4219 Google Scholar

2.

F. StukerJ. RipollM. Rudin, “Fluorescence molecular tomography: principles and potential for pharmaceutical research,” Pharmaceutics, 3 (2), 229 –274 (2011). http://dx.doi.org/10.3390/pharmaceutics3020229 PHARK5 1999-4923; Google Scholar

3.

S. Dufortet al., “Optical small animal imaging in the drug discovery process,” Biochim. Biophys. Acta., 1798 (12), 2266 –2273 (2010). http://dx.doi.org/10.1016/j.bbamem.2010.03.016 BBACAQ 0006-3002 Google Scholar

4.

J. Halleret al., “Visualization of pulmonary inflammation using noninvasive fluorescence molecular imaging,” J. Appl. Physiol., 104 (3), 795 –802 (2008). http://dx.doi.org/10.1152/japplphysiol.00959.2007 JAPYAA 0021-8987 Google Scholar

5.

S. Choet al., “Optical imaging techniques for the study of malaria,” Trends Biotechnol., 30 (2), 71 –79 (2012). http://dx.doi.org/10.1016/j.tibtech.2011.08.004 TRBIDM 0167-7799 Google Scholar

6.

M. J. Suteret al., “Intravascular optical imaging technology for investigating the coronary artery,” J. Am. Coll. Cardiol. Img., 4 (9), 1022 –1039 (2011). http://dx.doi.org/10.1016/j.jcmg.2011.03.020 1936-878X Google Scholar

7.

S. L. Luoet al., “A review of NIR dyes in cancer targeting and imaging,” Biomaterials, 32 (29), 7127 –7138 (2011). http://dx.doi.org/10.1016/j.biomaterials.2011.06.024 BIMADU 0142-9612 Google Scholar

8.

A. Da Silvaet al., “In vivo fluorescence molecular optical imaging: from small animal towards clinical applications,” in Proc. 16th Mediterranean Conf. on Control and Automation, 1335 –1340 (2008). Google Scholar

9.

M. HaaseH. Schafer, “Upconverting nanoparticles,” Angew. Chem. Int. Ed., 50 (26), 5808 –5829 (2011). http://dx.doi.org/10.1002/anie.v50.26 ACIEF5 1433-7851 Google Scholar

10.

C. T. Xuet al., “High-resolution fluorescence diffuse optical tomography developed with nonlinear upconverting nanoparticles,” ACS Nano, 6 (6), 4788 –4795 (2012). http://dx.doi.org/10.1021/nn3015807 1936-0851 Google Scholar

11.

H. C. LiuC. T. XuS. Andersson-Engels, “Multibeam fluorescence diffuse optical tomography using upconverting nanoparticles,” Opt. Lett., 35 (5), 718 –720 (2010). http://dx.doi.org/10.1364/OL.35.000718 OPLEDP 0146-9592 Google Scholar

12.

S. R. Arridge, “Methods in diffuse optical imaging,” Phil. Trans. R. Soc. A, 369 (1955), 4558 –4576 (2011). http://dx.doi.org/10.1098/rsta.2011.0311 PTRMAD 1364-503X Google Scholar

13.

Y. Luet al., “A parallel adaptive finite element simplified spherical harmonics approximation solver for frequency domain fluorescence molecular imaging,” Phys. Med. Biol., 55 (16), 4625 –4645 (2010). http://dx.doi.org/10.1088/0031-9155/55/16/002 PHMBA7 0031-9155 Google Scholar

14.

A. Ishimaru, Wave Propagation and Scattering in Random Media. Vol.1 Single Scattering and Transport Theory, Academic Press, New York (1978). Google Scholar

15.

S. R. ArridgeJ. C. Hebden, “Optical imaging in medicine: II. Modelling and reconstruction,” Phys. Med. Biol., 42 (5), 841 –853 (1997). http://dx.doi.org/10.1088/0031-9155/42/5/008 PHMBA7 0031-9155 Google Scholar

16.

P. F. Tianet al., “Monte Carlo simulation of the spatial resolution and depth sensitivity of two-dimensional optical imaging of the brain,” J. Biomed. Opt., 16 (1), 016006 (2011). http://dx.doi.org/10.1117/1.3533263 JBOPFO 1083-3668 Google Scholar

17.

J. ChenX. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys., 38 (10), 5788 –5798 (2011). http://dx.doi.org/10.1118/1.3641827 MPHYA6 0094-2405 Google Scholar

18.

G. Quanet al., “Monte Carlo-based fluorescence molecular tomography reconstruction method accelerated by a cluster of graphic processing units,” J. Biomed. Opt., 16 (2), 026018 (2011). http://dx.doi.org/10.1117/1.3544548 JBOPFO 1083-3668 Google Scholar

19.

A. D. KloseV. NtziachristosA. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys., 202 (1), 323 –345 (2005). http://dx.doi.org/10.1016/j.jcp.2004.07.008 JCTPAH 0021-9991 Google Scholar

20.

M. FrancoeurR. VaillonD. R. Rousse, “Theoretical analysis of frequency and time-domain methods for optical characterization of absorbing and scattering media,” J. Quant. Spectrosc. Radiat. Transfer, 93 (1–3), 139 –150 (2005). http://dx.doi.org/10.1016/j.jqsrt.2004.08.018 JQSRAE 0022-4073 Google Scholar

21.

A. Joshiet al., “Radiative transport-based frequency-domain fluorescence tomography,” Phys. Med. Biol., 53 (8), 2069 –2088 (2008). http://dx.doi.org/10.1088/0031-9155/53/8/005 PHMBA7 0031-9155 Google Scholar

22.

P. S. Mohanet al., “Variable order spherical harmonic expansion scheme for the radiative transport equation using finite elements,” J. Comput. Phys., 230 (19), 7364 –7383 (2011). http://dx.doi.org/10.1016/j.jcp.2011.06.004 JCTPAH 0021-9991 Google Scholar

23.

T. Tarvainenet al., “Coupled radiative transfer equation and diffusion approximation model for photon migration in turbid medium with low-scattering and non-scattering regions,” Phys. Med. Biol., 50 (20), 4913 –4930 (2005). http://dx.doi.org/10.1088/0031-9155/50/20/011 PHMBA7 0031-9155 Google Scholar

24.

G. BalY. Maday, “Coupling of transport and diffusion models in linear transport theory,” Math. Model. Numer. Anal., 36 (1), 69 –86 (2002). http://dx.doi.org/10.1051/m2an:2002007 RMMAEV 0764-583X Google Scholar

25.

T. Tarvainenet al., “Finite element model for the coupled radiative transfer equation and diffusion approximation,” Int. J. Numer. Meth. Engng., 65 (3), 383 –405 (2006). http://dx.doi.org/10.1002/(ISSN)1097-0207 IJNMBH 0029-5981 Google Scholar

26.

P. Y. Liu, “A new phase function approximating to Mie scattering for radiative transport-equations,” Phys. Med. Biol., 39 (6), 1025 –1036 (1994). http://dx.doi.org/10.1088/0031-9155/39/6/008 PHMBA7 0031-9155 Google Scholar

27.

D. GorpasD. YovaK. Politopoulos, “A three-dimensional finite elements approach for the coupled radiative transfer equation and diffusion approximation modeling in fluorescence imaging,” J. Quant. Spectrosc. Radiat. Transfer, 111 (4), 553 –568 (2010). http://dx.doi.org/10.1016/j.jqsrt.2009.11.006 JQSRAE 0022-4073 Google Scholar

28.

H. Liet al., “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo method,” Acad. Radiol., 11 (9), 1029 –1038 (2004). http://dx.doi.org/10.1016/j.acra.2004.05.021 1076-6332 Google Scholar

29.

W. P. Segarset al., “Development of a 4-D digital mouse phantom for molecular imaging research,” Mol. Imaging. Biol., 6 (3), 149 –159 (2004). http://dx.doi.org/10.1016/j.mibio.2004.03.002 1536-1632 Google Scholar

30.

W. HongkaiD. B. StoutA. F. Chatziioannou, “Estimation of mouse organ locations through registration of a statistical mouse atlas with micro-CT images,” IEEE Trans. Med. Imaging, 31 (1), 88 –102 (2012). http://dx.doi.org/10.1109/TMI.2011.2165294 ITMID4 0278-0062 Google Scholar

31.

B. Woutjanet al., “Murine cardiac images obtained with focusing pinhole SPECT are barely influenced by extra-cardiac activity,” Phys. Med. Biol., 57 (3), 717 (2012). http://dx.doi.org/10.1088/0031-9155/57/3/717 PHMBA7 0031-9155 Google Scholar

32.

L. Keum Silet al., “MR-based keyhole SPECT for small animal imaging,” Phys. Med. Biol., 56 (3), 685 (2011). http://dx.doi.org/10.1088/0031-9155/56/3/010 PHMBA7 0031-9155 Google Scholar

33.

D. ClarkG. A. JohnsonC. T. Badea, “Denoising of 4D cardiac micro-CT data using median-centric bilateral filtration,” Proc. SPIE, 8314 83143Z (2012). http://dx.doi.org/10.1117/12.911478 Google Scholar

34.

X. ZhangC. T. BadeaG. A. Johnson, “Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy,” J. Biomed. Opt., 14 (6), 064010 (2009). http://dx.doi.org/10.1117/1.3258836 JBOPFO 1083-3668 Google Scholar

35.

D. GorpasS. Andersson-Engels, “Dual coupled radiative transfer equation and diffusion approximation for the solution of the forward problem in fluorescence molecular imaging,” Proc. SPIE, 8225 822522 (2012). http://dx.doi.org/10.1117/12.907811 Google Scholar

36.

T. Tarvainenet al., “Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model,” J. Quant. Spectrosc. Radiat. Transfer, 112 (16), 2600 –2608 (2011). http://dx.doi.org/10.1016/j.jqsrt.2011.07.008 JQSRAE 0022-4073 Google Scholar

37.

R. AronsonN. Corngold, “Photon diffusion coefficient in an absorbing medium,” J. Opt. Soc. Am. A, 16 (5), 1066 –1071 (1999). http://dx.doi.org/10.1364/JOSAA.16.001066 JOAOD6 0740-3232 Google Scholar

38.

H. DehghaniB. Pogue, “Near infrared spectroscopic imaging: theory,” Alternative Breast Imaging, The International Series in Engineering and Computer Science, 183 –199 Springer(2005). Google Scholar

39.

L. G. HenyeyJ. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J., 93 70 –83 (1941). http://dx.doi.org/10.1086/144246 ASJOAB 0004-637X Google Scholar

40.

E. D. AydinC. R. E. de OliveirabA. J. H. Goddard, “A finite element-spherical harmonics radiation transport model for photon migration in turbid media,” J. Quant. Spectrosc. Radiat. Transfer, 84 (3), 247 –260 (2004). http://dx.doi.org/10.1016/S0022-4073(03)00180-8 JQSRAE 0022-4073 Google Scholar

41.

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 (2), R41 –R93 (1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

42.

M. Schweigeret al., “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys., 22 (11 Pt. 1), 1779 –1792 (1995). http://dx.doi.org/10.1118/1.597634 MPHYA6 0094-2405 Google Scholar

43.

P. P. SilvesterR. L. Ferrari, Finite Elements for Electrical Engineers, 3rd ed.Cambridge University Press, New York (1996). Google Scholar

44.

S. Richlinget al., “Radiative transfer with finite elements—I. Basic method and tests,” Astron. Astrophys., 380 (2), 776 –788 (2001). http://dx.doi.org/10.1051/0004-6361:20011411 AAEJAF 0004-6361 Google Scholar

45.

A. D. Klose, “The forward and inverse problem in tissue optics based on the radiative transfer equation: a brief review,” J. Quant. Spectrosc. Radiat. Transfer, 111 (11), 1852 –1853 (2010). http://dx.doi.org/10.1016/j.jqsrt.2010.01.020 JQSRAE 0022-4073 Google Scholar

46.

W. F. CheongS. A. PrahlA. J. Welch, “A review of the optical-properties of biological tissues,” IEEE J. Quantum. Electron., 26 (12), 2166 –2185 (1990). http://dx.doi.org/10.1109/3.64354 IEJQA7 0018-9197 Google Scholar

47.

R. SrinivasanD. KumarM. Singh, “Optical tissue-equivalent phantoms for medical imaging,” Trends Biomater. Artif. Organs, 15 (2), 42 –47 (2002). 0971-1198 Google Scholar

48.

H. Grissaet al., “Three-dimensional radiative transfer modeling using the control volume finite element method,” J. Quant. Spectrosc. Radiat. Transfer, 105 (3), 388 –404 (2007). http://dx.doi.org/10.1016/j.jqsrt.2006.12.008 JQSRAE 0022-4073 Google Scholar

49.

P. Coelho, “A hybrid finite volume/finite element discretization method for the solution of the radiative heat transfer equation,” J. Quant. Spectrosc. Radiat. Transfer, 93 (1–3), 89 –101 (2005). http://dx.doi.org/10.1016/j.jqsrt.2004.08.014 JQSRAE 0022-4073 Google Scholar

50.

Z. Chen, Finite Element Methods and Their Applications, Scientific Computation, 1st ed.Springer, Berlin, New York (2005). Google Scholar

51.

H. A. Van der Vorst, “Bi-CGStab—a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear-systems,” J. Sci. Stat. Comput., 13 (2), 631 –644 (1992). http://dx.doi.org/10.1137/0913035 SJOCE3 1064-8275 Google Scholar

52.

D. FoliniR. Walder, “3D radiative transfer under conditions of non local thermodynamic equilibrium: a contribution to the numerical solution,” Hyperbolic Problems: Theory, Numerics, Applications, International Series of Numerical Mathematics, 305 –314 Birkhauser Verlag, Basel, Switzerland (1999). Google Scholar

53.

H. ShenG. Wang, “A tetrahedron-based inhomogeneous Monte Carlo optical simulator,” Phys. Med. Biol., 55 (4), 947 –962 (2010). http://dx.doi.org/10.1088/0031-9155/55/4/003 PHMBA7 0031-9155 Google Scholar

Citation Download Citation

Dimitris S. Gorpas and Stefan Andersson-Engels "Evaluation of a radiative transfer equation and diffusion approximation hybrid forward solver for fluorescence molecular imaging," Journal of Biomedical Optics 17(12), 126010 (3 December 2012). https://doi.org/10.1117/1.JBO.17.12.126010

Published: 3 December 2012

Access the abstract

JOURNAL ARTICLE
19 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY