Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Opt Lett. Author manuscript; available in PMC 2013 January 14.
Published in final edited form as:
Opt Lett. 2012 November 15; 37(22): 4783–4785.
PMCID: PMC3544305

Direct Monte Carlo computation of time-resolved fluorescence in heterogeneous turbid media


We show that a multiexponential model for time-resolved fluorescence allows the use of an absorption-perturbation Monte Carlo (MC) approach based on stored photon path histories. This enables the rapid fitting of fluorescence yield, lifetimes, and background tissue absorptions in complex heterogeneous media within a few seconds, without the need for temporal convolutions or MC recalculation of photon path lengths. We validate this method using simulations with both a slab and a heterogeneous model of the mouse head.

Several methods have been developed for the efficient Monte Carlo (MC) computation of light transport with intrinsic [13] and fluorescence [47] contrast in turbid tissue. For fluorescence, the most efficient approach is the adjoint fluorescence MC (aFMC, or reverse-emission MC [4]), which is based on convolving the absorption and emission Green’s function distributions with the fluorescence decay. In addition to a temporal double convolution, aFMC requires recomputation of the distributions for changes in tissue absorptions and a fine spatial and temporal binning for accuracy [4]. These aspects make aFMC computationally intensive for tomographic applications. Here we present an alternative fluorescence MC approach using a multiexponential model for time-resolved diffuse fluorescence (TRF) [8], which incorporates fluorescence lifetimes into Green’s functions implicitly through a reduced absorption. This allows a direct MC calculation of the entire TRF using perturbation MC [1,3,7], using photon path histories in Beer–Lambert factors. The method avoids the convolutions and spatial binning errors inherent in aFMC, while allowing rapid recalculation for changes in fluorophore or tissue absorption.

Consider a diffuse medium, with background absorption and scattering { μax,μsx} at excitation (λx) and { μae,μse} at the emission (λe) wavelengths, and N fluorophores, with lifetimes τn = 1/Γn, yield distributions ηn(r)=qnμfnx(r), quantum yields qn, and absorption coefficients μfnx. In the Born approximation, and assuming τn>(vμax,e(r))-1 (which is widely held [9]), the TRF for a source at rs and detector at rd on the surface can be expressed as a multiexponential sum [8],


with the time-dependent decay amplitudes given by


Here Wnx,e=Gnx(rs,r,t)[multiply sign in circle]Gne(r,rd,t) is a convolution of the Green’s functions Gnx,e, which are just the Green’s functions Gx, Ge for light propagation at λx and λe but evaluated with a reduced background absorption: μax,e(r)-Γn/v [8]. Equation (1) compactly expresses the TRF as a multiexponential sum, with time-dependent amplitudes reflecting the temporal evolution of the intrinsic diffuse response of the medium. Note that, besides the exponential decay factors in Eq. (1), the dependence on τn is solely through the reduced absorption. Thus, the fluorophore properties within An are entirely equivalent to an absorption, which suggests a perturbation MC approach [1] to rapidly recalculate the amplitudes An (and hence the full TRF) for any given set of lifetimes and fluorescence yields (and background tissue absorptions), using the stored path history of photon packets exiting the medium. The photon paths depend only on the scattering coefficient [1,4] and can be computed using a single MC simulation and stored in memory.

A direct application of the path history approach to evaluate Gnx and Gne in Eq. (2) would, however, be computationally intensive since it requires the storage of path histories from each source/detector to every voxel within the medium. Alternately, we first recognize that the quantity in the square brackets in Eq. (2) is simply the Born approximation term for the time-resolved (TR) photon fluence for an absorption perturbation equivalent to ηn(r), but with the sensitivities Wnx,e depending on two wavelengths. Suppose we ignore the wave-length dependence of Wnx,e, given the slow variation of tissue optical properties in the near-IR wavelengths and the typically narrow (<50 nm) separation between fluorophore absorption and emission spectra. We can then set μax=μae=μam and μsx=μse=μsm, where μam and μsm are the mean absorption and scattering between λx and λe. The quantity within the square brackets in Eq. (2) is then the difference of the exact diffuse photon fluences evaluated at { μam,μsm}, without and with an ab- sorption perturbation of μfnx:qn[Φ(μam-(Γn/v);t)-Φ(μam+μfnx-(Γn/v);t)] are dropped for simplicity]. Note that the reduced absorption is still applied for both terms, the significance of which will be discussed below.

