Evaluation of path-history-based fluorescence Monte Carlo method for photon migration in heterogeneous media

Xu Jiang; Yong Deng; Zhaoyang Luo; Kan Wang; Lichao Lian; Xiaoquan Yang; Igor Meglinski; Qingming Luo

doi:10.1364/OE.22.031948

1. Introduction

Fluorescence imaging technology is important in biomedicine and is widely used in cancer diagnosis [1,2], drug discovery [3], and gene expression visualization. With the development of the fluorescence probe technique, fluorescence imaging technology has begun to be applied to tracking and measuring the routine activities of specific biological macromolecules inside the small animal model in vivo [4]. Under the action of an external excitation light source, fluorophores inside the biological tissue absorb specific wavelengths of light and then emit fluorescence; the detected fluorescence spectra can provide functional information on the molecule [5]. Therefore, in the past 20 years, fluorescence imaging has been combined with molecular probe techniques for vivo imaging of biological tissues, and it has attracted widespread interest [6–8].

Most biological tissues are highly scattering turbid media. It is essential to establish a highly accurate and efficient computational method for quantitative fluorescence tomography imaging [9].The Monte Carlo method is a discrete theoretical statistical method based on a random sampling procedure [10]. Unlike other methods, the Monte Carlo method can simulate photon transport processes with an arbitrary geometry, arbitrary boundary conditions, and arbitrary optical parameters [11–14]. Because of its broad applicability, the Monte Carlo method is considered the most direct, effective, and credible way to simulate photon transport processes [15]. For this reason, it is viewed as the gold standard method of evaluating other specific application methods such as the diffusion approximation method [16,17]. However, highly accurate Monte Carlo simulations are time-consuming [18]. As the light transmission distance increases, the intensity of the light transmission in the medium decays exponentially. To reach the same level of calculation accuracy, the time required for a Monte Carlo simulation increases exponentially with the size of the medium. The large number of calculations often places a greater burden on the computer.

Many groups have extended the Monte Carlo model to simulate fluorescence in biological tissue because of the growing interest in imaging for medical applications or cancer diagnosis [19–24]. Welch et al. described the laws of fluorescence excitation and propagation in biological tissue and proposed the standard fluorescence Monte Carlo method (sfMC) [19,20]. The simulation results for a semi-infinite layered turbid medium obtained by this method have proven to be accurate [25,26]. However, it often requires a complete sfMC simulation to obtain the results when the optical parameters in fluorescence tomography imaging reconstruction are changed. The calculation is always very large. Moreover, a large number of source–detector pairs must often be simulated in the forward Monte Carlo simulation of the reconstruction, which greatly increases the computation time. Therefore, the sfMC method is not suitable for widespread use [27]. A number of groups have focused on accelerating the fluorescence Monte Carlo method. In recent years, the efficient path-history-based fluorescence Monte Carlo method has attracted increasing attention [28–33]. This method can directly calculate the results of a Monte Carlo simulation by storing photon path histories when the optical parameters are changed, which greatly reduces the number of calculations.

Liebert et al. first proposed an efficient Monte Carlo algorithm for the simulation of time-resolved fluorescence in a layered turbid medium [29–31]. The analytical relationship between the transmitted photon weight and the tissue optical parameters is established by storing the photon path histories. Chen et al. derived a computationally efficient Monte-Carlo-based method of computing time-gated fluorescence Jacobians for the simultaneous imaging of two fluorophores with lifetime contrast [32]. The background weight matrix for the excitation source and detector at a given time can be easily calculated using the photon paths and absorption coefficients. Kumar et al. presented a multiexponential model for time-resolved fluorescence that allows the use of an absorption-perturbation Monte Carlo approach based on stored photon path histories [28,33]. The above three methods are collectively referred to as the perturbation fluorescence Monte Carlo method (pfMC). These pfMC methods involve some assumptions, such as that the fluorescence photons are launched anisotropically instead of isotropically and follow the random walk trajectories of the excitation photons after generation. The pfMC method could improve the computational efficiency of the forward Monte Carlo method when the optical parameters of the tissue change little.

In this work, we compare the theoretical results of a sfMC simulation with experimental results and validate the sfMC method in turbid media. We describe the decoupled fluorescence Monte Carlo (dfMC) method for direct computation of the fluorescence. The dfMC method uses the path information of the excitation light to calculate the fluorescence excitation and transmission process. We demonstrate that the pfMC methods proposed by Liebert et al., Chen et al., and Kumar et al. are equivalent to each other in Continuous Wave (CW) if the absorption coefficient of the background medium in the fluorescent zones is negligible. Furthermore, we compare the pfMC method with the dfMC method using the sfMC method as a reference and report the effect of the different tissue optical parameters and forward model parameters on the calculation accuracy and computational efficiency of the two methods, as well as the applicability of the two methods. We use an accelerated implementation of the sfMC, pfMC, and dfMC methods on GPU clusters.

2. Methods

2.1 Three types of fluorescence Monte Carlo method

2.1.1 sfMC method

The Monte Carlo method can only simulate light excitation and propagation in the past. Welch et al. introduced it into the field of fluorescence and proposed the sfMC model [19,20]. Monte Carlo methods of simulating the fluorescence photon transmission can be divided into four steps. In the first step, for each excitation photon, there is a fixed weight, which is initialized to 1. All the photons are launched from the surface of a medium with diffuse reflectance $R_{x}$ . In the second step, the transmission process of excitation photons is tracked, where the sampling step size is set to $s = - I n (ζ) / μ_{t}$ , and $μ_{t}$ is the extinction coefficient. $ζ$ is a uniform distribution random number in the range (0, 1). After a migration step length, $μ_{a}^{e x} / (μ_{a}^{e x} + μ_{s}^{e x})$ of the photons are absorbed, where $μ_{a}^{e x}$ and $μ_{s}^{e x}$ are the absorption coefficient and scattering coefficient at the excitation wavelength $λ_{x}$ , respectively. In the third step, the process of fluorescence excitation is considered to be isotropic; $μ_{a f} / (μ_{a}^{e x} + μ_{s}^{e x} + μ_{a f})$ of the photons are absorbed to excite fluorescence, and $μ_{a}^{e x} / (μ_{a}^{e x} + μ_{s}^{e x} + μ_{a f})$ of the photons are absorbed but do not excite fluorescence, where $μ_{a f}$ is the absorption coefficient of the fluorophore at the emission wavelength $λ_{m}$ . In the fourth step, the transmission of fluorescence photons is tracked; after each fluorescence photon makes a step, $μ_{a}^{e m} / (μ_{a}^{e m} + μ_{s}^{e m})$ of the photons are absorbed, where $μ_{a}^{e m}$ and $μ_{s}^{e m}$ are the absorption coefficient and scattering coefficient at the emission wavelength $λ_{m}$ , respectively. The remaining photons are tracked until termination or exit from the boundary.

2.1.2 dfMC method

In this section, we describe the dfMC model. By decoupling the excitation-to-emission and transport process of fluorescence from the path probability density function and associating the corresponding parameters involving the fluorescence process with the weight function, the dfMC model uses the path information of the excitation light to calculate the process of fluorescence excitation and transmission. The fluorescence intensity $D$ detected at the detector can be expressed as

