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

**|**HHS Author Manuscripts**|**PMC2700299

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Theoretical development
- 3. Singular value analysis of the time domain weight function
- 4. Tomography using direct TD and asymptotic approaches
- 5. Conclusions
- References and links

Authors

Related links

Opt Express. Author manuscript; available in PMC 2009 June 23.

Published in final edited form as:

PMCID: PMC2700299

NIHMSID: NIHMS86126

Anand T. N. Kumar

Department of Neurology and Department of Radiology, Massachusetts General Hospital, Harvard Medical School, Charlestown MA 02129 ; Email: ude.dravrah.hgm.rmn@ramukna

Scott B. Raymond

Harvard-MIT Division of Health Science and Technology, Boston MA, 02115

Gregory Boverman

Department of Biomedical Engineering, Rensselaer Polytechnic Institute, Troy NY, 12180

David A. Boas

Athinoula A Martinos center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown MA, 02129

Alzheimer’s Disease Research Unit, Department of Neurology, Massachusetts General Hospital, Charlestown MA, 02129

The publisher's final edited version of this article is available at Opt Express

See other articles in PMC that cite the published article.

A general linear model for time domain (TD) fluorescence tomography is presented that allows a lifetime-based analysis of the entire temporal fluorescence response from a turbid medium. Simulations are used to show that TD fluorescence tomography is optimally performed using two complementary approaches: A direct TD analysis of a few time points near the peak of the temporal response, which provides superior resolution; and an asymptotic multi-exponential analysis based tomography of the decay portion of the temporal response, which provides accurate localization of yield distributions for various lifetime components present in the imaging medium. These results indicate the potential of TD technology for biomedical imaging with lifetime sensitive targeted probes, and provide useful guidelines for an optimal approach to fluorescence tomography with TD data.

The development of disease-specific fluorescent markers and genomic reporters has prompted concurrent advances in optical tomography techniques for the non-invasive diagnosis of disease in a living animal or human subject [1, 2]. The most common optical tomographic techniques for fluorescence are based on continuous wave (CW) excitation [3] and frequency modulated (FD) excitation [4, 5]. Another measurement mode is in the time domain (TD) [7–15], where the excitation is performed using short laser pulses in conjunction with time resolved detection. Fluorescence lifetime reconstructions with turbid media have been discussed in previous works in conjunction with FD [4, 5] and TD [6, 11–14] measurements (CW measurements are incapable of distinguishing fluorescence lifetime from yield). A single TD measurement with a short laser pulse implicitly contains all modulation frequencies, and hence provides the most complete optical information. Moreover, the surface decay data from TD measurements can directly reveal the intrinsic fluorophore lifetime, without the need for reconstructions [11–13]. This feature could be of tremendous importance given the potential development of lifetime sensitive probes for in-vivo applications and the sensitivity of fluorescence lifetime to factors affecting local tissue environment such as pH, viscosity, oxygen concentration and tissue autofluorescence, in addition to molecular interactions such as Forster resonance energy transfer (FRET) [17]. Fluorescence lifetime imaging (FLIM) is already a well established microscopy technique that is used to probe lifetime contrast in thin tissue sections [18, 19].

Image reconstruction algorithms for TD fluorescence measurements have so far been primarily based on derived, or transformed, data types, such as the Laplace transform [8, 23], Mellin transform or moments [15] and the Fourier transform [27, 28]. One advantage of working in a transformed space of the TD data lies in the simplicity of the relationship between the fluorescence lifetime and the measured phase, as compared to the non-linear dependence on fluorophore lifetime (through the exponential decay factor) in the direct TD case. For instance, the phase measured in FD experiments is linearly related to the lifetime distribution [4]. Nevertheless, the lifetime is still in the form of a distribution, which can only be recovered using tomographic reconstructions to remove the contribution to the phase from diffusive propagation. Also, the measured phase at a certain modulation frequency (or imaginary frequency in the Laplace case) is an admixture of contributions due to all the lifetimes present in the medium. On the other hand, Laplace transforms have been applied to selectively reconstruct the early rise portion of the temporal response curve, since early arriving photons undergo minimal scattering, and are largely unaffected by the long lived fluorophore lifetimes [8].

