Home | About | Journals | Submit | Contact Us | Français |

**|**Biomed Opt Express**|**v.2(4); 2011 April 1**|**PMC3072127

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Methodology
- 3. In silico study
- 4. Experiment
- 5. Discussion
- 6. Conclusion
- References and links

Authors

Related links

Biomed Opt Express. 2011 April 1; 2(4): 871–886.

Published online 2011 March 14. doi: 10.1364/BOE.2.000871

PMCID: PMC3072127

* Email: ude.ipr@xsetni

Received 2011 January 4; Revised 2011 February 25; Accepted 2011 February 26.

Copyright ©2011 Optical Society
of America

This is an open-access article distributed under the terms of the Creative Commons Attribution-Noncommercial-No Derivative Works 3.0 Unported License, which permits download and redistribution, provided that the original work is properly cited. This license restricts the article from being modified or used commercially.

This article has been cited by other articles in PMC.

Time-resolved fluorescence optical tomography allows 3-dimensional localization of multiple fluorophores based on lifetime contrast while providing a unique data set for improved resolution. However, to employ the full fluorescence time measurements, a light propagation model that accurately simulates weakly diffused and multiple scattered photons is required. In this article, we derive a computationally efficient Monte Carlo based method to compute time-gated fluorescence Jacobians for the simultaneous imaging of two fluorophores with lifetime contrast. The Monte Carlo based formulation is validated on a synthetic murine model simulating the uptake in the kidneys of two distinct fluorophores with lifetime contrast. Experimentally, the method is validated using capillaries filled with 2.5nmol of ICG and IRDye™800CW respectively embedded in a diffuse media mimicking the average optical properties of mice. Combining multiple time gates in one inverse problem allows the simultaneous reconstruction of multiple fluorophores with increased resolution and minimal crosstalk using the proposed formulation.

Optical methods offer the possibility to image numerous molecular targets with multiple distinct agents similarly to immunofluorescence microscopy and fluorescence cytometry. Current *in vivo* fluorescence multiplexing studies are mainly conducted in pre-clinical settings with Fluorescence Reflectance Imaging (FRI) [1, 2]. However, due to the limited information collected by FRI, the technique is unable to resolve the signal depth and hence to quantify it. Moreover, due to the predominance of scattering, planar imaging suffers from low resolution. For these reasons, FRI is mainly limited to superficial observations of subcutaneous flank xenograft models, surgically exposed organs or for intra-operative use, but is not appropriate to study advanced disease models in internal organs. For such applications, Fluorescence Molecular Tomography (FMT) techniques are required.

FMT acquires tomographic data sets to retrieve the 3D bio-distribution of the molecular probe used [3]. Spectral multiplexing of fluorophores in optical techniques requires the collection of dense spectral data sets for efficient unmixing of independent signals [4]. Such data sets are acquired sequentially one wavelength at a time, leading to relatively long acquisition times, especially for whole body applications. Alternatively, fluorophore unmixing can be performed based on lifetime contrasts [5]. By recording time dependent fluorescence signals, it is possible to distinguish and estimate the fractional contribution of multiple fluorophores with a monochromatic data set, allowing for fast acquisition protocols. Hence, lifetime based multiplexing is the most promising approach to simultaneously image multiple biomarkers for whole body applications *in vivo*.

Lifetime sensing is performed by employing Time Domain (TD) instruments. For optical tomography applications, image reconstructions based on TD fluorescence measurements have been primarily studied using the equivalent datatype in the frequency domain (FD), due to the simplicity of the relationship between the fluorescence lifetime and the measured phase in FD data [6, 7]. Recently, several studies have investigated the potential of performing FMT directly using TD derived datatypes such as moments [8], or time gates [9]. In the case of lifetime multiplexed studies, unmixing algorithms applied to the decaying portion of the TPSFs have been employed to recover the fractional contribution of the different fluorescence components and perform FMT based on unmixed time-independent signals [10]. However, such methodology suffers from the same limitation as the continuous wave (CW) technique, i.e., limited resolution and requires an unmixing algorithm where performances are not robust with low photon counts [11]. Conversely, the use of discrete time gates spanning the full TPSF as the data set for FMT should alleviate these drawbacks. It is well established that resolution of the optical reconstructions can be improved by using the rising portion of the TPSFs (termed early-gates) [12, 13] and such technique has been recently applied to FMT [14–16]. Moreover, the use of time gates allows for simultaneous reconstruction of multiple fluorophores based on lifetime contrast without the use of unmixing algorithms [17]. However, when considering time-resolved studies in small animals based on full TPSF, it is critical to employ an accurate light propagation model. Especially, for the early rising portion of the TPSF, a light propagation model that can accurately simulate the transition between minimally scattered photons and diffuse photons is required [18–22].