\begin{array}{l} D = \int y_{m} (p) g (p) d p \\ = \sum_{m = 0}^{\infty} \int \dots \int τ_{m} (p) w_{m} (p) \prod_{k = 0}^{m} d p_{k} d p, \end{array}

where

p

is a point in six-dimensional phase space, which includes information on the spatial coordinates and direction. It represents any state of photon transmission in tissue. Further,

m

is the number of collisions experienced before the photon reaches a detector, and

g (p)

is the statistical function of the fluorescence detected at the detectors. The statistical function is known. After a photon has experienced

m

collisions,

y_{m} (p)

is the emission density at

p

,

τ_{m} (p)

is the probability density function of photons at

p

, and

w_{m} (p)

is the weight function of photons at

p

.

By using the radiative transfer equation, the emission density of excitation light $x (p)$ can be expressed as

x (p) = \int^{​} x (p') K (p' \to p, μ_{s}^{e x}, μ_{a}^{e x} + μ_{a f}, g^{e x}) d p' + S (p) .

The emission density of fluorescence $y (p)$ can be expressed as

y (p) = \int y (p') K (p' \to p, μ_{s}^{e m}, μ_{a}^{e m}, g^{e m}) d p' + \int x (p') K^{x m} (p' \to p, μ_{a f}) d p',

where

g^{e x}

is the anisotropy constant at the excitation wavelength

λ_{x}

,

g^{e m}

is the anisotropy constant at the emission wavelength

λ_{m}

,

S (p)

is the source emission density;

p'

and

p

are the states of photon transmission, and

p'

is before

p

. The kernel function

K

describes the transitions from state

p'

to

p

; it is the product of the collision kernel and the transport kernel [15,34];

K^{x m}

is the kernel function of the fluorescence excitation process.

K (p' \to p, μ_{s}, μ_{a}, g) = T (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}, μ_{a}) C (\hat{s}' \to \hat{s} | \vec{r}, μ_{s}, μ_{a}, g),

where

T

is the transport kernel, which describes the probability density of photons that propagate from

\vec{r}'

and then collide at

\vec{r}

. The transport kernel can be expressed as

T (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}, μ_{a}) = μ_{t} (\vec{r}) \exp (- \int^{​} μ_{t} (\vec{r}' + l \hat{s}) d l) δ (\hat{s} - (\vec{r} - \vec{r}') / | \vec{r} - \vec{r}' |) / | \vec{r} - \vec{r}' |^{2},

where

C

is the collision kernel, which describes the probability of photons moving in the direction

s^{'}

colliding at

\vec{r}

and then propagating along the

s

direction within solid angle element

d Ω

. The collision kernel can be expressed as

C (\hat{s}' \to \hat{s} | \vec{r}, μ_{s}, μ_{a}, g) = η μ_{s} (\vec{r}) P_{A} (\hat{s}' \cdot \hat{s}, g) / μ_{t} (\vec{r}),

where

P_{A}

is the anisotropy scattering phase function. It is the Henyey-Greenstein phase function.

η

is the quantum efficiency. In addition, we have

S (p) = S (p_{0})

.

According to the Neumann series solution of the integral equation,

\begin{array}{l} y_{m} (p) = \sum_{m = 0}^{\infty} \sum_{i = 0}^{m} \int \dots \int S (p_{0}) K (p_{0} \to p_{1}, μ_{s}^{e x}, μ_{a}^{e x} + μ_{a f}, g^{e x}) \dots K (p_{i - 1} \to p_{i}, μ_{s}^{e x}, μ_{a}^{e x} + μ_{a f}, g^{e x}) \\ \times K^{x m} (p_{i} \to p_{i + 1}, μ_{a f}) \times K (p_{i + 1} \to p_{i + 2}, μ_{s}^{e m}, μ_{a}^{e m}, g^{e m}) \dots K (p_{m - 1} \to p_{m}, μ_{s}^{e m}, μ_{a}^{e m}, g^{e m}) \prod_{k = 0}^{m} d p_{k} . \end{array}

In the dfMC method, the kernel function $K^{x m}$ of the fluorescence excitation process can be expressed as

\begin{array}{l} K^{x m} (p' \to p, μ_{a f}) = T (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}^{e x}, μ_{a}^{e x} + μ_{a f}) C^{x m} (\hat{s}' \to \hat{s} | \vec{r}, μ_{a f}) \\ = K (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}) \exp (- \int_{0}^{| \vec{r} - \vec{r}' |} μ_{a f} (\vec{r}' + l \hat{s}) d l) \\ \times \frac{μ_{t}^{e x} (\vec{r}) + μ_{a f} (\vec{r})}{μ_{t}^{e x} (\vec{r})} \frac{C^{x m} (\hat{s}' \to \hat{s} | \vec{r}, μ_{a f})}{C (\hat{s}' \to \hat{s} | \vec{r}, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x})} \\ = K (p' \to p, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}) \exp (- \int_{0}^{| \vec{r} - \vec{r}' |} μ_{a f} (\vec{r}' + l \hat{s}) d l) \frac{η μ_{a f} (\vec{r})}{μ_{s}^{e x} (\vec{r})} \frac{P_{I} (\hat{s}' \cdot \hat{s})}{P_{A} (\hat{s}' \cdot \hat{s}, g^{e x} (\vec{r}'))}, \end{array}

where

P_{I}

is the isotropic scattering phase function, and

μ_{t}^{e x}

is the extinction coefficient at the excitation wavelength

λ_{x}

. If the extinction coefficient at the excitation wavelength differs from that at the emission wavelength, we can perform the following conversion:

\begin{array}{l} K (p' \to p, μ_{s}^{e m}, μ_{a}^{e m}, g^{e m}) = T (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}^{e m}, μ_{a}^{e m}) C (\hat{s}' \to \hat{s} | \vec{r}, μ_{s}^{e m}, μ_{a}^{e m}, g^{e m}) \\ = T (\vec{r}' \to \vec{r} | \hat{s}', μ_{s}^{e x}, μ_{a}^{e x}) \exp (- \int_{0}^{| \vec{r} - \vec{r}' |} (μ_{t}^{e m} (\vec{r}' + l \hat{s}) - μ_{t}^{e x} (\vec{r}' + l \hat{s})) d l) \\ \times C (\hat{s}' \to \hat{s} | \vec{r}, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}) \frac{μ_{s}^{e m} (\vec{r})}{μ_{s}^{e x} (\vec{r})} \\ = K (p' \to p, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}) \exp (- \int_{0}^{| \vec{r} - \vec{r}' |} (μ_{t}^{e m} (\vec{r}' + l \hat{s}) - μ_{t}^{e x} (\vec{r}' + l \hat{s})) d l) \\ \times \frac{μ_{s}^{e m} (\vec{r})}{μ_{s}^{e x} (\vec{r})} \frac{P_{A} (\hat{s}' \cdot \hat{s}, g^{e m} (\vec{r}'))}{P_{A} (\hat{s}' \cdot \hat{s}, g^{e x} (\vec{r}'))}, \end{array}

where

μ_{t}^{e m}

is the extinction coefficient at the emission wavelength

λ_{m}

