Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates.

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.

1. Introduction

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.

2. Methodology

2.1. Time-resolved Monte Carlo forward model

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 λ is used to calculate the statistical accumulation of the absorbed excitation light A(r, r, t), defined as the absorbed photon weights by the fluorophore at position r at time t, illuminated by a source at r. The emission Monte Carlo simulation with the optical properties for the emission wavelength λ is computed to obtain the emission field E(r, r, t) at detector r and t. Here we define the total absorption coefficients as μ and μ at λ and λ, respectively. The total absorption coefficient at λ is defined as the sum of μ, the background absorption, and μ, the absorption coefficient contributed by the fluorophore, which is linearly related to the concentration of the fluorophore and the extinction coefficient. Then the effective quantum yield η(r) is defined as the probability of a photon to be emitted upon absorption of a photon by the total absorption coefficient:

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

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

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 over time, where τ is the lifetime of the fluorophore. We can then write the fluorescence intensity measured at r and time t for an impulsive excitation at r and t0 = 0 as:

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.

2.2. Calculation of the time-resolved Jacobians

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:

where W is the lifetime based weight function defined as:

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, r, t) and E(r, r, t) explicitly. This time convolution is extremely computationally intensive, and simulations with numerous photons at each position r in the region of interest are required to obtain statistically reliable calculation of E(r, r, t). Thus this algorithm is computationally inefficient especially in FMT applications where a large number of voxels are considered. However, a more manageable formulation can be derived by considering an assumption commonly employed in FMT, namely that the scattering coefficients are identical at λ and λ. In the NIR spectral range, the scattering coefficient is expected to vary less than 20% over the Stokes shift gap of organic fluorophores [30]. Furthermore, under the condition of isotropic scattering, it can be assumed that all fluorescence photons follow the trajectories of the excitation photons after generation [25]. Consider the i photon in an excitation simulation that propagates from a discrete source point r and then detected at the detector position r at time t. Following the White Monte Carlo (WMC) [31] approach, if the photon reaches r at time t, the weight of the i excitation photon w at r has decreased to:

where w is the initial weight of this photon, r(j = 1, ..., p) are the sub-regions that the photon consequently passes through from r to r, and l(r) is the path length that the photon passes at r. Note that the total number of p is photon dependent and the sub-regions for different photons are not necessarily the same. At r, the absorbed photon weight is given by:

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 as:

where t is the time that the photon passes from r to r and r(j = p + 1, ..., q) is the photon path from r to r. Note that sum of t and t is the total time t that the excitation photon travels from r till is detected, and that the sub-regions r(j = 1, ..., q) denote the total photon path. By summing up n detected emission photons for r and r, we can obtain U′ defined in Eq. (2). Comparing Eq. (9) to Eq. (2), we have:

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

According to Eq. (11), the background weight matrix for excitation source r and detector r at time t can be easily calculated using the photon paths and the absorption coefficients. This formulation should be employed in the case that the absorption coefficient at the excitation significantly differs from that at the emission wavelength, which might be the case for wavelengths below 700nm. However, the absorption optical spectra of bio-tissues is relatively flat in tissue in the NIR window [32]. Thus for simplicity, and following standard assumption in FMT, we assume that the absorption coefficients at the excitation and emission wavelengths are identical, that is, μ = μ. We therefore obtain:

where w(r, r, t) is the final weight of the i excitation photon from r and detected by r at t, defined as:

If μ(r)l(r) is small, which is the case in the voxelized geometry, from Taylor expansion, we can further simplify Eq. (12) to:

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 = 1, ..., N, where N is the number of species) for different lifetimes τ(k = 1, ..., N) according to Eq. (6). This allows for a fast and efficient computational implementation to simultaneously reconstruct multiple fluorophores based on one forward simulation. Note that for simplicity, all the above equations use r and r as discrete points, but they can be extended to any illumination or detection strategies applied to complex boundaries by collecting the photons generated or detected inside the areas. In this study, we employed an original broad field illumination strategy which is a new technique that we are currently developing [34], as described in Sec. 3 and 4.

2.3. Inverse problem

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:

where α incorporates unknown constants associated with wavelength-dependent gains and attenuations that can be measured once for every imaging system, M(r, r, t) is the total signal from all fluorophores with different lifetimes at the emission wavelength measured at the boundary position r and t excited from the source at r, M(r, r) and U(r, r) are the measured and simulated, respectively, total excitation flux measured by a photo-detector at r. This normalization efficiently mitigates the dependence of the detected fluorescent signal on the optical properties of the examined tissue [35, 36], thus we do not model the absorptive heterogeneities associated with the different organs in our synthetic murine model.

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 (N 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 k effective quantum yield for the entire imaging volume can be written as , where N is the total number of voxels in the region of interest. Overall, the forward problem can be expressed in a matrix form:

where WΩ is the weight function for the s(s = 1, ..., N) measurement and the k fluorophore species over the volume Ω, and β = α/U is the normalization factor for the s measurement and the corresponding weight function. The above equation has a general form A = b, then can be computed by solving this linear system of equations. Notice here the reconstructions for all fluorescent species provided herein are simultaneous based on monochromatic data and the measurements employed are direct measurements that are not derived from a pre-processed unmixing algorithm. Because of the ill-posedness of the inverse problem, the reconstruction method in fluorescence tomography usually sets up a least square optimization problem minimizing an objective function given by:

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

2.4. Reconstruction algorithm

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.

Fig. 1.

Block diagram summarizing the steps to perform fluorescence reconstruction.

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 λ and λ are assumed to be identical, can be used for computational efficiency. The lifetime-based weight matrices are then calculated by Eq. (6) in [C]. Knowing the experimental measurements, a conjugate gradient method is applied in [D] to minimize the objective function defined in Eq. (17) and obtain the fluorescence yield distribution. The fluorescence measurements in this study are either from an actual experiment or an emission Monte Carlo simulation at λ based on the absorption probability calculated in [A].

3. In silico study

3.1. Simulation setup

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 1mm3. The optical properties were set to μ = 0.3cm−1, μ′ = 25cm−1, μ = 0.36cm−1, μ′ = 20cm−1 over the entire body, to simulate the different optical properties of mouse tissues at different wavelengths in the NIR window [32]. The kidneys were considered to be labeled with two distinct fluorophores with lifetime of 0.5ns (black) and 1ns (red), respectively. The μ for both kidneys were set to 0.1 times of the background absorption coefficient μ, and the quantum yield was set to 1. 36 bar-shaped patterns were used as the source with each pattern illuminating half of the imaged surface along the x and y axes independently (Fig. 2(b)). A grid of 64 point detectors in transmission geometry were arranged covering a 31mm×27mm×18mm volume spanning the abdomen. The sources and detectors are projected to the surface along the z axis. The initial positions of the photons were randomly spanned (that is, uniformly distributed) over the illuminated area to accommodate the broad field illumination. The detectors had a 2mm separation along the y axis and a 2.2mm separation along the x axis with a radius of 1mm.

Fig. 2.

(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, 109 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 λ and λ (Eq. (14)), the time-resolved Jacobians for FMT were simultaneously computed by convolving the excitation field Jacobians by the lifetime decay. The temporal profile was recorded over a 6ns time window with 20ps resolution and a gate-width of 200ps (total of 120 gates). These parameters were selected to replicate the operating characteristics of our pre-clinical imaging platform [16] and correspond to the typical minimum gate width. Note that we do not explicitly include noise in the simulation data. However, the Poisson noise inherent in time-resolved studies also exists in Monte Carlo simulations. As the excitation measurements and fluorescence measurements are both directly simulated in this work (not computed through the multiplication of a Jacobian with a synthetic phantom), both measurements have an uncorrelated Poisson noise that replicates experimental conditions.

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.

Fig. 3.

(a) Representative simulated signals from the fluorophores with τ1 = 0.5ns (blue) and τ2 = 1ns (red), and the total fluorescence signal (black) used in the reconstructions. All the signals are generated using the mouse phantom and a single transmittance S-D pair as shown in (b) and (c). (b) Weight matrices for τ = 1ns at the 3 gates (gray) presented in (a). (c) The direct detector readings from the region of interest (black boxes) from the simulation on the mouse phantom. The bars in (b) and white boxes in (c) indicate the illumination pattern and the crosses in (b) and (c) indicate the detector used in (a).

3.2. Gate information content for fluorophore multiplexing

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 Q1 = max[η1(Ω1)]/H1(Ω1), where Ω1 denotes the known location of the inclusions with lifetime τ1 = 0.5ns and H1 denotes the expected value of η1. The crosstalk is defined as X1 = 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, 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 y = Ax and reconstructed the image. The cube root of the volume of voxels within half the peak reconstructed value is defined as FWHM. Plots of the quantification, crosstalk and FWHM of PSFs are shown in Fig. 4. Note that for all these reconstructions, the inverse problem size was identical and that the same iteration number (150) was used in the CG algorithm. This ensured consistency between the different reconstructions.

Fig. 4.

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 gate (t = 0.2ns) does not follow the trend due to the poor statistics of the weight matrix at the early gates, where the number of detected photons is insufficient to generate reliable statistics. It results in increased artifacts in the reconstructions when using very early gates. From the above analysis, we conclude that multiple fluorophores with lifetime contrast cannot be efficiently separated with minimal crosstalk if a single gate of the TPSF is used.

3.3. Lifetime multiplexing tomography based on multiple gates

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.

Fig. 5.

(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.

Table 1.

Quantitative Comparison of the Multi-Gate Reconstructions

Set	Quantif. (%)		Crosstalk (%)		FWHM (mm)
	τ1	τ2	τ1	τ2	τ1	τ2
1	44.53	45.65	39.74	50.43	2.33	2.44
2	45.89	49.22	35.69	36.50	2.52	2.64
3	46.98	52.19	32.25	33.42	2.69	2.78
4	49.04	53.96	29.20	30.66	2.74	2.85
5	49.77	63.05	24.92	24.10	3.21	3.16
A	47.64	45.72	37.19	49.19	2.19	2.32
B	48.59	50.21	33.34	33.98	2.43	2.46
C	49.35	53.40	31.01	31.45	2.63	2.67
D	50.47	54.39	27.34	30.05	2.72	2.76
E	55.76	66.45	21.77	23.19	3.17	3.11

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.

4. Experiment

4.1. Diffuse optical tomography system

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.

Fig. 6.

(a) Schematic of the TD fluorescence tomography system. LS: laser; PC: power control; FBC: fiber-beam coupler; BE: beam expander; DMD: DmD-based DLP chip; IL: imaging lens; IC: imaging chamber; CCD: intensified CCD camera; HRI: high rate imager; TDU: trigger delay unit; TG: trigger generator. (b) Phantom setup. The red box shows the pattern illumination area. All dimensions are in mm.

4.2. Phantom setup

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 (μ = 0.16cm−1 and μ′ = 17cm−1) of murine models.

4.3. Experimental reconstruction

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 X1 = 28.51% and X2 = 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).

Fig. 7.

Experiment results using (a) TD data and (b) CW data. The images are normalized to the maximum of the reconstruction.

5. Discussion

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 μ and μ can be taken into account, noticing that there are no significant increase of computational burden using Eqs. (11) and (14) thanks to the WMC method [31]. Such flexibility in incorporating differences in spectral absorption allows its use for compounds with relatively large Stokes shift (e.g. quantum dots). Here we investigated a 20% change in the in silico study and an improvement less than 5% was achieved using the corrected μ, implying the proposed approach is robust for reconstructions in the NIR window. On the other hand, the scattering coefficient cannot be rescaled like the absorption coefficient due to the transmittance geometry employed in this work [31]. While an inaccurate estimation of the scattering coefficient will result in errors in the reconstruction, μ varies more slowly (seldom more than 20%) over the spectral range in the NIR window [30]. Previous studies have reported that 20% change in μ will result in less than 15% error in diffuse optical tomography [42] and similar results are expected in FMT.

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.

6. Conclusion

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.