Consider the imaging medium as divided into subregions [based for, e.g., on anatomically segmented images from magnetic resonance imaging (MRI) or computed tomography (CT)] indexed by “j.” The regionwise background and fluorophore absorptions are denoted as μamj and μfnxj. If Lkj denotes the path length of the kth photon in the jth region as computed by the MC simulation using the mean scattering coefficient, μsm, Eq. (2) can be approximated as (using Beer–Lambert law to calculate Φ)


where ‘k’ runs over the N(t) photons arriving between times 0 and t at detector rd. Note that the summation over k eliminates the t′ time integral in Eq. (2). Equation (3) is the central result of this Letter, and, together with Eq. (1), provides a fast analytical recipe for the calculation of TRF through a turbid medium of arbitrary shape and heterogeneity, provided the photon path lengths Lkj have been computed and stored in memory. [A similar approach was presented in [7] but without the use of a reduced absorption model as in Eqs. (1) and (2).]

We tested the accuracy of the history-FMC (hFMC) method [Eqs. (1) and (3)] by comparing it with the aFMC method. For the aFMC, Green’s functions from each rs and rd to all voxels were computed (108 photons, 512 time bins, 25 ps) and convolved with the fluorescence decay in the Fourier domain. For hFMC, 108 photons were launched at rs, and the total path lengths for each photon reaching the detector at rd were stored for the background ( Lk1) and inclusion ( Lk2), for a given μsmj and with μamj=0. Consider first the case of a homogeneous, isotropic (g = 0) diffusive slab (60 mm × 40 mm × 20 mm, 1 mm3 voxels), with a single source and 11 detectors (3 mm spacing), located on the surface, in the transmission geometry. A single 5 mm3 fluorescent inclusion is assigned near the center of the slab ( μf1x2=10-4) with no fluorescence in the background ( μf1x1=0). To quantify the influence of the key approxi- mation in Eq. (3), viz., equal optical properties at λx and λm, we recovered the fluorescence lifetime τ1 and absorption ( μf1x1 and μf1x2) by fitting the normalized TRF calculated using hFMC to that predicted using the aFMC approach [4] for a range of emission optical properties and with fixed excitation optical properties ( μax=0.01/mm,μsx=1/mm). The fitting was performed using the MATLAB “fminsearch” function. The error in the recovered lifetime (Fig. 1) was less than 2% for up to a 50% variation between μax and μae and decreased with longer lifetimes. The lifetime estimation error with respect to the scattering variations μsx-μae of up to 60% (with μse<=μsx as per Mie scattering theory) was also found to be well below 2% (not shown). For all the fits, the recovered fluorescence ratio of background to inclusion, μnfx1/μnfx2, was negligible (<10−5) as expected. The lifetime error for varying fluorophore absorption μf1x from 10−4 mm−1 up to the background absorption of 10−2 mm−1, with μf1e=0, was also estimated and found to be less than 1%. These errors are well below typical experimental uncertainties for lifetime estimation (>10%) [10]. Further, the ranges of variation of μa, μs used here are well above typical values, indicating the feasibility of the proposed model for practical applications. Note that, although the mean absorption, μam, was assumed to be known, it can also be included with the unknown fit parameters, μfnx and τn in Eq. (3), with a more complete tomographic data set.

Fig. 1
Error in lifetime recovered by fitting the hFMC approach [Eqs. (1),(3)] to the aFMC, for a slab model with a single fluorescent inclusion at the center (see text). The aFMC simulation assumed μax=0.01 and a varying μae between 0.005 ...

The accuracy of the entire TRF as predicted by hFMC was quantified as the rms error: E=[(1/N)k(Ukhis-Ukadj)2/Ukadj]1/2, which was less than 5% across the entire range of detectors and optical properties studied. Figure 2 shows a representative TRF for a single source–detector pair for the slab case, computed using the aFMC (solid blue) and hFMC (red circles). Also shown is the prediction of hFMC without the use of the reduced absorption in Eq. (3) (dotted black curve), which results in E > 15% and a lifetime error of >20% (Fig. 1, dashed-dotted curve), with the mismatch in the early time points of the TRF particularly noteworthy (Fig. 2). These results indicate the nontrivial role played by the reduced absorption in accurately predicting entire TRF. Equation (3) is thus also applicable for early photon tomography [11].