. The conversion between Eqs. (8) and (9) replaces the kernel function of the fluorescence excitation and fluorescence transmission processes with that of the excitation light transmission process. At the same time, we can obtain

τ_{m} (p) = S (p_{0}) K (p_{0} \to p_{1}, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}) \dots K (p_{m - 1} \to p_{m}, μ_{s}^{e x}, μ_{a}^{e x}, g^{e x}),

\begin{array}{l} w_{m} (p) = g (p) \sum_{i = 0}^{m} (\prod_{j = 0}^{i} \exp (- \int_{0}^{| {\vec{r}}_{j + 1} - {\vec{r}}_{j} |} (μ_{a f} ({\vec{r}}_{j} + l {\hat{s}}_{j + 1}) d l)) \frac{η μ_{a f} ({\vec{r}}_{i})}{μ_{s}^{e x} ({\vec{r}}_{i})} \\ \times \frac{P_{I} ({\hat{s}}_{i} \cdot {\hat{s}}_{i + 1})}{P_{A} ({\hat{s}}_{i} \cdot {\hat{s}}_{i + 1}, g^{e x} ({\vec{r}}_{i}))} \prod_{j = i + 1}^{m} \frac{μ_{s}^{e m} ({\vec{r}}_{j})}{μ_{s}^{e x} ({\vec{r}}_{j})} \frac{P_{A} ({\hat{s}}_{j} \cdot {\hat{s}}_{j + 1}, g^{e m} ({\vec{r}}_{j}))}{P_{A} ({\hat{s}}_{j} \cdot {\hat{s}}_{j + 1}, g^{e x} ({\vec{r}}_{j}))} \\ \times \exp (- \int_{0}^{| {\vec{r}}_{j + 1} - {\vec{r}}_{j} |} (μ_{t}^{e m} ({\vec{r}}_{j} + l {\hat{s}}_{j + 1}) - μ_{t}^{e x} ({\vec{r}}_{j} + l {\hat{s}}_{j + 1})) d l) . \end{array}

2.1.3 pfMC method

Liebert et al. first proposed applying the pfMC method along the path of excitation light using the Beer–Lambert theorem to calculate the fluorescence [29–31]. This method is based mainly on two assumptions: 1) All the excitation and fluorescence trajectories from a given location within the medium to the detector are virtually the same. 2) The scattering is isotropic and can be characterized by a single reduced scattering coefficient [31].

Excitation light is perpendicularly incident on the surface of tissue with the initial weight $w_{i 0}$ . According to the Beer–Lambert theorem, the weight of the $i^{t h}$ photon having a source at $r_{s}$ and fluorescence excitation position at $r$ can be expressed as

w_{i} (r_{s}, r) = w_{i 0} \exp (- \sum_{j = 1}^{p_{i}} (μ_{a}^{e x} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})),

where

r_{j}

(j = 1, \dots, p_{i})

are the subregions that the excitation photon passes through from

r_{s}

to

r

, and

l_{i} (r_{j})

is the path length that the photon travels at

r_{j}

. The excitation-to-emission conversion probability

p_{i} (r)

can be expressed as

p_{i} (r) = Φ (1 - \exp (- μ_{a f} l_{i} (r))),

where

Φ

is the quantum efficiency. The weight of the emitted fluorescence photon at

r

is set to

w_{i 0}^{'}

.The weight of the fluorescence photon detected at

r_{d}

has decreased to

w_{i} (r, r_{d}) = w_{i 0}^{'} \exp (- \sum_{j = p_{i} + 1}^{q_{i}} (μ_{a}^{e m} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})),

where

r_{j}

(j = p_{i} + 1, \dots, q_{i})

are the subregions that the fluorescence photon passes through from

r

to

r_{d}

. Calculating all the photons, we obtain the sum of the weight of the fluorescence photons that originate at

r_{s}

and are detected at

r_{d}

.

\begin{array}{l} W (r_{s}, r_{d}, r) = \sum_{i = 1}^{n} w_{i, 0} \exp (- \sum_{j = 1}^{p_{i}} (μ_{a}^{e x} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) (1 - \exp (- μ_{a f} (r) l_{i} (r))) \\ \times \exp (- \sum_{j = p_{i} + 1}^{q_{i}} (μ_{a}^{e m} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) Φ . \end{array}

Chen et al. derived the pfMC model on the basis of the Born approximation [32]. The fluorescence intensity $U_{F} (r_{s}, r_{d})$ detected by the detector at $r_{d}$ for a source position at $r_{s}$ and excitation position at $r$ can be expressed as

U_{F} (r_{s}, r_{d}) = \int_{Ω} d r G^{x} (r_{s}, r) G^{m} (r, r_{d}) η (r),

where

G^{x}

and

G^{m}

are the Green’s functions for light propagation at the excitation and emission wavelengths, respectively.

η (r)

is the effective quantum yield,

η (r) = \frac{μ_{a f} (r) Φ}{μ_{a}^{x} (r)},

where

μ_{a}^{x} (r)

is the total absorption coefficient of the subregions at the excitation wavelength

λ_{x}

. We define a weighting function

W' (r_{s}, r_{d}, r)

:

U_{F} (r_{s}, r_{d}) = \int_{Ω} d r W' (r_{s}, r_{d}, r) η (r) .

We then calculate the weight of the fluorescence along the path of the excitation light,

\begin{array}{l} W (r_{s}, r_{d}, r) = W' (r_{s}, r_{d}, r) η (r) \\ = \sum_{i = 1}^{n} w_{i, 0} \exp (- \sum_{j = 1}^{p_{i}} μ_{a}^{x} (r_{j}) l_{i} (r_{j})) (1 - \exp (- μ_{a}^{x} (r) l_{i} (r)) \exp (- \sum_{j = p_{i} + 1}^{q_{i}} μ_{a}^{m} (r_{j}) l_{i} (r_{j})) η (r), \end{array}

where

μ_{a}^{m} (r)

is the total absorption coefficient of the subregions at the emission wavelength

λ_{m}

. By using a Taylor expansion, Eqs. (15) and (19) can be written as

\begin{array}{l} W (r_{s}, r_{d}, r) = \sum_{i = 1}^{n} w_{i, 0} \exp (- \sum_{j = 1}^{p_{i}} (μ_{a}^{e x} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) (1 - \exp (- μ_{a f} (r) l_{i} (r))) \\ \times \exp (- \sum_{j = p_{i} + 1}^{q_{i}} (μ_{a}^{e m} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) Φ \\ = \sum_{i = 1}^{n} w_{i, 0} Φ \exp (- \sum_{j = 1}^{p_{i}} (μ_{a}^{e x} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) \exp (- \sum_{j = p_{i} + 1}^{q_{i}} (μ_{a}^{e m} (r_{j}) + μ_{a f} (r_{j})) l_{i} (r_{j})) \\ \times (μ_{a f} (r) l_{i} (r) - \frac{{(μ_{a f} (r) l_{i} (r))}^{2}}{2!} \dots - (^{- 1) n} \frac{{(μ_{a f} (r) l_{i} (r))}^{n}}{n!}), \end{array}

\begin{array}{l} W (r_{s}, r_{d}, r) = W' (r_{s}, r_{d}, r) η (r) \\ = \sum_{i = 1}^{n} w_{i, 0} \exp (- \sum_{j = 1}^{p_{i}} μ_{a}^{x} (r_{j}) l_{i} (r_{j})) \exp (- \sum_{j = p_{i} + 1}^{q_{i}} μ_{a}^{m} (r_{j}) l_{i} (r_{j})) (1 - \exp (- μ_{a}^{x} (r) l_{i} (r)) \frac{μ_{a f} (r) Φ}{μ_{a}^{x} (r)} \\ = \sum_{i = 1}^{n} w_{i, 0} Φ \exp (- \sum_{j = 1}^{p_{i}} μ_{a}^{x} (r_{j}) l_{i} (r_{j})) \exp (- \sum_{j = p_{i} + 1}^{q_{i}} μ_{a}^{m} (r_{j}) l_{i} (r_{j})) \\ \times (μ_{a f} (r) l_{i} (r) - \frac{μ_{a f} (r) μ_{a}^{x} (r) l_{i} {(r)}^{2}}{2!} \dots - (^{- 1) n} \frac{μ_{a f} (r) μ_{a}^{x} {(r)}^{n - 1} l_{i} {(r)}^{n}}{n!}) . \end{array}

From Eqs. (20) and (21), we can see that if $μ_{a f}^{} (r) = μ_{a}^{x} (r)$ , that is, the absorption coefficient of the background medium in the fluorescent zones is negligible, Eq. (15) is equivalent to Eq. (19). Thus, the pfMC model proposed by Chen et al. [32] is equivalent to that proposed by Liebert et al. [31]. The absorption-perturbation Monte Carlo model [33] proposed by Kumar is equivalent to the pfMC model proposed by Liebert et al. in Continuous Wave (CW). Therefore, the three above-mentioned pfMC models can be unified. Note that the fluorescence photons are launched anisotropically instead of isotropically in these three pfMC models.

2.2 Code and configuration of the GPU cluster scheme

2.2.1 Photon transport process

A large number of photons propagate in the medium, excite fluorescence, and produce a sequence of scattering histories and fluence distributions until they exit the medium. This process can be briefly described as:

The incident light source is characterized for the collection of a specific number of photons. A direction and incident light source position are assigned to each photon. The initial weight of each photon is set to 1 [12].
The scattering length, i.e., the distance to the next scattering event site, is computed using the extinction coefficient (in the sfMC method) or the scattering coefficient (in the pfMC and dfMC methods) of the current position on the basis of an exponential distribution.
The photons are moved one substep along the scattering trajectory.
In the sfMC method, after each substep, a new random number $ζ$ which is uniformly distributed in (0, 1) is generated. Whether the photon is absorbed by tissue or a fluorescence photon is excited depends on $ζ$ . In the transmission process of the excitation photon, if $ζ < μ_{a}^{e x} / (μ_{a}^{e x} + μ_{s}^{e x})$ , the excitation photons is absorbed to excited fluorescence. In the process of fluorescence excitation, if $ζ < μ_{a f} / (μ_{a}^{e x} + μ_{s}^{e x} + μ_{a f})$ , the excitation photon is absorbed to excite fluorescence; if $μ_{a f} / (μ_{a}^{e x} + μ_{s}^{e x} + μ_{a f}) < ζ < (μ_{a}^{e x} + μ_{a f}) / (μ_{a}^{e x} + μ_{s}^{e x} + μ_{a f})$ , the excitation photon is absorbed but does not excite fluorescence. In the transmission process of the fluorescence photon, if $ζ < μ_{a}^{e m} / (μ_{a}^{e m} + μ_{s}^{e m})$ , the fluorescence photon is absorbed. If the photon is absorbed but does not excite a fluorescence photon, it is terminated with the weight set to zero. In the pfMC and dfMC methods, the weight of the photon packet is reduced by the sum of the absorption coefficients along the excitation light trajectory. After each substep, the weight of the photon packet is calculated. At the fluorescence excitation position, $μ_{a f} (r) l_{i} (r)$ is calculated in the pfMC method, and $\frac{μ_{a f} (\vec{r})}{μ_{s}^{e x} (\vec{r})} \frac{P_{I}}{P_{A}}$ is calculated in the dfMC method.
Repeat from step 3 until the photon has traveled the total scattering length [18].
At the scattering point, a new scattering direction vector is calculated on the basis of the Henyey–Greenstein phase function.
Repeat from step 2 until all of the photons terminate or exit the boundary.

2.2.2 Computational settings

An accelerated process of tracking the photons is implemented on GPU clusters with the Compute Unified Device Architecture (CUDA) [35–39]. The GPU clusters support Message Passing Interface (MPI) and Open Multi-Processing for parallel computing.

The entire frame is implemented in the C programming language. Computing tasks are assigned to each child node for parallel computing by the MPI communication protocol [37]. To maximize the parallel computing capacity, reduce the computation time of the fluorescence Monte Carlo methods, and take into account the fact that the stored path information limits the display memory size, there are 64 blocks in each GPU, with 622 threads in each block. A sequence of photon migration processes is implemented on each thread.

The GPU cluster we set up includes two child nodes. One child node has a 2.93 GHZ Intel Xeon X5570 CPU with four cores and 12 GB of memory. The child node contains two GPUs. One GPU is a Tesla C1060 with 2 GB of display memory and 240 CUDA cores; the other GPU has a Quadro FX4800 with 1.5 GB of display memory and 192 CUDA cores. The other child node is 3.4 GHZ Intel Core i7 2600 CPU with four cores and 32 GB of memory. The GPU is a GTX580 with 1.5 GB of display memory and 512 CUDA cores.

3. Results

3.1 Experimental verification of the sfMC

3.1.1 Experimental system

The sfMC method is accurate and describes the true process of fluorescence excitation and transmission [25,26]. Welch et al. verified the sfMC method by one- and two-layer simulations [20]. To validate the sfMC method with sufficient accuracy in turbid media, we use the FMT&mCT combined experimental system developed in our lab [40–42]. Figure 1 shows a sketch of the experimental system. The FMT subsystem is mounted orthogonally to the mCT subsystem. The mCT subsystem consists mainly of an x-ray source (UltraBright, Oxford Instruments, USA) and a flat panel x-ray detector (PaxScan2520V, Varian Medical Systems, USA). In a typical mCT scan, the x-ray source is operated with 60 kVp and 40 W. Four hundred projections are acquired to reconstruct the image by the Feldkmap’s algorithm. The FMT subsystem consists mainly of a 748-nm continuous-wave (CW) diode laser (B&W TEK, Newark, Delaware), a dual-axis galvanometer scanner (Galvo Scanner ST8061, Shiji Tuotian, China), an f-theta lens (F160, custom-made, China) and a CCD camera (DU-897, Andor, UK). The dual-axis galvanometer scanner is used to implement the laser beam scanning on the surface of the phantom within a certain range. The f-theta lens is used to focus the laser beam on the surface of the phantom [41]. The phantom is placed on the fixed object stage, which can be rotated $360^{°}$ by the rotating stage (ADRS-100, Aerotech, USA). The CCD camera is fixed on the other side of the galvanometer. A band pass filter (750FS10-25, Andover Corp., Salem, USA) and a long pass filter (HQ770LP, Chroma, USA) are placed in front of the CCD camera to capture the excitation light signal and fluorescence signal, respectively. The image taken by the CCD camera has 512*512 pixels. The computer is used to store and process the data.

Fig. 1 Sketch of the experimental system.

Download Full Size | PDF

The laser beam from the 748-nm continuous-wave diode laser is collimated, expanded, and then conducted into the dual-axis galvanometric scanner. In the dual-axis galvanometer scanner, the mirrors can rotate different angles with changeable input voltages. The laser beam is incident on the center of X-axis mirror, and then it is reflected by the X-axis mirror and directed into Y-axis mirror. Finally, it is reflected by Y-axis mirror and focused on the surface of the phantom by f-theta lens. Since the distance between the X-axis and the Y-axis mirrors is much smaller than the working distance of the f-theta lens, the center of the Y-axis mirror is approximately considered the starting point of the laser beam. The direction vector of the laser beam depends on the input voltage. Combined the starting point and the direction vector of the laser beam with the surface profile of the phantom provided by the mCT subsystem, the spatial coordinates of the laser spot on the surface of the phantom can be determined by a ray tracing method in the code [42].

The spatial coordinates of the light intensity distribution on the surface of the phantom can be obtained by the mCT subsystem, and the projection of the light intensity distribution on the CCD camera can be obtained by the FMT subsystem. To determine the mapping relationship between the intensity distribution on the surface of the phantom and the light intensity distribution on the CCD camera, four coordinate systems are established [42]. Through the coordinate transformation and the pinhole imaging principle, the mapping relationship between the light intensity distribution on the surface of the phantom and the light intensity distribution on the CCD camera can be built to simulate the imaging process in the code.

3.1.2 Cylindrical phantom

The cylindrical phantom consists of a glass box filled with 1% intralipid $(μ_{s} = 100.0 c m^{- 1}, μ_{a} = 0.02 c m^{- 1})$ [43]. One transparent glass tube filled with DiR-BOA fluorescent dye [44] is placed in the glass box. The DiR-BOA fluorescent dye used in the experiment is diluted to a certain concentration. The reflectivity and transmissivity of a cuvette of the DiR-BOA fluorescent dye can be measured by the spectrophotometer (Lamda 950, PerkinElmer, USA) with a single-integrating sphere. Using adding-doubling (IAD) method, the fluorescence coefficient of the DiR-BOA fluorescent dye is calculated from the measured data. The fluorescence coefficient of DiR-BOA fluorescent dye is $1.000 \pm 0.001 c m^{- 1}$ . The cylindrical phantom used in the experiment is shown in Fig. 2. Its diameter and height are 3 cm and 4.2 cm, respectively. In the sfMC simulation, the cylindrical phantom is numerical. It is dissected into voxels each 0.0515 cm in size. We set the refractive index to $n = 1.37$ and use the anisotropic factor $g = 0.9$ for both the fluorescent target areas and the background. The optical parameters used for the simulation are shown in Table 1. The excitation and emission wavelengths usually differ by a few 10 nm only [31]. Therefore, we set the same absorption coefficient and scattering coefficient at the excitation and emission wavelengths.

Fig. 2 The cylindrical phantom is composed of a glass box and fluorescent tube.

Download Full Size | PDF

Table 1. Optical Properties Used in the Cylindrical Phantom

View Table | View all tables in this article

3.1.3 Experimental results

The number of simulated photons is $10^{9}$ in the sfMC simulation. According to the transformation from the FMT coordinate system to the surface pixel coordinate system of the CCD camera images, the fluorescence photons that exit the surface of the cylindrical model are projected onto the CCD camera imaging surface. Figures 3(a) and 3(b) compare the normalized fluorescence intensity of the sfMC simulation and the experiment. Figures 3(c) and 3(d) compare the intensity in more detail along the center line of the CCD camera and as contours.

Fig. 3 Normalized fluorescence intensity on the CCD camera in (a) the sfMC simulation and (b) the experiment. (c) Comparison of the fluorescence intensity along the center line of the CCD camera. (d) Comparison of the contour lines of the fluorescence intensity on the CCD camera in the sfMC simulation and experiment.

Download Full Size | PDF

The comparison reveals that the result of the sfMC simulation agrees with the experimental result. When the number of simulated photons reaches $10^{9}$ , we can obtain a better result. Therefore, we select the simulation result with $10^{9}$ photons as a reference.

3.2 Comparison of the pfMC and dfMC methods

3.2.1 Statistical criteria

The calculation accuracy is a key index for evaluating the fluorescence Monte Carlo method. The spatial distribution of the fluorescence intensity on the detector determines the method’s calculation accuracy. To quantitatively evaluate the path-history-based fluorescence Monte Carlo method, we define an objective error metric $e (r)$ [9], where the error $e$ is at position $r$ .

e (r) = | \frac{v (r) - v_{r e f} (r)}{v_{r e f} (r)} | \times 100 %,

where

v (r)

and

v_{r e f} (r)

are the simulated and observed fluorescence intensity on the detector at

r

, respectively. The average of e for all the detectors is denoted as

\bar{e}

. In previous research, the sfMC method is generally considered to have sufficient precision. We also verified that the result of the sfMC simulation agrees with the experimental result. Therefore,

v_{r e f} (r)

in Eq. (22) can be replaced by the result of the sfMC simulation.

3.2.2 Cylindrical model

We evaluate the path-history-based fluorescence Monte Carlo method through a simulation study. The designed cylindrical model [45–47] is shown in Fig. 4. We consider a medium with four different types of tissue (muscle, bone, lungs, and heart), where the fluorophore is located at the lungs. The fluorescent zones consist of lungs, and the non-fluorescent zones consist of muscle, bone, and heart. The model has 301,401 voxels, each 0.05 cm in size. The cylindrical model is 3.05 cm × 3.05 cm × 4.05 cm in size. The position and direction of the source are shown in Fig. 4. The detectors are placed opposite to the source within a range of $120^{°}$ , which is distributed evenly in 10 layers with 15 detectors in each layer. Various parameters affect the calculation accuracy and computational efficiency of the fluorescence Monte Carlo method, such as the number of photons, fluorescence coefficient, absorption coefficient, and scattering coefficient. These parameters are independent of each other. In each simulation, we change one parameter while keeping the other parameters constant and identify the effect of the change in that parameter on the calculation accuracy and computational efficiency. We list the fixed values of the tissue parameters [46] in Table 2. The other parameters are set to typical values for biological tissue in the near-infrared spectral region: $g$ is set to 0.9, the refractive index $n$ is set to 1.37, and $Φ$ is set to 1.

Fig. 4 Cylindrical model with a source composed of bone (yellow), lungs (blue), heart (red), and muscle (green). The fluorophore is located at the lungs.

Download Full Size | PDF

Table 2. Fixed Value of the Tissue Optical Parameters

View Table | View all tables in this article

3.2.3 Efficiency comparison $\bar{e}$

Changes in the number of photons have important effects on the statistical distribution of the fluorescence intensity at the detector. As shown in Fig. 5, the mean error $\bar{e}$ of the two methods decreases as the number of simulated photons increases. When the number of simulated photons increases from a relatively small value of $(10^{5} - 10^{6})$ to $10^{7}$ , the mean error of the two methods decreases rapidly. However, when the number of photons continues to increase, the mean error $\bar{e}$ of the two methods begins to stabilize. When the number of simulated photons reaches $10^{7}$ , the mean error $\bar{e}$ of the pfMC method becomes constant around 38.3%. The mean error $\bar{e}$ of the dfMC method decreases as the number of simulated photons increases. The dfMC method requires $10^{6}$ photons, compared to $10^{7}$ for the pfMC method, to reach the same level of statistics. Therefore, the dfMC method has a higher computational efficiency than the pfMC method under the same conditions. To reduce the calculation burden, we simulate $10^{8}$ photons in the following comparison.

Fig. 5 Mean error at the detectors versus number of photons for pfMC and dfMC methods.

Download Full Size | PDF

3.2.4 Accuracy comparison

An increase in the fluorescence coefficient will improve the fluorescence excitation efficiency. The changes in the fluorescence intensity on the detector affect the accuracy of the calculation. We changed the fluorescence coefficient of the fluorophore at the lungs; the results are shown in Fig. 6. As the fluorescence coefficient increases, the mean error $\bar{e}$ of the pfMC method decreases, whereas that of the dfMC method changes slightly. We can conclude that an increase in the fluorescence coefficient results in an increase in the calculation accuracy of the pfMC method, and the dfMC method is not sensitive to the fluorescence coefficient. From Fig. 6, we can see that the dfMC method has a higher calculation accuracy than the pfMC method for a wide range of fluorescence coefficients.

Fig. 6 Mean error at the detectors versus fluorescence coefficient for pfMC and dfMC methods.

Download Full Size | PDF

The path-history-based fluorescence Monte Carlo method uses the Beer–Lambert theorem to calculate the fluorescence along the excitation light trajectory. As the absorption coefficient increases, the fluorescence weight decays more rapidly, and the fluorescence intensity detected by the detectors becomes weaker. We changed the absorption coefficients of the lungs, muscle, and bone separately; the results are shown in Fig. 7. As the absorption coefficient increases, the mean error $\bar{e}$ of the dfMC method increases. The mean error $\bar{e}$ of the pfMC method increases if the absorption coefficient of the fluorescent zones (lungs) increases, and the mean error $\bar{e}$ of the pfMC method decreases if the absorption coefficient of the non-fluorescent zones (muscle and bone) increases. We can conclude that changes in the absorption coefficient of the fluorescent zones and that of the non-fluorescent zones have opposite effects on the calculation accuracy of the pfMC method; however, an increase in the absorption coefficient decreases the calculation accuracy of the dfMC method. From Fig. 7, we can see that the dfMC method has a higher calculation accuracy than the pfMC method for a wide range of absorption coefficients.

Fig. 7 Mean error at the detectors versus absorption coefficient of (a) lungs, (b) muscle, and (c) bone for pfMC and dfMC methods.

Download Full Size | PDF

In the path-history-based fluorescence Monte Carlo method, the scattering coefficient determines the size of the sampling step; thus, changes in the scattering coefficient have an important effect on the calculation accuracy. We changed the scattering coefficients of the lungs, muscle, and heart separately, and the results are shown in Fig. 8. As the scattering coefficient increases, the trend in the mean error is uncertain in the pfMC and dfMC methods. We can conclude that changes in the scattering coefficient of the fluorescent zones (lungs) strongly affect the calculation accuracy; however, changes in the scattering coefficient of the non-fluorescent zones (muscle and heart) have a slight effect on the calculation accuracy in the pfMC method. From Fig. 8, we can see that the dfMC method has a higher calculation accuracy than the pfMC method for a wide range of scattering coefficients.

Fig. 8 Mean error at the detectors versus scattering coefficient of (a) lungs, (b) muscle, and (c) heart for pfMC and dfMC methods.

Download Full Size | PDF

4. Discussion and conclusions

In this paper, we verify the calculation accuracy of the sfMC method in a three-dimensional turbid medium by experimenting with a cylindrical phantom. We describe the dfMC method and its implementation. We also investigate the pfMC method. The pfMC methods with different forms of expression are equivalent to each other in Continuous Wave (CW) when the absorption coefficient of the background medium in the fluorescent zones is negligible. We compare the calculation accuracy and computational efficiency of the pfMC and dfMC methods using that of the sfMC method as a reference.

The dfMC and sfMC methods are both derived from the radiative transfer equation [1]. The dfMC method does not make any assumptions. Thus, if the number of simulated photons is sufficient, the results of a dfMC simulation are close to those of a sfMC simulation. Regardless of changes in the optical parameters, the calculation accuracy of the dfMC method remains essentially stable. The pfMC method makes some assumptions in calculating the fluorescence. These assumptions are satisfied under certain conditions. Thus, the calculation accuracy of the pfMC method is limited. The pfMC method always has a lower calculation accuracy than the dfMC method. In the pfMC method, we use the Beer–Lambert theorem to calculate the fluorescence. As the fluorescence coefficient increases, the fluorescence intensity detected by the detectors will increase. The calculation accuracy of the pfMC method will be better for a high fluorescence coefficient. In addition, it is found that as the scattering coefficient of the fluorescent zones varies, the mean error at the detectors is fluctuant in the pfMC method. This is mainly because (i) the fluorescence photons are assumed to be launched anisotropically instead of isotropically in the pfMC method, (ii) the scattering coefficient of the fluorescent zones directly affects the scattering step, the scattering direction and the corresponding scattering probability of the fluorescence photons, and (iii) the scattering step, the scattering direction and the corresponding scattering probability of the fluorescence photons are dependent on the launching mode. If the difference of the scattering step, the scattering direction and as well as the corresponding scattering probability between the fluorescence photons launched anisotropically and isotropically is small, the mean error at the detectors will be small, otherwise the large error will occur. Therefore, such the change of the scattering coefficient of the fluorescent zones will cause a random fluctuation of the mean error at the detectors.

In a previous study, Liebert et al. validated the pfMC method by comparing it with the diffusion approximation method for a homogeneous semi-infinite turbid medium [31]. Chen et al. applied the pfMC method to fluorescence tomography imaging reconstruction and successfully reconstructed two tubes filled with fluorescent dyes in a square box [32]. Kumar validated the pfMC method by simulating the transmission process of photons in both a slab and a heterogeneous model of the mouse head [33]. Altogether, these researchers compared the pfMC method with the diffusion approximation method, the adjoint fluorescence Monte Carlo method [33], or an experiment using the normalized fluorescence intensity. However, using the sfMC method as a reference, we compared the pfMC method with the dfMC method using the fluorescence intensity. It is helpful to study the absolute quantification of fluorescence tomography imaging.

In the pfMC and dfMC methods, the photon’s trajectory is generated only once. When the fluorescence coefficient is changed, we can directly calculate the fluorescence by applying the relations for the stored photon histories. The computation time required to calculate the fluorescence for different optical properties can be greatly reduced.

In conclusion, the dfMC method requires fewer photons than the pfMC method to reach the same level of statistics. The dfMC method has a higher computational efficiency than the pfMC method. The calculation accuracy of the pfMC method is high if the absorption coefficient of the fluorescent zones and the fluorescence coefficient are high. Changes in the scattering coefficient of the fluorescent zones affect the calculation accuracy of the pfMC method more strongly than changes in that of the non-fluorescent zones. However, the dfMC method can remain quite accurate over a wide range of optical parameters. In most cases, the calculation accuracy of the dfMC method is significantly higher than that of the pfMC method. In the future, the computational efficiency of the path-history-based fluorescence Monte Carlo method can be further improved by combining it with the controlled Monte Carlo method [48] or two-step Monte Carlo method [49]. The path-history-based fluorescence Monte Carlo method is potentially useful for fluorescence-based in vivo tomography.

Acknowledgments

This project was supported by the National Major Scientific Research Program of China (No. 2011CB910401), National Key Technology R&D Program (Grant No. 2012BAI23B02), Science Fund for Creative Research Group (Grant No. 61121004), and National Natural Science Fund (Grant No. 61078072).

References and links

1. J. Lakowicz, Principle of Fluorescence Spectroscopy (Springer, 2006).

2. A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R. Arridge, M. D. Schnall, and A. G. Yodh, “Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express 15(11), 6696–6716 (2007). [CrossRef] [PubMed]

3. V. Ntziachristos, E. A. Schellenberger, J. Ripoll, D. Yessayan, E. Graves, A. Bogdanov Jr, L. Josephson, and R. Weissleder, “Visualization of antitumor treatment by means of fluorescence molecular tomography with an annexin V-Cy5.5 conjugate,” Proc. Natl. Acad. Sci. U.S.A. 101(33), 12294–12299 (2004). [CrossRef] [PubMed]

4. A. Becker, C. Hessenius, K. Licha, B. Ebert, U. Sukowski, W. Semmler, B. Wiedenmann, and C. Grötzinger, “Receptor-targeted optical imaging of tumors with near-infrared fluorescent ligands,” Nat. Biotechnol. 19(4), 327–331 (2001). [CrossRef] [PubMed]

5. V. Ntziachristos, C.-H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8(7), 757–761 (2002). [CrossRef] [PubMed]

6. R. Y. Tsien, “Building and breeding molecules to spy on cells and tumors,” FEBS Lett. 579(4), 927–932 (2005). [CrossRef] [PubMed]

7. Q. le Masne de Chermont, C. Chanéac, J. Seguin, F. Pellé, S. Maîtrejean, J.-P. Jolivet, D. Gourier, M. Bessodes, and D. Scherman, “Nanoprobes with near-infrared persistent luminescence for in vivo imaging,” Proc. Natl. Acad. Sci. U.S.A. 104(22), 9266–9271 (2007). [CrossRef] [PubMed]

8. F. A. Jaffer, P. Libby, and R. Weissleder, “Optical and multimodality molecular imaging: insights into atherosclerosis,” Arterioscler. Thromb. Vasc. Biol. 29(7), 1017–1024 (2009). [CrossRef] [PubMed]

9. J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011). [CrossRef] [PubMed]

10. N. Metropolis and S. Ulam, “The Monte Carlo method,” J. Am. Stat. Assoc. 44(247), 335–341 (1949). [CrossRef] [PubMed]

11. S. T. Flock, M. S. Patterson, B. C. Wilson, and D. R. Wyman, “Monte Carlo modeling of light propagation in highly scattering tissue--I: Model predictions and comparison with diffusion theory,” IEEE Trans. Biomed. Eng. 36(12), 1162–1168 (1989). [CrossRef] [PubMed]

12. L. Wang, S. L. Jacques, and L. Zheng, “MCML--Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef] [PubMed]

13. J. Barton, T. Pfefer, A. Welch, D. Smithies, J. Nelson, and M. Van Gemert, “Optical Monte Carlo modeling of a true portwine stain anatomy,” Opt. Express 2(9), 391–396 (1998). [CrossRef] [PubMed]

14. R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo Method (John Wiley & Sons, 2011).

15. I. Lux and L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations (CRC Press, 1991).

16. T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. 19(4), 879–888 (1992). [CrossRef] [PubMed]

17. S. E. Skipetrov and S. S. Chesnokov, “Analysis, by the Monte Carlo method, of the validity of the diffusion approximation in a study of dynamic multiple scattering of light in randomly inhomogeneous media,” Quantum Electron. 28(8), 733–737 (1998). [CrossRef]

18. D. Boas, J. Culver, J. Stott, and A. Dunn, “Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head,” Opt. Express 10(3), 159–170 (2002). [CrossRef] [PubMed]

19. A. J. Welch and R. Richards-Kortum, “Monte Carlo simulation of the propagation of fluorescent light,” in Laser-induced Interstitial Thermotherapy, G. J. Müller, and A. Roggan ed. (SPIE Press, 1995).

20. A. J. Welch, C. Gardner, R. Richards-Kortum, E. Chan, G. Criswell, J. Pfefer, and S. Warren, “Propagation of fluorescent light,” Lasers Surg. Med. 21(2), 166–178 (1997). [CrossRef] [PubMed]

21. T. J. Pfefer, K. T. Schomacker, M. N. Ediger, and N. S. Nishioka, “Light propagation in tissue during fluorescence spectroscopy with single-fiber probes,” IEEE J. Sel. Top. Quantum Electron. 7(6), 1004–1012 (2001). [CrossRef]

22. J. Swartling, A. Pifferi, A. M. Enejder, and S. Andersson-Engels, “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A 20(4), 714–727 (2003). [CrossRef] [PubMed]

23. D. Churmakov, I. Meglinski, S. A. Piletsky, and D. Greenhalgh, “Analysis of skin tissues spatial fluorescence distribution by the Monte Carlo simulation,” J. Phys. D Appl. Phys. 36(14), 1722–1728 (2003). [CrossRef]

24. D. Y. Churmakov, I. V. Meglinski, and D. A. Greenhalgh, “Amending of fluorescence sensor signal localization in human skin by matching of the refractive index,” J. Biomed. Opt. 9(2), 339–346 (2004). [CrossRef] [PubMed]

25. K. Vishwanath, B. Pogue, and M.-A. Mycek, “Quantitative fluorescence lifetime spectroscopy in turbid media: comparison of theoretical, experimental and computational methods,” Phys. Med. Biol. 47(18), 3387–3405 (2002). [CrossRef] [PubMed]

26. E. Péry, W. C. Blondel, C. Thomas, and F. Guillemin, “Monte Carlo modeling of multilayer phantoms with multiple fluorophores: simulation algorithm and experimental validation,” J. Biomed. Opt. 14(2), 024048 (2009). [CrossRef] [PubMed]

27. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 050902 (2013). [CrossRef] [PubMed]

28. A. T. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and B. J. Bacskai, “Time resolved fluorescence tomography of turbid media based on lifetime contrast,” Opt. Express 14(25), 12255–12270 (2006). [CrossRef] [PubMed]

29. A. Liebert, H. Wabnitz, H. Obrig, R. Erdmann, M. Möller, R. Macdonald, H. Rinneberg, A. Villringer, and J. Steinbrink, “Non-invasive detection of fluorescence from exogenous chromophores in the adult human brain,” Neuroimage 31(2), 600–608 (2006). [CrossRef] [PubMed]

30. A. Liebert, H. Wabnitz, R. Macdonald, and H. Rinneberg, “Monte Carlo simulations of time-resolved fluorescence excited in a layered turbid medium,” in Biomedical Topical Meeting, Technical Digest (Optical Society of America, 2006), paper ME7. [CrossRef]

31. A. Liebert, H. Wabnitz, N. Zołek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express 16(17), 13188–13202 (2008). [CrossRef] [PubMed]

32. J. Chen, V. Venugopal, and X. Intes, “Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates,” Biomed. Opt. Express 2(4), 871–886 (2011). [CrossRef] [PubMed]

33. A. T. Kumar, “Direct Monte Carlo computation of time-resolved fluorescence in heterogeneous turbid media,” Opt. Lett. 37(22), 4783–4785 (2012). [CrossRef] [PubMed]

34. C. K. Hayakawa and J. Spanier, Perturbation Monte Carlo Methods for the Solution of Inverse Problems (Springer, 2004).

35. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef] [PubMed]

36. N. Ren, J. Liang, X. Qu, J. Li, B. Lu, and J. Tian, “GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues,” Opt. Express 18(7), 6811–6823 (2010). [CrossRef] [PubMed]

37. G. Quan, H. Gong, Y. Deng, J. Fu, and Q. Luo, “Monte Carlo-based fluorescence molecular tomography reconstruction method accelerated by a cluster of graphic processing units,” J. Biomed. Opt. 16(2), 026018 (2011). [CrossRef] [PubMed]

38. F. Cai and S. He, “Using graphics processing units to accelerate perturbation Monte Carlo simulation in a turbid medium,” J. Biomed. Opt. 17(4), 040502 (2012). [CrossRef] [PubMed]

39. X. Yi, W. Chen, L. Wu, W. Ma, W. Zhang, J. Li, X. Wang, and F. Gao, “GPU-accelerated Monte-Carlo modeling for fluorescence propagation in turbid medium,” Proc. SPIE 8216, 82160U (2012). [CrossRef]

40. X. Yang, H. Gong, G. Quan, Y. Deng, and Q. Luo, “Combined system of fluorescence diffuse optical tomography and microcomputed tomography for small animal imaging,” Rev. Sci. Instrum. 81(5), 054304 (2010). [CrossRef] [PubMed]

41. J. W. Fu, X. Q. Yang, G. T. Quan, and H. Gong, “Fluorescence molecular tomography system for in vivo tumor imaging in small animals,” Chin. Opt. Lett. 8(11), 1075–1078 (2010). [CrossRef]

42. J. Fu, X. Yang, K. Wang, Q. Luo, and H. Gong, “A generic, geometric cocalibration method for a combined system of fluorescence molecular tomography and microcomputed tomography with arbitrarily shaped objects,” Med. Phys. 38(12), 6561–6570 (2011). [CrossRef] [PubMed]

43. X. Liu, D. Wang, F. Liu, and J. Bai, “Principal component analysis of dynamic fluorescence diffuse optical tomography images,” Opt. Express 18(6), 6300–6314 (2010). [PubMed]

44. Z. Zhang, W. Cao, H. Jin, J. F. Lovell, M. Yang, L. Ding, J. Chen, I. Corbin, Q. Luo, and G. Zheng, “Biomimetic nanocarrier for direct cytosolic drug delivery,” Angew. Chem. Int. Ed. Engl. 48(48), 9171–9175 (2009). [CrossRef] [PubMed]

45. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006). [CrossRef] [PubMed]

46. T. Correia, N. Ducros, C. D’Andrea, M. Schweiger, and S. Arridge, “Quantitative fluorescence diffuse optical tomography in the presence of heterogeneities,” Opt. Lett. 38(11), 1903–1905 (2013). [CrossRef] [PubMed]

47. J. Ye, C. Chi, Z. Xue, P. Wu, Y. An, H. Xu, S. Zhang, and J. Tian, “Fast and robust reconstruction for fluorescence molecular tomography via a sparsity adaptive subspace pursuit method,” Biomed. Opt. Express 5(2), 387–406 (2014). [CrossRef] [PubMed]

48. N. Chen, “Controlled Monte Carlo method for light propagation in tissue of semi-infinite geometry,” Appl. Opt. 46(10), 1597–1603 (2007). [CrossRef] [PubMed]

49. A. Sassaroli, “Fast perturbation Monte Carlo method for photon migration in heterogeneous turbid media,” Opt. Lett. 36(11), 2095–2097 (2011). [CrossRef] [PubMed]

Area	Excitation wavelength		Fluorescence wavelength
Area	$μ_{a}^{e x} (c m^{- 1})$	$μ_{s}^{e x} (c m^{- 1})$	$μ_{a}^{e m} (c m^{- 1})$	$μ_{s}^{e m} (c m^{- 1})$
Muscle	0.18	100	0.18	100
Lungs	0.33	200	0.33	200
Heart	0.19	76	0.19	76
Bone	0.02	180	0.02	180

Area	Excitation wavelength		Fluorescence wavelength
Area	$μ_{a}^{e x} (c m^{- 1})$	$μ_{s}^{e x} (c m^{- 1})$	$μ_{a}^{e m} (c m^{- 1})$	$μ_{s}^{e m} (c m^{- 1})$
Muscle	0.18	100	0.18	100
Lungs	0.33	200	0.33	200
Heart	0.19	76	0.19	76
Bone	0.02	180	0.02	180

Evaluation of path-history-based fluorescence Monte Carlo method for photon migration in heterogeneous media

Abstract

1. Introduction

2. Methods

2.1 Three types of fluorescence Monte Carlo method

2.1.1 sfMC method

2.1.2 dfMC method

2.1.3 pfMC method

2.2 Code and configuration of the GPU cluster scheme

2.2.1 Photon transport process

2.2.2 Computational settings

3. Results

3.1 Experimental verification of the sfMC

3.1.1 Experimental system

3.1.2 Cylindrical phantom

3.1.3 Experimental results

3.2 Comparison of the pfMC and dfMC methods

3.2.1 Statistical criteria

3.2.2 Cylindrical model

3.2.3 Efficiency comparison $\bar{e}$

3.2.4 Accuracy comparison

4. Discussion and conclusions

Acknowledgments

References and links

Cited By

Figures (8)

Tables (2)

Equations (22)

Optics Express

Abstract

1. Introduction

2. Methods

2.1 Three types of fluorescence Monte Carlo method

2.1.1 sfMC method

2.1.2 dfMC method

2.1.3 pfMC method

2.2 Code and configuration of the GPU cluster scheme

2.2.1 Photon transport process

2.2.2 Computational settings

3. Results

3.1 Experimental verification of the sfMC

3.1.1 Experimental system

3.1.2 Cylindrical phantom

3.1.3 Experimental results

3.2 Comparison of the pfMC and dfMC methods

3.2.1 Statistical criteria

3.2.2 Cylindrical model

3.2.3 Efficiency comparisone¯

3.2.4 Accuracy comparison

4. Discussion and conclusions

Acknowledgments

References and links

Cited By

Figures (8)

Tables (2)

Equations (22)

Optics Express

3.2.3 Efficiency comparison $\bar{e}$