We recently demonstrated experimentally [13] that analyzing the asymptotic decay portion of the diffuse fluorescence temporal response (DFTR) can by itself have distinct advantages: The yield distribution(s) for multiple lifetime component(s) present within the medium can be localized separately using the surface decay amplitudes extracted from multi-exponential fits. In what follows we will simply label this the “asymptotic” approach. The asymptotic approach reduces a cumbersome analysis of a large temporal data set in the decay portion of the DFTR into a multi-exponential curve fitting followed by simple CW reconstructions. This approach can be viewed as an application of the *inverse* Laplace transform (which is equivalent to multi-exponential fitting for a few discrete lifetimes) to reconstruct the decay portion of the DFTR. However, the restriction of this method to the decay portion excludes the information from the earlier portion of the DFTR, which is characterized by a better signal-to-noise (SNR) ratio, and may also contain useful spatial information. It is thus imperative to seek an approach that incorporates the rising and peak portions of the DFTR data into the analysis. In this work, we develop a theoretical formalism that allows a lifetime-based separation of fluorescence yield distributions using the entire TD data. In order to evaluate the optimal choice of temporal measurements for tomography using the direct TD approach, we use a singular value decomposition analysis [29, 30] of the TD weight matrix. To our knowledge, this optimization has not been carried out previously for TD fluorescence measurements, although optimization of multi-frequency data in FD has been studied previously [27, 28]. The relative merits of the optimized direct TD approach and the asymptotic approach are then compared using simulations and the advantages of TD data over CW, in the presence of lifetime contrast, are demonstrated.

This paper consists of three central parts. In Section 2, key integral expressions are presented that enable lifetime-specific tomographic reconstructions of the entire DFTR, along with a discussion of the conditions when the results are valid. In Section 3, the optimization of temporal measurements for fluorescence tomography is addressed numerically, using a SVD analysis. In Section 4, the theoretical formalism developed and the results of the SVD analysis are applied to simulated noisy data to more specifically determine the imaging performance of the rise and decay portions of the DFTR.

In this section, we revisit the basic expressions involved in TD fluorescence tomography, and present a simplified expression under the specific condition of long lifetimes. Consider a turbid imaging medium embedded with fluorophores, described by yield and lifetime distributions (at one wavelength) η (**r**) and τ (**r**) = 1/Γ(**r**), respectively. For tomography, optical sources and detectors are arranged on the surface of the imaging medium. The fluorescence intensity measured at a detector point **r*** _{d}* at time

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

(1)

where the weight function, or sensitivity function is given by

$$W({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t)={\int}_{0}^{t}{\mathit{dt}}^{\prime}{\int}_{0}^{{t}^{\prime}}{\mathit{dt}}^{\prime \prime}{G}^{m}({\mathbf{r}}_{d}-\mathbf{r},t-{t}^{\prime}){e}^{-\Gamma \left(\mathbf{r}\right)({t}^{\prime}-{t}^{\prime \prime})}{G}^{x}(\mathbf{r}-{\mathbf{r}}_{s},{t}^{\prime \prime}),$$

(2)

with *G ^{x}* and

As it stands, the TD fluorescence weight function in Eq. (2) is a double convolution in time and is computationally intensive, especially for a tomographic measurement setup with a large number of sources and detectors. But a closer inspection reveals that Eq. (2) can be rewritten in a more manageable form. First, we define a background weight function as:

$${W}^{B}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},{t}^{\prime})={\int}_{0}^{{t}^{\prime}}{\mathit{dt}}^{\prime \prime}{G}^{x}({\mathbf{r}}_{s},\mathbf{r},{t}^{\prime}-{t}^{\prime \prime}){G}^{m}(\mathbf{r},{\mathbf{r}}_{d},{t}^{\prime \prime}),$$

(3)

which depends only on the medium optical properties and reduces to the weight function for an absorption perturbation when the excitation and emission wavelengths coincide. Using the commutativity of the convolution operation, we can now re-write Eq. (2) as,

$$W({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t)={\int}_{0}^{t}{\mathit{dt}}^{\prime}{W}^{B}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},{t}^{\prime}){e}^{-\Gamma \left(\mathbf{r}\right)(t-{t}^{\prime})}.$$

(4)

Since *W ^{B}* can be pre-calculated with a knowledge of background optical properties, the advantage of Eq. (4) over Eq. (2) is that only a single time integral is left for the tomographic recovery of the yield and lifetime distributions [20]. A more useful form of this expression is realized if the fluorophores within the medium are described as multiple distributions,

$${W}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t)={\int}_{0}^{t}{\mathit{dt}}^{\prime}{W}^{B}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},{t}^{\prime}){e}^{-{\Gamma}_{n}(t-{t}^{\prime})},$$

(5)

so that the total fluorescence signal is expressed as

$${U}_{F}({\mathbf{r}}_{s},{\mathbf{r}}_{d},t)=\sum _{n}{\int}_{\Omega}{d}^{3}{\mathit{rW}}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t){\eta}_{n}\left(\mathbf{r}\right).$$

(6)

If it is further assumed that the lifetimes are longer than the absorption timescale, i.e., *τ _{n}* >

$$\begin{array}{cc}\hfill (\mathbf{s}\cdot \nabla & +\frac{1}{v}\frac{\partial}{\partial t}+{\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)+{\mu}_{s}^{(x,m)}\left(\mathbf{r}\right)){G}^{(x,m)}(\mathbf{r},\mathbf{s},t)\hfill \\ \hfill & ={\mu}_{s}^{(x,m)}\left(\mathbf{r}\right){\int}_{\Omega}\Theta (\mathbf{s},{\mathbf{s}}^{\prime}){G}^{(x,m)}(\mathbf{r},{\mathbf{s}}^{\prime},t)d{\mathbf{s}}^{\prime},\hfill \end{array}$$