To date, analytical models, such as the radiative transport equation (RTE) and its approximation, the diffusion equation (DE), have been developed to solve the forward problem in FMT [21, 22]. However, the small finite volumes and their associated complex boundary conditions, the broad span of optical properties and potential low-scattering properties of certain organs may limit the validity of DE based models in small animals [23]. Most importantly, the early rising part of the measured TPSF which contains both minimally scattered photons and diffuse photons [18, 20, 24] is not accurately modeled by the DE [15, 19, 21, 22]. Alternatively, the Monte Carlo (MC) method is considered an accurate light propagation model without any of the limitations of the DE to simulate pulse propagations in diffuse media [23]. A number of studies focusing on Monte Carlo methods for fluorescence signal prediction have demonstrated its accuracy using synthetic and experimental data [25–27]. However, time-resolved Monte Carlo approaches have focused on spectroscopic applications but not tomographic ones. While inversion strategies based on Monte Carlo models in CW and FD have been reported [10, 28], no formulation for time-resolved fluorescence reconstruction based on Monte Carlo models has been yet reported.

In this work, we investigate the feasibility of performing time-resolved FMT of multiple fluorescence compounds based on Monte Carlo forward model and Time Gates. We derive a computationally efficient model that enables the calculation of the weight functions for both absorption perturbation and fluorophore distribution while simulating only the propagation of excitation photons in the tissue. We apply this new model to the simultaneous reconstruction of two fluorophores with lifetime contrast.We investigate for this purpose the information content imparted by single time gates and perform optical reconstructions based on multiple time gates to simultaneously image two fluorophores. The relative merits of different gate sampling strategies are compared and the validity of the model is experimentally established using a slab phantom with objects containing two commercial fluorophores.

The Monte Carlo method is capable of modeling light propagation in diffuse media over a broad spectral range, for small finite complex volume, for all the optical properties encountered in bio-tissue and for non-scattering or low-scattering regimes [23]. In this work, we follow the conventional forward-excitation and forward-emission model outlined by Welch et al. [29]. The model has been extended to the time domain by assuming that the emission occurs immediately after the absorption of the excitation photon. In this approach, two simulations of a complete set of photons are required: one is the propagation of the excitation photons, the other one is the isotropic emission and propagation of the emitted photons. The excitation simulation based on the optical properties at the excitation wavelength *λ _{x}* is used to calculate the statistical accumulation of the absorbed excitation light

$$\mathit{\eta}\left(r\right)=\frac{{\mu}_{af}^{x}\left(r\right)\mathrm{\Phi}}{{\mu}_{a}^{x}\left(r\right)}$$

(1)

where Φ is the quantum yield. The time-resolved fluorescence signal without delay caused by the lifetime of the fluorophore is then expressed as:

$${U}_{F}^{\prime}({r}_{s},{r}_{d},t\prime )={\int}_{\mathrm{\Omega}}d{r}^{3}W\prime ({r}_{s},{r}_{d},r,t\prime )\eta \left(r\right),$$

(2)

where the integration domain Ω is defined as the entire imaging volume, and the background weight function *W*′ is given by:

$$W\prime ({r}_{s},{r}_{d},r,t\prime )={\int}_{0}^{t\prime}dt\u2033A({r}_{s},r,t\u2033)E(r,{r}_{d},t\prime -t\u2033).$$

(3)

Assuming single exponential time-decay for the fluorophore, the fluorescence emission delay can be taken into account by convolving the temporal fluorescence signals with the decay e^{−t/τ} over time, where *τ* is the lifetime of the fluorophore. We can then write the fluorescence intensity measured at *r _{d}* and time

$${U}_{F}({r}_{s},{r}_{d},t)={\int}_{0}^{t}dt\prime {U}_{F}^{\prime}({r}_{s},{r}_{d},t\prime ){e}^{\frac{-(t-t\prime )}{\tau}}.$$

(4)

This general expression is used to compute the fluorescence forward model. In the case of multiple fluorophores, the different fluorescent components can be independently computed and summed to obtain the complex fluorescence time resolved measurements. In all the *in silico* studies herein, the synthetic fluorescence measurements were based on this approach.

The forward model described above can be casted into a tomography inverse problem to reconstruct the effective quantum yield. In order to compute the time-resolved Jacobians, or the spatial and temporal sensitivity maps with respect to *η*(*r*), Eq. (4) is written as:

$${U}_{F}({r}_{s},{r}_{d},t)={\int}_{\mathrm{\Omega}}d{r}^{3}W({r}_{s},{r}_{d},r,t)\eta \left(r\right),$$

(5)

where *W* is the lifetime based weight function defined as:

$$W({r}_{s},{r}_{d},r,t)={\int}_{0}^{t}dt\prime W\prime ({r}_{s},{r}_{d},r,t\prime ){e}^{-(t-t\prime )\tau},$$

(6)