Fig. 2
(Color online) TRF curves calculated using the hFMC approach (red circles) with the aFMC (solid blue line) and the error, aFMC – hFMC (dash-dotted red line), for the same slab geometry as used in Fig. 1, with {μax,μae}={0.01, ...

We also tested the accuracy of hFMC for a heterogeneous, anisotropic medium with complex boundaries and multiple fluorophores. We used a publicly available segmented mouse atlas [12] [Fig. 3(a)], retaining only the head region (29 mm × 17 mm × 17 mm, 0.3 mm voxels) and assigned optical properties { μax,μsx,g} ranging from {0.001 mm−1, 1 mm−1, 0} (eyes) to {0.02 mm−1, 12.5 mm−1, 0.9} (brain). For each tissue segment, μae=0.8μax, and μse=0.8μsx. Four brain regions were assigned fluorophores with distinct lifetimes: cerebellum (0.6 ns), cerebrum (0.8 ns), striatum (1 ns) and rest of the brain (1.5 ns), all with μfx=10-4mm-1. Figure 3(b) shows representative TRFs computed using hFMC and aFMC (108 photons used for both) for a single source (S) and two detectors (D1 and D2) located on the mouse head. The two methods showed excellent agreement (E < 3%) across a range of detectors on the surface [open circles in Fig. 3(a)]. While aFMC required more than 3 h/source or detector using two quad-core Xeon 5472 CPUs, hFMC took less than 1 s for recalculation of the full TRF for any set of background absorption, fluorophore yield, and lifetimes.

Fig. 3
(Color online) (a) Heterogeneous absorption map of a digitized mouse, with four brain regions assigned fluorophores with distinct lifetimes as indicated. (b) Representative TRF curves predicted by aFMC (solid blue) and hFMC (red dots) for a source S(x) ...

The simulations presented here illustrate the feasibility of the hFMC method for rapid fitting of fluorescence yield and lifetime in heterogeneous turbid media. The increasing availability of high-resolution CT or MRI anatomical maps in conjunction with optical tomography makes hFMC particularly relevant for multimodality imaging. The hFMC can allow a rapid and accurate estimation of fluorescence in anatomical segments (e.g., fluorophore-labeled plaques or cancer cells within organs) rather than recovering full distributions, which are limited by the ill posedness of diffuse optical tomography. In addition, the approach can naturally take advantage of lifetime multiplexing [10] when lifetime sensitive probes are available, while also being applicable to other TR techniques such as early photon tomography [11]. In combination with hardware-accelerated MC implementations [13], the hFMC approach could be generalized to include fitting of tissue scattering, thereby making real-time recovery of in vivo fluorescence yield and lifetimes a possibility.


This work was supported by the National Institutes of Health grant R01 EB015325.


1. Hayakawa CK, Spanier J, Bevilacqua F, Dunn AK, You JS, Tromberg BJ, Venugopalan V. Opt Lett. 2001;26:1335. [PubMed]
2. Wang L, Jacques SL, Zheng L. Comput Methods Programs Biomed. 1995;47:131. [PubMed]
3. Boas DA, Culver JP, Stott JJ, Dunn AK. Opt Express. 2002;10:159. [PubMed]
4. Swartling J, Pifferi A, Enejder AMK, Andersson-Engels S. J Opt Soc Am A. 2003;20:714. [PubMed]
5. Liebert A, Wabnitz J, Zolek N, Macdonald R. Opt Express. 2008;16:13188. [PubMed]
6. Churmakov DY, Meglinski IV, Greenhalgh DA. J Biomed Opt. 2004;9:339. [PubMed]
7. Chen J, Venugopal V, Intes X. Biomed Opt Express. 2011;2:871. [PMC free article] [PubMed]
8. Kumar ATN, Raymond SB, Boverman G, Boas DA, Bacskai BJ. Opt Express. 2006;14:12255. [PMC free article] [PubMed]
9. Kumar ATN, Skoch J, Bacskai BJ, Boas DA, Dunn AK. Opt Lett. 2005;30:3347. [PubMed]
10. Raymond SB, Boas DA, Bacskai BJ, Kumar ATN. J Biomed Opt. 2010;15:046011. [PubMed]
11. Li Z, Niedre M. Biomed Opt Express. 2011;2:665. [PMC free article] [PubMed]
12. Dogdas B, Stout D, Chatziioannou AF, Leahy RM. Phys Med Biol. 2007;52:577. [PMC free article] [PubMed]
13. Fang Q, Boas DA. Opt Express. 2009;17:20178. [PMC free article] [PubMed]