(7)

where Θ(**s,s’**) is the scattering phase function. The source terms are dropped from the excitation RTE on the basis that the fluence is calculated away from the source location, and similarly from the emission RTE, given that only a single fluorophore emission event is considered in accordance with the Born approximation initially made in Eq. 2. (Multiple absorption of the excitation and emission light by the fluorophore is still incorporated in the total absorption at these wavelengths viz., ${\mu}_{a}^{x}$ and ${\mu}_{a}^{x}$.) Suppose now that we write (dropping the angular dependence for simplicity):

$${G}^{(x,m)}(\mathbf{r},t)={G}_{0}^{(x,m)}(\mathbf{r},t){e}^{-v{\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)t},$$

(8)

it can be verified by substituting the above solution into Eq. (7) that the functions ${G}_{0}^{(x,m)}$ are dependent only on the *gradient* of the absorption coefficient, $\nabla {\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)$, and independent of ${\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)$ itself. Thus, ${G}_{0}^{(x,m)}$ are invariant to constant shifts in the absorption. If we define the Green’s functions ${G}_{n}^{(x,m)}$ evaluated using a reduced background medium absorption, ${\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)-{\Gamma}_{n}\u2215v$, which is positive under the long lifetime condition viz., ${\Gamma}_{n}<v{\mu}_{a}^{(x,m)}\left(\mathbf{r}\right)$, it is then easily verified using Eq. (8) that

$${G}_{n}^{(x,m)}(\mathbf{r},t)={G}^{(x,m)}(\mathbf{r},t){\mid}_{{\mu}_{a}-{\Gamma}_{n}\u2215v}={G}^{(x,m)}(\mathbf{r},t){\mid}_{{\mu}_{a}}{e}^{{\Gamma}_{{n}^{t}}}.$$

(9)

Now, writing e^{Γ nt’} = *e*^{Γn(t’—t”)eΓnt”} , we can use Eq. (9) in Eq. (5) to show that:

$${W}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t)={e}^{-{\Gamma}_{{n}^{t}}}{\int}_{0}^{t}{\mathit{dt}}^{\prime}{W}_{n}^{B}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t{}^{\prime}),$$

(10)

where ${W}_{n}^{B}$ is given by Eq. (3) but with the reduced absorption Green’s functions ${G}_{n}^{x,m}$.

The form of the weight function in Eq. (10) allows the fluorescence signal to be expressed as a multi-exponential sum, analogous to fluorescence lifetime imaging [18] (FLIM):

$${U}_{F}({\mathbf{r}}_{s},{\mathbf{r}}_{d},t)=\sum _{n}{A}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},t){e}^{-{\Gamma}_{{n}^{t}}},$$

(11)

where the decay amplitudes *A _{n}* depend on time, in addition to the source and detector locations. The amplitudes define a linear inversion problem for the yield distributions of each lifetime component:

$${A}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},t)=\int {d}^{3}r\left[{\int}_{0}^{t}d{t}^{\prime}{W}_{n}^{B}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},{t}^{\prime})\right]{\eta}_{n}\left(\mathbf{r}\right).$$

(12)

For fluorophores with lifetimes comparable to optical diffusion time scales in biological media (≈ nanoseconds), *A _{n}(t)* has a non-trivial time evolution that is determined by the size and optical properties of the imaging medium. In Figure 1, the temporal evolution of

From Eq. (10), it is clear that the weight function for each lifetime component is an average over a time *t* of the background sensitivity function ${W}_{n}^{B}$. Let *τ _{D}* denote the timescale for the evolution of ${W}_{n}^{B}$, which will depend on absorption, scattering and medium boundaries (see section 2.1 below). For

$$\underset{t>{\tau}_{D}}{\mathrm{lim}}{W}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r},t)\to {e}^{-{\Gamma}_{{n}^{t}}}{\stackrel{\u2012}{W}}_{n}({\mathbf{r}}_{s},{\mathbf{r}}_{d},\mathbf{r}).$$

(13)

Therefore, Eqs. (11-12) along with Eq. (3) constitute a TD generalization of asymptotic lifetime-based tomography, that includes the early arriving photons in addition to the asymptotic decay portion corresponding to the late arriving fluorescent photons. Note from Fig. 1 that the amplitude of the asymptotic fluorescence decay (shown as dashed red lines) equals the long time value of *A(t)*, which is related to the CW sensitivity function .

The results presented in this section can be summarized as follows. With the lifetimes calculated from the asymptotic decay of the TD signal, Eq. (11) can be used to reconstruct the yield distribution for each lifetime using TD data. Since the amplitude for each lifetime is in general time dependent and cannot be separated, the reconstruction is performed directly on the measured data:

$$\left(\begin{array}{c}\hfill .\hfill \\ \hfill {U}_{F}^{j}\left({t}_{k}\right)\hfill \\ \hfill {U}_{F}^{j}\left({t}_{k+1}\right)\hfill \\ \hfill .\hfill \end{array}\right)=\left(\begin{array}{cccc}\hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill {W}_{n}^{j}\left({t}_{k}\right)\hfill & \hfill {W}_{n+1}^{j}\left({t}_{k}\right)\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill {W}_{n}^{j}\left({t}_{k+1}\right)\hfill & \hfill {W}_{n+1}^{j}\left({t}_{k+1}\right)\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \end{array}\right)\left(\begin{array}{c}\hfill .\hfill \\ \hfill {\eta}_{n}\left(\mathbf{r}\right)\hfill \\ \hfill {\eta}_{n+1}\left(\mathbf{r}\right)\hfill \\ \hfill .\hfill \end{array}\right),$$

(14)

where for simplicity of notation we have dropped the source-detector co-ordinates and have instead used a single superscript *j* to denote measurement index that labels each source-detector (S-D) pair. For times longer than *τ _{D}*, the decay amplitude becomes time-independent, so that the amplitude for each lifetime component can be recovered asymptotically using multi-exponential fits. These amplitudes constitute a derived data set (inverse Laplace transform) for the inversion of the yield distributions:

$$\left(\begin{array}{c}\hfill .\hfill \\ \hfill {A}_{n}^{j}\hfill \\ \hfill {A}_{n+1}^{j}\hfill \\ \hfill .\hfill \end{array}\right)=\left(\begin{array}{cccc}\hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill {\stackrel{\u2012}{W}}_{n}^{j}\hfill & \hfill 0\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill 0\hfill & \hfill {\stackrel{\u2012}{W}}_{n+1}^{j}\hfill & \hfill .\hfill \\ \hfill .\hfill & \hfill .\hfill & \hfill .\hfill & \hfill .\hfill \end{array}\right)\left(\begin{array}{c}\hfill .\hfill \\ \hfill {\eta}_{n}\left(\mathbf{r}\right)\hfill \\ \hfill {\eta}_{n+1}\left(\mathbf{r}\right)\hfill \\ \hfill .\hfill \end{array}\right).$$

(15)

The key difference between Eqs. (14) and (15) is that the asymptotic weight function in Eq. (15) is block diagonal, whereas the TD weight function in Eq. (14) has off-diagonal terms. Thus, the direct TD reconstruction will be characterized by significantly more cross-talk between the lifetime distributions than the asymptotic reconstruction. At least two questions immediately arise, related to the practical application of the above results. Firstly, what is the optimal choice of time points for the TD reconstruction using Eq. (14)? Secondly, what are the relative merits of the direct TD and asymptotic approaches? We will address these questions in Sections 3 and 4.

A basis for the theoretical work presented in this paper is the direct estimation of the intrinsic fluorophore lifetimes from the decay of TD fluorescence signals. There are two different time scales involved in determining whether fluorophore lifetimes are revealed in the measured decay on the surface of the turbid medium. First is the intrinsic absorption time scale *τ _{a}* = (

The strong condition is easily satisfied for nanosecond lifetime fluorophores in biomedical applications (μ* _{a}* > 0.1

In this section, we will present an optimization study of the number and location of time points for a direct TD reconstruction. The general optimization of source-detector (S-D) pairs and time points is a complex multi-dimensional problem since each S-D pair could ideally be associated with a different time gate. It is, however, reasonable to view the temporal points and S-D arrangements as independent dimensions in the optimization, since for biomedical imaging applications, the length scales involved are not too large and the correspondingly small variations in the temporal response along the measurement surface can be assumed not to affect the results in a significant way. Therefore, in this paper, we consider a fixed S-D geometry and focus on optimizing the temporal measurements for fluorescent tomography. The optimization of S-D configuration has been discussed in previous works for CW fluorescence [29] and nonfluorescent [30] tomography, using a singular value decomposition (SVD) analysis. Here, we will apply SVD to the TD weight matrix *W _{n}* defined in Eq. (10) for optimization of the time points within the DFTR. SVD of a matrix

The weight matrix as defined in Eq. (10) was simulated for a diffusive slab medium of thickness 2cm in the transmission geometry, with a S-D arrangement as shown in Fig. 2, with 21 sources and 29 detectors arranged in a honeycomb pattern (yielding 609 S-D measurement pairs). The medium consisted of 3564 voxels of size 2*mm*^{3} (1mm × 1mm × 2mm). The temporal points were chosen to be 200*ps* apart, corresponding to the typical minimum gate width in time-gated detection techniques [10, 13], and spread across a time range of 6*ns*. To begin with, consider performing tomography using Eq. 14 with all S-D pairs and a single time point. What is the location for this time point for an optimal reconstruction? To answer this question, the singular value spectra for *W _{n}* evaluated at various time points were calculated. The five spectra with the highest values are plotted in Fig 3(a), and the number of singular values,

Measurement geometry and arrangement of sources (*) and detectors (o) used for the simulations. The medium was assumed to be a diffusive slab of thickness 2cm. The targets used in the tomographic reconstructions are shown as gray shaded areas. (a) shows **...**

Next, SVD was performed on the weight matrix calculated for all possible pairs of time points, and the pair of time points with the highest number of singular values, *N*_{σ}, was determined. It turns out that one of the time points was again near the peak of the DFTR and the other was located near the rise portion of the DFTR. In the same way, the location and *N*_{σ} for multiple combinations of time points were determined. In Figure 3(b), *N*_{σ} is plotted as a function of the number of time points used. It seen that the proportional increase in *N*_{σ} diminishes rapidly after the first 3 or 4 temporal measurements. Also shown in the inset in Fig 3(b) are the first five optimal time points on a representative DFTR on the surface, which are located near the initial portion of the DFTR before the beginning of the fluorescence decay. It was determined that additional time points were located near the decay portion of the DFTR and added little to *N*_{σ}.

While the exact location of the optimal time points might vary slightly depending on the specific medium geometry, S-D arrangement and the location of the heterogeneity, it is generally clear from the results in Fig. 3 that the most useful time points of the DFTR for a direct TD reconstruction are located near the rise and peak portions. This result is consistent with Eq. 13, which shows that the weight function is asymptotically factorized into a spatial and temporal component so that multiple time points in the decay region are redundant for tomography. In other words, a brute force direct TD approach is *not* ideal for tomography with the long time decay data. (When the lifetimes are widely separated, the shorter lived components may be suppressed by reconstructing later delays [16].) Instead, the asymptotic approach based on a derived data type, viz., the inverse laplace transform (i.e., multi-exponential fit) is more appropriate. In the next section, we will perform tomographic reconstructions with simulated data using realistic noise levels to more quantitatively study the imaging performance of the direct TD and asymptotic approaches.

The results presented so far in the paper suggest that time domain fluorescence tomography can be comprehensively performed in three simple steps. (1) Obtain the *intrinsic* fluorescence lifetime(s) and the corresponding decay amplitude(s) from the asymptotic tail. (2) Reconstruct the individual yield distributions for each decay component using the decay amplitudes for all S-D pairs. (3) Reconstruct the yield distributions for each lifetime component using a few time points near the rise and the peak of the DFTR. These three steps reduce the computational complexity involved in a brute-force reconstruction of a large temporal data set, while retaining the most complete information possible from a TD experiment.

The question immediately arises as to the relative performance of the direct TD and asymptotic approaches. To address this, tomography was performed on simulated noisy data. The simulations employed the same slab geometry used in the SVD analysis above, with two fluorescent inclusions positioned symmetrically with respect to the medium geometry and S-D arrangement (Fig. 2) to remove any intrinsic bias due to the point spread function. The forward data was simulated with a shot noise model, which is characteristic of photon counting detection schemes, for three laterally placed inclusions (Fig. 2(a)) with center-to-center separations of 7mm, 5mm and 3mm. Also considered was the case where the inclusions had non-zero axial separation of 4mm, with zero lateral separation (Fig. 2(b)). The inclusions had distinct lifetimes of 1ns and 1.5ns. The regularization was carried out using a Moore-Penrose inversion algorithm, using the pre-calculated SVD matrices, *U,S,V* of the weight matrix *W _{n}*. Denoting the measurement vector by

$$X=VS{({S}^{2}+\alpha \lambda I)}^{-1}{U}^{T}y$$

(16)

where $\alpha =\text{max}\left\{\text{diag}\left({W}_{n}^{T}{W}_{n}\right)\right\}$ and the regularization parameter λ is typically between 10^{−5} and 10^{−3}. Three different reconstructions were performed, namely, CW, direct TD, and asymptotic. The CW reconstructions were performed using the time integrated TD data. The direct TD reconstructions used Eq. (14) with a set of 4 time points on the rising edge of the DFTR, following the SVD analysis results of Fig. 3. The asymptotic TD reconstruction was performed using the amplitudes obtained from a multi-exponential analysis of the decay portion in Eq. 15.

The sensitivity of the reconstructed images to measurement noise was quantified by simulating 100 data sets with noise for each S-D pair and time gate. The contrast-to-noise ratio (CNR), and the full-width-half maximum (FWHM) were then calculated as a function of λ . The CNR is estimated as the ratio of the peak image intensity in a region of interest surrounding the known location of the inclusion, to the mean standard deviation of the voxels in that region. The FWHM was estimated as the cube-root of the total volume of the voxels within half the peak intensity. In addition, for the lifetime based TD and asymptotic reconstructions, which provide separate yield reconstructions η_{1} and η_{2} for the 1ns and 1.5ns lifetimes, the cross-talk X was estimated to quantify the separability of the two inclusions based on lifetimes. If Ω_{1} denotes a chosen region-of-interest around the known location of the 1ns inclusion, then X^{1ns} = max[*η*_{2}(Ω_{1})]/max[*η*_{1}(Ω_{1})]. The yield cross talk for the 1.5ns component was similarly evaluated. The CNR vs FWHM plot is shown in Fig. 4 for the 1ns lifetime inclusion, for the case with 7mm lateral separation between the inclusions. It is clear that the TD reconstruction shows a dramatic improvement in the CNR and FWHM over the asymptotic reconstruction, and an improvement over the resolution of the CW case. The CNR improvement is evidently due to the better SNR of the peak portion of the DFTR compared to the asymptotic tail. The FWHM improvement of the TD over CW is due to the tomographic separation of the yield distributions for the lifetime components. Thus, for fixed CNR, the lifetime based TD reconstruction will have superior resolution compared to the asymptotic and CW reconstructions. However, the cross-talk, (which is the reconstructed amplitude of the 1.5ns inclusion at the location of the 1ns inclusion) is significantly higher for TD than the asymptotic case, and is attributable to the non-diagonal nature of the TD portion of the forward problem in Eq. (14). We note that the crosstalk for the asymptotic approach will depend on the separability of the lifetimes from the multi-exponential fits of raw experimental data, an aspect that will be explored in future work.

Plot of the contrast-to-noise ratio (CNR) vs full-width-half-maximum (FWHM) for reconstructions using CW (solid black), optimal direct TD using 4 most significant time points near the rise and peak of the DFTR (red), and the asymptotic approach (blue **...**

The effect of the crosstalk can be more clearly seen in the reconstructed tomographic images shown in Fig. 5, where the X-Z slices of the 3-D reconstructions for all three data sets and separations are displayed. Also shown are the plots of the yield at a fixed depth (Z) where the yield is maximum. It should be noted that the regularization λ was not identical for all the reconstructions but was rather determined by the condition that *log*_{10}(CNR) was near 1. This is necessary to properly account for the difference in the noise characteristics of the different methods. (For example, CW has the best SNR, and should thus be the least regularized.) To visualize crosstalk easily, the yield images for the 1ns and 1.5ns components for the TD and asymptotic approaches were assigned red and green colors in an RGB colormap of a single image. The degree of crosstalk is thus revealed as a mixture of these two colors (e.g., yellow implies 100% crosstalk). Thus, CW reconstructions have no lifetime information so that they are shaded in yellow. It is clear from Fig. 5 that the TD reconstruction has superior resolution but significantly more cross-talk than the asymptotic reconstructions, as can also be seen in the intensity plots in the bottom panel. For small target separations, the cross talk of the TD method proves detrimental to its accuracy, whereas the asymptotic case recovers the localizations accurately even for 3mm separation. Thus it can be said that the direct TD approach using optimal time points provides more precise (better-resolved) reconstructions, and is useful when the targets are well separated, whereas the asymptotic reconstructions are more accurate but less precise (less-resolved). In Figure 6, the reconstructions are shown with the targets located axially, i.e., along the S-D axis. The advantage of the lifetime based asymptotic reconstruction is even more evident in this case, as the localizations of the two lifetimes are not reproduced either for the CW or the direct TD reconstructions.

X-Z slices of the 3-D Reconstructions of two targets with separations of 7mm, 5mm and 3mm, and with lifetimes of 1ns and 1.5ns, located transverse to the S-D axis, using CW (a-c), direct TD (d-f) and asymptotic (g-i) data sets. The measurement geometry **...**

Although characterized by cross-talk, the presence of lifetime contrast enhances the images reconstructed using direct TD data. To delineate this effect clearly, Fig. 7 shows a comparison of the reconstructions of the two laterally placed targets using the direct TD approach, with and without lifetime contrast between the targets. The effect of lifetime contrast on the direct TD reconstruction is clearly seen as a significant reduction in the point-spread-function of the reconstructed yield distributions for the two lifetimes. Of course, this effect will depend on the difference between the fluorophore lifetimes and the diffuse propagation time *τ _{D}*, which is near 0.4ns for the present simulation. (corresponding to μ

We have presented a theoretical formalism for TD fluorescence tomography with turbid media that allows a lifetime-based reconstruction of yield distributions using the entire TD data. Besides providing a comprehensive understanding of TD fluorescence signals from diffuse media, a key advance of this work from the previously presented asymptotic lifetime based tomography is an algorithm for lifetime based tomographic separation using the peak and rise portions of the temporal diffuse fluorescence response. The formalism is generally valid provided the intrinsic fluorescence lifetimes are revealed in the long time decay, a condition well satisfied [11–13] given the typical optical response times of diffuse tissue and fluorophore lifetimes used in molecular imaging. This is important since the measured average lifetime on the surface can by itself provide useful diagnostic information, without the need for tomography. This shows the potential for lifetime sensing in diagnostic imaging and extends the application of this work to a wide range of biomedical imaging problems.

The results presented here can be viewed as an inevitable consequence of the long lifetime condition: The longer lived fluorescence decay effectively convolves over the intrinsic diffuse material response, resulting in a decay tail that is separated in space and time. This implies that for tomography with the long time data, a direct use of multiple time points in the decay portion is redundant. The optimal approach is to perform tomography using the amplitudes recovered from multi-exponential fits. This result was shown to be consistent with a SVD analysis of the time-dependent Born weight functions, which showed that the optimal time points to use in a direct TD reconstruction are located near the peak and rise portions of the DFTR.

Tomographic reconstructions with simulated noisy data also revealed the relative merits of optimized direct TD and the asymptotic approaches. It was found that the direct TD and asymptotic approaches yield complementary information: The asymptotic approach provides superior localization ability due to minimal cross-talk between the images for multiple lifetimes, while the optimal direct TD reconstructions yield better resolution due to superior SNR near the peak of the DFTR. Thus, when no lifetime contrast is present in the medium, the direct TD analysis should be the method of choice. For targets located along the S-D axis, it was shown that the asymptotic analysis is superior to both direct TD and CW in its ability to accurately localize the targets, provided they have different lifetimes. Axially located targets could occur for example in small animal brain imaging, where transillumination may be the preferred geometry when depth resolution along the various brain regions is desired.

The simulations presented here considered two distinct lifetime inclusions placed both lateral and axial to the measurement axis. Although simplistic, this example has highlighted key aspects of TD fluorescence tomography with lifetime contrast that can potentially be extended to more complicated spatial distributions of lifetimes. This work is also based on the assumption that lifetimes present in the medium are few and discrete, or at least can be described as sharp distributions centered around a mean lifetime. In the more general case when lifetimes are broadly distributed, a numerical inverse Laplace transform can be used to recover amplitude distributions [26]. The simulation analysis presented here was not meant to optimize for any particular TD detection technique (e.g., wide-field time gated, time correlated detection schemes), but was rather an attempt to explore the information content in a TD signal and to provide a recipe for an optimal approach for TD fluorescence tomography with turbid media. The purpose of the numerical simulations was also not to make a statement about the absolute resolution achievable by TD methods. This quantity can be optimized using better S-D arrangement and adjusting actual experimental conditions. Indeed, sub-mm resolution has recently been demonstrated using CW measurements [3]. Such optimization will enhance the resolution of all three approaches viz., CW, direct TD and asymptotic, so that the main results obtained here will not be affected.

In future work, we will attempt to extend the formalism developed here to a hybrid model that incorporates both the direct TD and the asymptotic approaches into a single inverse problem in a self-consistent fashion. This hybrid TD-asymptotic approach is expected to provide optimal localization and resolution. It should be reiterated that although a diffusive slab model was assumed for the simulations, the formalism developed here can readily incorporate Green’s functions calculated as solutions of either the diffusion or transport equations, as appropriate, for heterogenous media with complex boundaries. We are currently engaged in applying the formalism developed here to imaging complex shaped mouse phantoms and mouse models of disease.

Optical molecular imaging can immensely benefit from the use of biochemical reporter probes that are not merely disease specific, but also provide specific molecular signatures such as spectral and lifetime shifts that help isolate the disease from background tissue [2]. The unique advantages of time domain technology as explored in this work strongly motivates the development of fluorescent contrast agents that exhibit target specific lifetime shift upon binding. We also hope that this work will provide useful guidelines for biological imaging using time domain fluorescence tomography.

This work was supported by National Institutes of Health Grants EB000768, AG026240 and P41-RR14075.

1. Hintersteiner In-vivo detection of amyloid-beta deposits by near-infrared imaging using an oxazine-derivative probe. Nat. Biotechnol. 2005;23:577. [PubMed]

2. Ntziachristos V, Ripoll J, Wang LV, Weissleder R. Looking and listening to light. Nat. Biotechnol. 2005;23:314. [PubMed]

3. Graves EE, Ripoll J, Weissleder R, Ntziachristos V. A submillimeter resolution for small animal imaging. Med. Phys. 2003;30:901. [PubMed]

4. Oleary MA, Boas DA, Li XD, Chance B, Yodh AG. Fluorescence lifetime imaging in turbid media. Opt. Lett. 1996;21:158. [PubMed]

5. Godavarty A, Sevick-Muraca EM, Eppstein MJ. Three-dimensional fluorescence lifetime tomography. Med. Phys. 2003;48:1701–1720. [PubMed]

6. Das BB, Liu F, Alfano RR. Time-resolved fluorescence and photon migration studies in biomedical and model random media. Rep. Prog. Phys. 1997;60:227.

7. Chen K, Perelman LT, Zhang QG, Dasari RR, Feld MS. Optical computed tomography in a turbid medium using early arriving photons. J. Biomed. Opt. 2000;5:144. [PubMed]

8. Wu J, Perelman L, Dasari RR, Feld MS. Fluorescence tomographic imaging in turbid media using early-arriving photons and Laplace transforms. Proc. Natl. Acad. Sci. USA. 1997;94:8783. [PubMed]

9. Hall D, Ma G, Lesage F, Wang Y. Simple time-domain optical method for estimating the depth and concentration of a fluorescent inclusion in a turbid medium. Opt. Lett. 2004;29:2258. [PubMed]

10. Turner GM, Zacharakis G, Sourbet A, Ripoll J, Ntziachristos V. Complete-angle projection diffuse optical tomography by use of early photons. Opt. Lett. 2005;30:409. [PubMed]

11. Ma G, Mincu N, Lesage F, Gallant P, McIntosh L. System irf impact on fluorescence lifetime fitting in turbid medium. Proc. SPIE. 2005;5699:263.

12. Bloch S, Lesage F, Mackintosh L, Gandjbakche A, Liang K, Achilefu S. Whole-body fluorescence lifetime imaging of a tumor-targeted near-infrared molecular probe in mice. J. Biomed. Opt. 2005;10:054003–1. [PubMed]

13. Kumar ATN, Skoch J, Bacskai BJ, Boas DA, Dunn AK. Fluorescence lifetime-based tomography for turbid media. Opt. Lett. 2005;30:3347–3349. [PubMed]

14. Hattery D, Chernomordik V, Loew M, Gannot I, Gandjbakhche A., J. Opt. Soc. Am. A, Analytical solutions for time-resolved fluorescence lifetime imaging in a turbid medium such as tissue J. Opt. Soc. Am 2001. 181523–1530.1530 [PubMed]

15. Lam X, Lesage F, Intes X. Time domain fluorescent diffuse optical tomography:analytical expressions. Opt. Express. 2005;13:2263–2275. [PubMed]

16. Apreleva SV, Wilson DF, Vinogradov SA. Feasibility of diffuse optical imaging with long-lived luminescent probes. Opt. Lett. 2006;31:1082–1084. [PMC free article] [PubMed]

17. Sevin PR. The renaissance of fluorescence resonance energy transfer. Nat. Struct. Biol. 2000;7:730–734. [PubMed]

18. Bastiaens PIH, Squire A. Fluorescence lifetime imaging microscopy: spatial resolution biochemical processes in the cell. Trends Cell. Biol. 1999;9:48–52. [PubMed]

19. Bacskai BJ, Skoch J, Hickey GA, Allen R, Hyman BT. Fluorescence resonance energy transfer determinations using multiphoton fluorescence lifetime imaging microscopy to characterize amyloid-beta plaques. J. Biomed. Opt. 2003;8:368–375. [PubMed]

20. Equation (4) is a generalized form of a semi-analytical expression for an infinite medium given in Reference 9, as can be checked by setting *t*’ → *t* — *t*’ and using the analytical expression for the Greens functions with ${\mu}_{a,s}^{x}={\mu}_{a,s}^{m}$

21. Arridge SR. Optical tomography in medical imaging. Inverse Problems. 1999;15:R41.

22. Chandrasekhar S. Radiative Transfer. Dover; New York: 1960.

23. Gao F, Zhao H, Tanikawa Y, Yamada Y. A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography. Opt. Express. 2006;14:7109–7124. [PubMed]

24. Patterson MS, Chance B, Wilson BC. Time-resolved reflectance and transmittance for the non-invasive measurement of tissue optical properties. Appl. Opt. 1989;28:2331–2336. [PubMed]

25. Haselgrove JC, Schotland JC, Leigh JS. Long-time behavior of photon diffusion in an absorbing medium: application to time-resolved spectroscopy. Appl. Opt. 1992;31:2678–2683. [PubMed]

26. Kumar ATN, Zhu L, Christian JF, Demidov AA, Champion PM. On the Rate Distribution Analysis of Kinetic Data Using the Maximum Entropy Method: Applications to Myoglobin Relaxation on the Nanosecond and Femtosecond Timescales. J. Phys. Chem. B. 2001;105:7847–7856.

27. Milstein A, Stott JJ, Oh S, Boas DA, Millane RP, Bouman CA, Webb KJ. Fluorescence optical diffusion tomography using multiple-frequency data. J. Opt. Soc. Am A. 2004;21:1035–1049. [PubMed]

28. Kumar ATN, Skoch J, Hammond FL, Dunn AK, Boas DA, Bacskai BJ. Time resolved fluorescence imaging in diffuse media. Proc. of SPIE. 2005;6009:60090Y.

29. Graves EE, Culver JP, Ripoll J, Wessleder R, Ntziachristos V. Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography. J. Opt. Soc. Am. A. 2004;21:231–241. [PubMed]

30. Culver JP, Ntziachristos V, Holboke MJ, Yodh AG. Optimization of optode arrangements for diffuse optical tomography: A singular-value analysis. Opt. Lett. 2001;26:701–703. [PubMed]

31. Viswanath K, Mycek M. Do fluorescence decays remitted from tissues accurately reflect intrinsic fluorophore lifetimes? Opt. Lett. 2004;29:1512–1514. [PubMed]

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