Equations (3) and (6) give the general expression of the lifetime based weight matrix using Monte Carlo simulations. As Eq. (3) stands, the background weight function *W*′ can be calculated knowing *A*(*r _{s}*,

$${w}_{i}^{x}({r}_{s},r,{t}_{s,r})={w}_{i,0}\mathrm{exp}(-\underset{j=1}{\overset{{p}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{x}\left({r}_{j}\right)l\left({r}_{j}\right)),$$

(7)

where *w*_{i,0} is the initial weight of this photon, *r _{j}*(

$${A}_{i}({r}_{s},r,{t}_{s,r})={w}_{i}^{x}({r}_{s},r,{t}_{s,r})(1-\mathrm{exp}(-{\mu}_{a}^{x}\left(r\right){l}_{i}\left(r\right))).$$

(8)

The absorbed weight multiplied by *η*(*r*) is then converted to the initial weight of the fluorescence photon. We can calculate the final weight of the fluorescence photon detected at *r _{d}* as:

$${w}_{i}^{m}(r,{r}_{d},{t}_{r,d})={A}_{i}({r}_{s},r,{t}_{s,r})\eta \left(r\right)\mathrm{exp}(-\underset{j={p}_{i}+1}{\overset{{q}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{m}\left({r}_{j}\right){l}_{i}\left({r}_{j}\right)),$$

(9)

where *t*_{r,d} is the time that the photon passes from *r* to *r _{d}* and

$$w\prime ({r}_{s},{r}_{d},r,t)=\underset{i=1}{\overset{n}{\mathrm{\Sigma}}}{A}_{i}({r}_{s},r,t)\mathrm{exp}(-\underset{j={p}_{i}+1}{\overset{{q}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{m}\left({r}_{j}\right){l}_{i}\left({r}_{j}\right)).$$

(10)

Inserting Eqs. (7) and (8) to Eq. (10), then we have:

$$W\prime ({r}_{s},{r}_{d},r,t)=\underset{i=1}{\overset{n}{\mathrm{\Sigma}}}{w}_{i,0}\mathrm{exp}(-\underset{j=1}{\overset{{p}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{x}\left({r}_{j}\right){l}_{i}\left({r}_{j}\right))\times (1-\mathrm{exp}(-{\mu}_{a}^{x}\left(r\right){l}_{i}\left(r\right))\times \mathrm{exp}(-\underset{j={p}_{i}+1}{\overset{{q}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{m}\left({r}_{j}\right){l}_{i}\left({r}_{j}\right)).$$

(11)

According to Eq. (11), the background weight matrix for excitation source *r _{s}* and detector

$$W\prime ({r}_{s},{r}_{d},r,t)=\underset{i=1}{\overset{n}{\mathrm{\Sigma}}}{w}_{i}^{x}({r}_{s},{r}_{d},t)(1-\mathrm{exp}(-{\mu}_{a}^{x}\left(r\right){l}_{i}\left(r\right)),$$

(12)

where *w ^{x}_{i}*(

$${w}_{i}^{x}({r}_{s},{r}_{d},t)={w}_{i,0}\mathrm{exp}(-\underset{j=1}{\overset{{q}_{i}}{\mathrm{\Sigma}}}{\mu}_{a}^{x}\left({r}_{j}\right){l}_{i}\left({r}_{j}\right)).$$

(13)

If *μ ^{x}_{af}*(

$$W\prime ({r}_{s},{r}_{d},r,t)=\underset{i=1}{\overset{n}{\mathrm{\Sigma}}}{w}_{i}^{x}({r}_{s},{r}_{d},t){\mu}_{a}^{x}\left(r\right){l}_{i}\left(r\right).$$

(14)

Integrating Eq. (14) over time results in a CW reconstruction algorithm similar to the method proposed by Zhang et al. [28], that the photon paths weighted by their final fluence rate seen by the detector is the sensitivity function. Note also that the background weight matrix for fluorescence calculated by Eq. (14) is precisely the weight function for absorption perturbations [33]. This indicates that the absorption and fluorophore distribution can be estimated by using the same weight matrix calculated by one excitation Monte Carlo simulation. Moreover, it is noteworthy that if there are different fluorophore species existing in the region of interest, only one excitation simulation is required to obtain the weight matrices *W _{k}* (

To cast the inverse problem, we employed a Born formulation in which the emission field is normalized by the excitation field [35] by normalizing the experimental time-domain emission measurements to the CW excitation flux at the same position. We therefore have:

$$\frac{{M}^{m}({r}_{s},{r}_{d},t)}{{M}^{x}({r}_{s},{r}_{d})}=\frac{\alpha}{{U}^{x}({r}_{s},{r}_{d})}\underset{k=1}{\overset{{N}_{F}}{\mathrm{\Sigma}}}\int {d}^{3}r{W}_{k}({r}_{s},{r}_{d},r,t){\eta}_{k}\left(r\right),$$

(15)

where *α* incorporates unknown constants associated with wavelength-dependent gains and attenuations that can be measured once for every imaging system, *M ^{m}*(

The normalized fluorescence signals corresponding to different source-detector (S-D) pairs and time gates (left term in Eq. (15)) can be stacked up as a vector
$\gamma \in {R}^{{N}_{m}}$
(*N _{m}* is the total number of measurements, that is, the product of the number of S-D pairs and the number of time gates) where it can be used as part of an inverse problem. Similarly, the

$$\left[\begin{array}{c}{\gamma}_{1}\\ \vdots \\ {\gamma}_{{N}_{m}}\end{array}\right]=\left[\begin{array}{ccc}{\beta}_{1}{W}_{1,1}^{\mathrm{\Omega}}& \cdots & {\beta}_{1}{W}_{1,{N}_{F}}^{\mathrm{\Omega}}\\ \vdots & \ddots & \vdots \\ {\beta}_{{N}_{m}}{W}_{{N}_{m},1}^{\mathrm{\Omega}}& \cdots & {\beta}_{{N}_{m}}{W}_{{N}_{m},{N}_{F}}^{\mathrm{\Omega}}\end{array}\right]\left[\begin{array}{c}{\eta}_{1}^{\mathrm{\Omega}}\\ \vdots \\ {\eta}_{{N}_{F}}^{\mathrm{\Omega}}\end{array}\right],$$

(16)

where *W*^{Ω}_{s,k} is the weight function for the *s ^{th}*(

$$\mathrm{\Psi}={\parallel Ax-b\parallel}^{2}+{\parallel \lambda \left(r\right)x\parallel}^{2},$$

(17)

where *λ*(*r*) is a spatially variant regularization parameter to compensate for the spatial dependence of the contrast and resolution in the reconstruction [37].

In summary, our approach to reconstruct the effective quantum yield is schematically depicted in Fig. 1. The capital letters shown in the square brackets below refer to the individual steps in this scheme. [A] represents a time-resolved Monte Carlo simulation for excitation photons, following the WMC approach to analytically compute the light transport in tissue. This simulation can be adapted to different illumination and detection strategies.

To calculate the background weight matrix in [B], we directly add the weighted paths to the weight matrix along with the Monte Carlo simulation, instead of storing all the photon trajectories. This reduces the post processing time for weight matrix calculation as well as the storage space for the huge number of photon paths. Note that any of Eqs. (11), (12) and (14) can be used in [B]. It should be noted that Eq. (14), wherein the absorption coefficient at *λ _{x}* and

We first validated the proposed algorithm *in silico* in a model replicating small animal fluorescence imaging. The Monte Carlo simulations were performed on a mouse geometry with the kidneys and the skin shape extracted from a whole-body atlas [38] (Fig. 2(a)). The entire field of view consisted of 89 × 33 × 19 voxels with size 1mm^{3}. The optical properties were set to *μ ^{x}_{ab}* = 0.3cm

(a) The mouse model geometry and arrangement of sources (top) and detectors (bottom) used for the simulations. The fluorescent kidneys are marked as red and black, according to different lifetimes. (b) The illumination patterns (red: on, black: off).

In the Monte Carlo simulation, 10^{9} photons were consecutively launched for each pattern source on a massively parallel environment (blue gene, CCNI at RPI [39]) to calculate the excitation and fluorescence measurements as described in Sec. 2.1. Under the assumption that the optical properties are equal at *λ _{x}* and

The average calculation time for an excitation or emission MC simulation for one pattern source and all detectors and gates was less than 8 minutes (on 4096 nodes). Figure 3 shows typical simulated fluorescence signals and the background weight matrices calculated using Eq. (14) at different time gates.

In order to quantitatively assess the information content at single gates, we need to establish quantitative performance metrics. We investigated the quantification, crosstalk and resolution based on the reconstructions using one gate from 0.2ns to 2.9ns (opening of the gate) at an interval of 0.1ns. The quantification of the reconstructed *η*_{1} is defined as *Q*_{1} = max[*η*_{1}(Ω_{1})]/*H*_{1}(Ω_{1}), where Ω_{1} denotes the known location of the inclusions with lifetime *τ*_{1} = 0.5ns and *H*_{1} denotes the expected value of *η*_{1}. The crosstalk is defined as *X*_{1} = max[*η*_{2}(Ω_{1})]/max[*η*_{1}(Ω_{1})] to quantify the separability of the two inclusions with lifetime contrast. The quantification and yield crosstalk for the 1ns component was similarly evaluated. Note that *Q* and *X* are overestimated using the definitions based on only maximum value in the region of interest. The spatial resolution is measured as the full width at half maximum (FWHM) of the point spread functions (PSF). PSFs are created by simulating a perturbation at a single point. The simulated perturbation, *x _{sim}*, is a vector of zeros except for one target pixel at the center of the reconstructed volume with a value of one. We then generated simulated measurements

Information content evaluation based on single gate reconstructions in terms of (a) quantification, (b) crosstalk and (c) resolution.

From the plots of performance metrics, it is clear that for the reconstructed object with the shorter lifetime, the quantitative accuracy is decreasing with time while the crosstalk is increasing. Conversely, for the reconstructed object with the longer lifetime, an opposite behavior is observed with equivalent quantification and crosstalk around the maximum gate. It is noteworthy that the reconstruction using the latest gate (*t* = 2.8ns) results in the most accurate quantification and the least crosstalk for the compound with the longer lifetime, implying the necessity of the late gates for effective separation of the objects with lifetime contrast. It is impossible to achieve the best quantification and minimum crosstalk using only the early gates because of the comparable contribution of the two fluorophores at the early time points (Fig. 3). As expected, the steadily increasing FWHM of PSFs with time indicates that the early gates provide improved resolution when compared to the late gates [7, 15]. The better resolution at early gates is due to the narrower probing volume by photons detected at the early time points, as seen in Fig. 3(b). However, there is no significant reduction in the resolution of the reconstructed yield distributions for the two lifetimes when using gates later than the maximum gate. The reconstruction at the 10* _{th}* gate (

To utilize the various information carried by different gates, we investigated reconstructions using multiple gates simultaneously. In this regard, 5 sets of time gated measurements consisting of 5 gates each were considered. The number of gates was limited to 5 to reduce the memory burden as all the gates are simultaneously used in one inverse problem. The 5 different sets simulated in this study are illustrated in Fig. 5(a). Sets 1, 3, 5 have gates spaced at uniform intervals covering different time ranges, corresponding to the early to the maximum gate, full span, and the maximum to the end of the TPSFs. Sets 2 and 4 span the full time range, with non-uniform spacings exponentially scaled towards the rising and the decaying edges of the TPSFs respectively. This is a heuristic but straightforward approach to investigate the reconstruction performance using multiple combined gates. The reconstructions using combined gates of set 1 and set 5 are displayed in Fig. 5(b) and 5(c), respectively.

(a) Sets of time gates investigated in multi-gate reconstructions. (b) and (c) display the reconstructions using combined gates of set 1 and set 5, respectively.

Table 1 summarizes the quantitative metrics described in the previous section for the 5 sets investigated. It should be noted that the introduction of late gates improves the quantification of both fluorescence yields. Moreover, significantly less crosstalk than any of the single-gate reconstructions was achieved by using combined gates. The reconstruction using early gates (Set 1) has superior resolution compared to the reconstruction using late gates. However, the crosstalk for this case, is significantly higher than those applying late gates. For sets 2–4, the quantification and crosstalk are improving with the gate selection shifting to late gates, with similar performance in resolution due to the use of an early gate. It should be noted that the most accurate quantification and minimum crosstalk is obtained for both lifetimes when using the maximum and late gates (Set 5). However, without applying any early gates in reconstruction, the resolution is worsened by 20% compared to the other cases.

In order to ascertain the effect of assuming equal optical properties at excitation and emission wavelengths, a second set of time-resolved Jacobians were computed for the same set of gates using Eq. (11) (A–E). The reconstruction results when using sets are also given in Table 1. Comparing sets 1–5 and sets A–E, it is determined that the assumption of equal optical properties at excitation and emission wavelengths results in less than 5% quantification error for all three criteria.

The lifetime based Monte Carlo algorithm was experimentally validated by using a time-domain small animal molecular imaging system, which is schematically depicted in Fig. 6(a) [16, 40]. A tunable (690nm - 1020nm) Ti-Sapphire laser is used as the light source in the system. The laser source is expanded and projected on a digital micro-mirror device (DmD) based DLP board (Discovery 1100, Texas Instruments). 36 binary sliding-bar patterns similar to Fig. 2(b) are reflected on the programmable DmD chip and projected on the imaging chamber to provide an illumination area of 35mm × 20mm (matching the dimensions of a small animal torso). The transmitted signal is detected by a time-gated intensified CCD (ICCD) camera (Picostar HR, LaVision) with a time resolution of 40ps in the interest of a shorter acquisition time (24 minutes total using 36 patterns at excitation and emission wavelengths). A gate width of 300ps was selected to avoid temporal errors due to jitter when employing 200ps gates [41]. Photon measurements at the excitation and fluorescence wavelengths were acquired consecutively using bandpass filters. The incident fields for excitation and the fluorescence data for excitation were collected at 780nm and 832nm, respectively. The laser power was independently optimized for both the excitation and fluorescence measurements to maintain a maximum signal level of 4000 counts. Further, the acquired images are post-processed to generate 1mm × 1mm detector measurements by averaging over 7×7 pixels.

A polycarbonate slab phantom (Fig. 6(b)) that consisted of two tubes having an inner diameter of 1.5mm was used to experimentally assess the performance of the algorithm. To achieve lifetime contrast, tube I, II were filled with 2.5nmol of Indocyanine Green (*τ*_{1} = 0.45ns) and IRDye 800CW (*τ*_{2} = 0.8ns) dissolved in 180*µ*L ethanol, respectively. The effective quantum yields of the two tubes was externally calibrated and tube I was found to have 1.5 times the yield of tube II. The tank was filled with a mixture of Intralipid-20% (Sigma-Aldrich) and a water-soluble NIR Dye (Epolight 2717, Epolin) diluted in water to simulate the average optical properties (*μ _{a}* = 0.16cm

For the MC computation, the phantom was digitized into a 28 × 40 × 14 volume image with 1mm × 1mm × 1mm voxels. A background weight matrix as defined in Eq. (14) covering the field of view of 24 × 34 × 14 area, resulting in a total of 11,424 voxels, was calculated. The accuracy of the inverse model with experimental data is affected by the heterogeneous spatial distribution of the experimental pattern on the phantom. The illumination spatial heterogeneity is taken into account in the Monte Carlo simulation by scaling the initial photon weight to the illumination intensity at the injection position using experimental calibration patterns. The system impulse response functions (IRF) for all sources was convolved into the model before tomographic inversion, to match the model prediction with the experimental measurements. Since the lifetimes are known in advance in this controlled phantom experiment, the weight functions are generated using the shorter lifetime at *τ*_{1} = 0.45ns (ICG) and the longer lifetime at *τ*_{2} = 0.8ns (IRdye™) using Eq. (6). In *in vivo* applications, the lifetime of the dye within a control animal should be separately estimated due to the possible fluctuation of the fluorophore lifetime under different micro environments.

To achieve better separation between the two fluorophores while maintaining resolution, we selected the measurements at 20% of the rising edge, 75%, 35%, 20% and 15% of the decaying edge of the convolved TPSFs (mimicking gate set 4 in Sec. 3) to recover the fluorescence yield distribution based on the analysis of the simulated result. Gates with less than 400 photons (10% of the max for transmittance detector) were discarded in the reconstruction process to maintain an appropriate signal-to-noise ratio (SNR). This SNR led to a reduced size of the dataset by 8%. For comparison, we also reconstructed the fluorophore distribution using CW data at the same emission wavelength. The CW weight functions for the inversion were calculated by integrating the corresponding TD weight functions over time and the TPSFs were integrated over time to obtain the CW datatype.

The reconstruction results are shown in Fig. 7(a), similar as the 3D visualization of the *in silico* results. The two objects are separated with crosstalk *X*_{1} = 28.51% and *X*_{2} = 26.24%. The ratio of the mean reconstructed yields within the 50% isovolume of tube I and II was found to be 1.52. The average dimension of the reconstructed tubes was 34.5% larger along the *x* axis and 85% larger along the *z* axis (depth). This difference in resolution is due to the transmittance arrangement of the S-D pairs, which smears the resolution along the *z* axis. Hence, the diameter of the reconstructed isosurface is approximately 2 times bigger than the actual diameter due to the loss in the depth dimension, which can be attributed to the high scattering, and to the position of the two objects in the middle of the phantom, where the weight matrix has the least sensitivity. The CW reconstruction using one wavelength presented in Fig. 7(b), shows poor resolution (103% increase in depth) and no separation (100% crosstalk).

In this work we have proposed a new formulation to compute time-resolved fluorescence Jacobians for FMT. We validated this MC formulation for time-gated datatype to simultaneously image two fluorophores using *in silico* and experimental studies. The formulation is computationally efficient and offers three distinct advantages. First, once the excitation simulation is computed, it can be used to rapidly calculate the weight functions for multiple fluorophores with different lifetimes. Thus this method becomes attractive for lifetime multiplexed studies. Second, this scheme can be extended to multi-spectral imaging where the absorption coefficients can be modified in Eq. (14) to accommodate the changes in extinction coefficients at different wavelengths without the need of simulating another set of absorption coefficients. Third, using Eq. (11) the difference between *μ ^{x}_{a}* and

We validated the proposed formulation by reconstructing two different compounds simultaneously based on one excitation/emission wavelength pair without using an lifetime based unmixing algorithm on the diffuse measurements. We employed the time-gated datatype to retain the appeal of time-resolved measurements, improved resolution and lifetime contrast unmixing. While using early gates increases the resolution, the crosstalk between fluorophore is significant due to their equal contribution to the signal when used alone. The introduction of late gates in reconstruction is required to improve separation and quantification of fluorescent probes with lifetime contrast. The combination of early and late gates allows one to trade off the advantages and disadvantages of different time gates and thereby obtain enhanced resolution and crosstalk simultaneously. Employing 5 gates, the fluorescent inclusions were accurately resolved in the region of interest with minimal crosstalk. However, not all the individual metrics are optimal as an emphasis was put on reduction of crosstalk in this work due to the particular application considered. In the case of applications with only one compound, it will be probably possible to identify a set of gates that allow superior resolution and quantification [15]. Moreover, a reconstruction scheme, in which the spatial a priori information derived from an early gate reconstruction is incorporated in reconstruction based only on late gates (or even CW datatype for enhanced SNR), could provide superior results. Also, in the experiment conducted herein, the two compounds were spatially separated. The method proposed here however does allow the reconstruction of the concentration of two fluorophores that are co-localized. This cannot be achieved when employing methods that reconstruct both quantum yield and lifetime. And this provides a new method to monitor the multiple markers *in vivo*.

In this work we have developed a new formulation to perform time resolved FMT based on the Monte Carlo forward model. The formulation derived allows to compute fluorescence Jacobians within the fluorescence Born formulation in a computationally efficient manner. Based on standard assumption in the field of FMT, i.e. that the optical properties are similar at the excitation and emission wavelength, the time-resolved fluorescence Jacobians for multiple fluorophores can be computed from one excitation simulation.We have investigated the use of this new formulation in the context of simultaneous tomographic imaging of multiple fluorophores without using any unmixing algorithm. We considered the use of time-gated datatypes spanning the full TPSFs, including the early gates for which the diffusion equation is known to be inaccurate when compared to the MC. The *in silico* and experimental studies demonstrate that this new formulation can simultaneously reconstruct two fluorophores with improved resolution and reduced crosstalk when multiple gates are employed, whereas the two fluorophores cannot be resolved if only early gates are employed.

We gratefully acknowledge the technical support of the Computational Center for Nanotechnology Innovations (CCNI) at Rensselaer Polytechnic Institute.

1. Ntziachristos V., Ripoll J., Wang L. V., Weissleder R., “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005).10.1038/nbt1074 [PubMed] [Cross Ref]

2. Sevick-Muraca E. M., Rasmussen J. C., “Molecular imaging with optics: primer and case for near-infrared fluorescence techniques in personalized medicine,” J. Biomed. Opt. 13(4), 041303 (2008).10.1117/1.2953185 [PMC free article] [PubMed] [Cross Ref]

3. Hielscher A. H., “Optical tomographic imaging of small animals,” Curr. Opin. Biotechnol. 16(1), 79–88 (2005).10.1016/j.copbio.2005.01.002 [PubMed] [Cross Ref]

4. Mansfield J. R., “Distinguished photons: a review of *in vivo* spectral fluorescence imaging in small animals,” Curr. Pharm. Biotechnol. 11(6), 628–638 (2010).10.2174/138920110792246474 [PubMed] [Cross Ref]

5. Cubeddu R., Comelli D., D’Andrea C., Taroni P., Valentini G., “Time-resolved fluorescence imaging in biology and medicine,” J. Phys. D Appl. Phys. 35(9), R61–R76 (2002).10.1088/0022-3727/35/9/201 [Cross Ref]

6. Zhang L., Gao F., He H., Zhao H., “Three-dimensional scheme for time-domain fluorescence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express 16(10), 7214–7223 (2008).10.1364/OE.16.007214 [PubMed] [Cross Ref]

7. Soloviev V. Y., D’Andrea C., Valentini G., Cubeddu R., Arridge S. R., “Combined reconstruction of fluorescent and optical parameters using time-resolved data,” Appl. Opt. 48(1), 28–36 (2009).10.1364/AO.48.000028 [PubMed] [Cross Ref]

8. Lam S., Lesage F., Intes X., “Time domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express 13(7), 2263–2275 (2005).10.1364/OPEX.13.002263 [PubMed] [Cross Ref]

9. Kumar A. T. N., Raymond S. B., Boverman G., Boas D. A., Bacskai B. J., “Time resolved fluorescence tomography of turbid media based on lifetime contrast,” Opt. Express 14(25), 12255–12270 (2006).10.1364/OE.14.012255 [PMC free article] [PubMed] [Cross Ref]

10. Kumar A. T. N., Raymond S. B., Dunn A. K., Bacskai B. J., Boas D. A., “A time domain fluorescence tomography system for small animal imaging,” IEEE Trans. Med. Imaging 27(8), 1152–1163 (2008).10.1109/TMI.2008.918341 [PMC free article] [PubMed] [Cross Ref]

11. Köllner M., Wolfrum J., “How many photons are necessary for fluorescence-lifetime measurements?” Chem. Phys. Lett. 200(1-2), 199–204 (1992).10.1016/0009-2614(92)87068-Z [Cross Ref]

12. Liu F., Yoo K. M., Alfano R. R., “Ultrafast laser-pulse transmission and imaging through biological tissues,” Appl. Opt. 32(4), 554–558 (1993).10.1364/AO.32.000554 [PubMed] [Cross Ref]

13. Wu J., Perelman L., Dasari R. R., Feld M. S., “Fluorescence tomographic imaging in turbid media using early-arriving photons and Laplace transforms,” Proc. Natl. Acad. Sci. U.S.A. 94(16), 8783–8788 (1997).10.1073/pnas.94.16.8783 [PubMed] [Cross Ref]

14. Niedre M. J., de Kleine R. H., Aikawa E., Kirsch D. G., Weissleder R., Ntziachristos V., “Early photon tomography allows fluorescence detection of lung carcinomas and disease progression in mice *in vivo*,” Proc. Natl. Acad. Sci. U.S.A. 105(49), 19126–19131 (2008).10.1073/pnas.0804798105 [PubMed] [Cross Ref]

15. Niedre M., Ntziachristos V., “Comparison of fluorescence tomographic imaging in mice with early-arriving and quasi-continuous-wave photons,” Opt. Lett. 35(3), 369–371 (2010).10.1364/OL.35.000369 [PubMed] [Cross Ref]

16. Venugopal V., Chen J., Lesage F., Intes X., “Full-field time-resolved fluorescence tomography of small animals,” Opt. Lett. 35(19), 3189–3191 (2010).10.1364/OL.35.003189 [PubMed] [Cross Ref]

17. Suhling K., French P. M., Phillips D., “Time-resolved fluorescence microscopy,” Photochem. Photobiol. Sci. 4(1), 13–22 (2005).10.1039/b412924p [PubMed] [Cross Ref]

18. Mitra K., Kumar S., “Development and comparison of models for light-pulse transport through scattering-absorbing media,” Appl. Opt. 38(1), 188–196 (1999).10.1364/AO.38.000188 [PubMed] [Cross Ref]

19. Yaroslavsky I. V., Yaroslavsky A. N., Tuchin V. V., Schwarzmaier H.-J., “Effect of the scattering delay on time-dependent photon migration in turbid media,” Appl. Opt. 36(25), 6529–6538 (1997).10.1364/AO.36.006529 [PubMed] [Cross Ref]

20. Sakami M., Mitra K., Vo-Dinh T., “Analysis of short-pulse laser photon transport through tissues for optical tomography,” Opt. Lett. 27(5), 336–338 (2002).10.1364/OL.27.000336 [PubMed] [Cross Ref]

21. Xu M., Cai W., Lax M., Alfano R. R., “Photon migration in turbid media using a cumulant approximation to radiative transfer,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 65(6), 066609 (2002).10.1103/PhysRevE.65.066609 [PubMed] [Cross Ref]

22. Cai W., Xu M., Alfano R. R., “X. M, and R. Alfano, “Three-dimensional radiative transfer tomography for turbid media,” IEEE J. Sel. Top. Quantum Electron. 9(2), 189–198 (2003).10.1109/JSTQE.2003.813312 [Cross Ref]

23. J. C. Rasmussen, T. Pan, A. Joshi, T. Wareing, J. McGhee, and E. M. Sevick-Muraca, “Comparison of radiative transport, Monte Carlo, and diffusion forward models for small animal optical tomography,” in“2007 IEEE ISBI,” (2007), pp. 824–827.

24. Wang L., Liang X., Galland P., Ho P. P., Alfano R. R., “True scattering coefficients of turbid matter measured by early-time gating,” Opt. Lett. 20(8), 913–915 (1995).10.1364/OL.20.000913 [PubMed] [Cross Ref]

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

26. Swartling J., Pifferi A., Enejder A. M. K., Andersson-Engels S., “Accelerated Monte Carlo models to simulate fluorescence spectra from layered tissues,” J. Opt. Soc. Am. A 20(4), 714–727 (2003).10.1364/JOSAA.20.000714 [PubMed] [Cross Ref]

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

28. X. Zhang, C. Badea, M. Jacob, and G. A. Johnson, “Development of a noncontact 3-d fluorescence tomography system for small animal *in vivo* imaging,” SPIE **7171**, 71910D (2009). [PMC free article] [PubMed]

29. Welch A. J., Gardner C. M., Richards-Kortum R., Chan E., Criswell G., Pfefer J., Warren S., “Propagation of fluorescent light,” Lasers Surg. Med. 21(2), 166–178 (1997).10.1002/(SICI)1096-9101(1997)21:2<166::AID-LSM8>3.0.CO;2-O [PubMed] [Cross Ref]

30. Mourant J. R., Freyer J. P., Hielscher A. H., Eick A. A., Shen D., Johnson T. M., “Mechanisms of light scattering from biological cells relevant to noninvasive optical-tissue diagnostics,” Appl. Opt. 37(16), 3586–3593 (1998).10.1364/AO.37.003586 [PubMed] [Cross Ref]

31. Alerstam E., Andersson-Engels S., Svensson T., “White Monte Carlo for time-resolved photon migration,” J. Biomed. Opt. 13(4), 041304 (2008).10.1117/1.2950319 [PubMed] [Cross Ref]

32. A. Roggan, K. Dorschel, O. Minet, D. Wolff, and G. Muller, “The optical properties of biological tissue in the near infrared wavelength range : review and measurements,” in *Laser-induced Interstitial Thermotherapy* (SPIE, Bellingham WA, Etats-Unis, 1995), pp. 10–44.

33. Chen J., Intes X., “Time-gated perturbation Monte Carlo for whole body functional imaging in small animals,” Opt. Express 17(22), 19566–19579 (2009).10.1364/OE.17.019566 [PubMed] [Cross Ref]

34. Chen J., Venugopal V., Lesage F., Intes X., “Time-resolved diffuse optical tomography with patterned-light illumination and detection,” Opt. Lett. 35(13), 2121–2123 (2010).10.1364/OL.35.002121 [PubMed] [Cross Ref]

35. Soubret A., Ripoll J., Ntziachristos V., “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging 24(10), 1377–1386 (2005).10.1109/TMI.2005.857213 [PubMed] [Cross Ref]

36. Vinegoni C., Razansky D., Figueiredo J.-L., Nahrendorf M., Ntziachristos V., Weissleder R., “Normalized Born ratio for fluorescence optical projection tomography,” Opt. Lett. 34(3), 319–321 (2009).10.1364/OL.34.000319 [PMC free article] [PubMed] [Cross Ref]

37. Pogue B. W., McBride T. O., Prewitt J., Osterberg U. L., Paulsen K. D., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt. 38(13), 2950–2961 (1999).10.1364/AO.38.002950 [PubMed] [Cross Ref]

38. Dogdas B., Stout D., Chatziioannou A. F., Leahy R. M., “Digimouse: a 3D whole body mouse atlas from CT and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007).10.1088/0031-9155/52/3/003 [PMC free article] [PubMed] [Cross Ref]

39. “Computational Center for Nanotechnology Innovations (CCNI),” http://www.rpi.edu/research/ccni/

40. Venugopal V., Chen J., Intes X., “Development of an optical imaging platform for functional imaging of small animals using wide-field excitation,” Biomed. Opt. Express 1(1), 143–156 (2010).10.1364/BOE.1.000143 [PMC free article] [PubMed] [Cross Ref]

41. D’Andrea C., Comelli D., Pifferi A., Torricelli A., Valentini G., Cubeddu R., “Time-resolved optical imaging through turbid media using a fast data acquisition system based on a gated CCD camera,” J. Phys. D Appl. Phys. 36(14), 1675–1681 (2003).10.1088/0022-3727/36/14/304 [Cross Ref]

Articles from Biomedical Optics Express are provided here courtesy of **Optical Society of America**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |