A Time Domain Fluorescence Tomography System for Small Animal Imaging

doi:10.1109/TMI.2008.918341

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 August 11.

Published in final edited form as:

IEEE Trans Med Imaging. 2008 August; 27(8): 1152–1163.

doi: 10.1109/TMI.2008.918341

PMCID: PMC2920137

NIHMSID: NIHMS221199

A Time Domain Fluorescence Tomography System for Small Animal Imaging

Anand T. N. Kumar, Scott B. Raymond, Andrew K. Dunn, Brian J. Bacskai, and David A. Boas

Anand T. N. Kumar, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Charlestown, MA 02129 USA;

Contributor Information.

Corresponding author.

Anand T. N. Kumar: ankumar/at/nmr.mgh.harvard.edu

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

Abstract

We describe the application of a time domain diffuse fluorescence tomography system for whole body small animal imaging. The key features of the system are the use of point excitation in free space using ultrashort laser pulses and noncontact detection using a gated, intensified charge-coupled device (CCD) camera. Mouse shaped epoxy phantoms, with embedded fluorescent inclusions, were used to verify the performance of a recently developed asymptotic lifetime-based tomography algorithm. The asymptotic algorithm is based on a multiexponential analysis of the decay portion of the data. The multiexponential model is shown to enable the use of a global analysis approach for a robust recovery of the lifetime components present within the imaging medium. The surface boundaries of the imaging volume were acquired using a photogrammetric camera integrated with the imaging system, and implemented in a Monte-Carlo model of photon propagation in tissue. The tomography results show that the asymptotic approach is able to separate axially located fluorescent inclusions centered at depths of 4 and 10 mm from the surface of the mouse phantom. The fluorescent inclusions had distinct lifetimes of 0.5 and 0.95 ns. The inclusions were nearly overlapping along the measurement axis and shown to be not resolvable using continuous wave (CW) methods. These results suggest the practical feasibility and advantages of a time domain approach for whole body small animal fluorescence molecular imaging, particularly with the use of lifetime as a contrast mechanism.

Index Terms: Fluorescence diffuse optical tomography, life-time-based sensing, molecular imaging, time-resolved imaging

I. Introduction

Fluorescence diffuse optical tomography [1]–[35] based on extrinsic contrast agents offers several advantages over intrinsic absorption contrast diffuse optical tomography. Most important of these is the use of the red-shifted fluorescence emission to isolate signal arising from the fluorophores inside the tissue, while rejecting the incident light at the excitation wavelength. Recent advances in the engineering of disease-targeted fluorescent markers have significantly contributed to the success of this technique by enabling specific in vivo fluorescence labeling of molecular events and tumors [36]–[56].

Fluorescence tomographic imaging can be broadly classified based on the excitation mechanism employed: continuous-wave (CW) methods use steady state excitation sources [1]–[11], frequency-domain (FD) methods use frequency modulated (~ MHz) sources [12]–[26], and time-domain (TD) methods [27]–[35] use short (~ fs–ps) laser pulses in conjunction with time resolved detection. For detection, the use of wide-field, noncontact imaging using charge-coupled devices (CCD), as opposed to fiber-based approaches, offers several advantages as exploited previously [8]–[10], [18], [19]. Noncontact detection offers the flexibility of optimizing detector choice and positioning on the imaging surface, while also avoiding fiber-tissue contact effects. Noncontact imaging can also do away with the need for matching fluids to simplify the imaging geometry, provided the surface boundaries of the imaging subject can be determined.

The two primary quantities of interest in fluorescence tomography are the fluorescence yield and the fluorescence lifetime [57]. While CW, FD, and TD approaches can recover the fluorescence yield, only TD and FD systems are capable of measuring fluorescence lifetime. In addition, TD measurements offer the most comprehensive information among the three techniques, since a single laser pulse implicitly contains all modulation frequencies, including the CW (zero frequency) component. The use of high-density CCD array detectors and time gated multichannel image intensifiers allows for multiple temporal measurements to be rapidly acquired over the surface of a biological imaging subject.

The combination of time-resolved measurements with CCD based detection thus implies a potentially large amount of information. Each source-detector (S–D) measurement in a TD experiment also contains the temporal response (on the subnanosecond time scales) of the turbid imaging medium. An experimental implementation of TD fluorescence tomography therefore has to address the question of how this vast information can be optimally utilized. We recently showed that, provided the fluorescence lifetime is longer than the intrinsic absorption time scales present within the medium, the diffuse fluorescence temporal response (DFTR) separates into an early diffuse portion and an asymptotic fluorescent decay portion [32], [33]. Also, the spatio–temporal evolution of the DFTR in the asymptotic regime is factorized into a product of CW spatial sensitivity functions and pure exponential decays. This observation implies that fluorescence lifetimes can be measured directly from surface measurements using multiexponential fits to the asymptotic decay. It also simplifies the analysis of a large portion of tomographic TD data and allows for an optimal way of localizing multiple lifetime components [32]. We note that FD techniques [12]–[26] and derived data types from TD measurements such as the Fourier [58], Laplace [28], or Mellin transforms (i.e., moments) [59] do not achieve the separation of multiple lifetimes directly from the raw data. Also, each derived data set (e.g., Fourier component) of the DFTR is related to an integral over the entire DFTR, so that the asymptotic decay portion is used in a suboptimal fashion.

While fluorescence lifetime tomography has been discussed previously from FD measurements [12]–[26], this paper demonstrates the application of a time-resolved imaging system to recover lifetime-based tomographic images. The boundary surface of the imaging subject is acquired using photogrammetric methods and employed in Monte-Carlo models of photon propagation. We show that the use of noncontact imaging with TD data offers certain advantages in addition to those explored in previous work [8]–[10], [18], [19]. For instance, the system impulse response function and the time origin can be measured by imaging the sources directly. This procedure avoids the necessity of determining these quantities from complicated fitting procedures. We consider the application of the system for small animal imaging using a mouse shaped phantom with embedded fluorescence inclusions and with tissue mimicking optical properties. The detailed experimental and numerical modeling steps necessary for the application of asymptotic lifetime-based-tomography are presented. Finally, the accuracy of the phantom results are verified by comparison with computed tomography (CT) images of the mouse phantom.

II. Theory: Review of the TD Fluorescence Forward Problem

A theory of TD fluorescence in turbid media was recently presented [33]. For completeness, we review the key results in this section. Consider a turbid imaging medium embedded with fluorophores, assumed to have discrete lifetime components τ_n = 1/Γ_n with associated yield distributions η_n(r). For tomography, optical sources, and detectors are arranged on the surface of the imaging medium. The fluorescence intensity measured at a wavelength λ_m, at detector position r_d and time t, for an impulsive excitation at wavelength λ_x, position r_s and t = 0, can be expressed as a multiexponential sum

(1)

It should be noted that the decay amplitudes A_n depend on time, in addition to the source and detector locations. Within the Born approximation, which assumes a single fluorescence emission event, the A_n’s define a linear inverse problem for the yield distributions of each lifetime component¹

(2)

where

(3)

The Green’s functions equation M4

used in (3) are simply the time-dependent background Green’s functions for light propagation (described by the transport or diffusion quations [67]) at the excitation and emission wavelengths, G^x and G^m, but evaluated with reduced absorptions, equation M5

and

, where v denoting the velocity of light in the medium.

The set of expressions (1)–(3) offers a compact and rigorous description of fluorescent light transport within turbid media. These expressions directly incorporate the lifetimes obtained from the decay of TD data into the reconstruction for the yield distributions. It should be noted that the Born approximation is only made for the excitation-fluorescence emission process: multiple interactions of the light with the background medium along the excitation and emission paths are implicitly accounted for through the dependence of the Green’s functions equation M7

and

on the respective total background absorptions equation M9

and

. The background absorptions can in general include the unknown fluorophore absorption in a nonlinear treatment. Also, equation M11

and

can be most generally obtained as solutions to the diffusion or transport equations at the excitation and emission wavelengths for arbitrary heterogenous media bounded by complex surfaces.

The theory developed in (1)–(3) can be applied to the entire TD data. This paper will however be concerned with the asymptotic or long time portion of the DFTR. At times longer than the evolution of the diffusive response of the background medium τ_d (which is smaller than the intrinsic absorption timescale τ_a = (vμ_a)⁻¹ [33], [60], [61]), the TD expression in (1) smoothly goes over to an asymptotic limit where the amplitudes A_n become time independent [33]. In this limit, the fluorescence decay factorizes into spatial and temporal components

(4)

The amplitudes, a_n, are the long time values of A_n, obtained by extending the upper limit of the time-integral on the right-hand side of (2) to ∞ with the understanding that the equation M14

has completed its evolution and is constant in the asymptotic region [33]

(5)

where the CW weight functions equation M16

are given by

(6)

Realizing from (6) that equation M18

is just the zero-frequency component of the Fourier transform of equation M19

, which according to (3) is a convolution, it is readily shown that

(7)

where

and

are the CW counterparts of equation M23

and

. The weight functions equation M25

are therefore just the CW, Born sensitivity functions for an absorption perturbation, but evaluated with the background absorption reduced by Γ_n/v. The key point to note from the above expressions is that the amplitudes a_n that form the measurement set for the linear inverse problem in (5), are identified directly with the coefficients of multiexponential fits recovered from the surface decay profiles. We thus have a scenario where the forward problem has been completely separated into multiple lifetime components, consistent with the asymptotic portion of the raw TD data itself. This allows for the clean tomographic separation of multiple lifetime components present within the imaging medium. In this work we will focus on the practical implementation of the asymptotic TD theory as given in (5) for complex media, using an experimental tomographic imaging setup designed for small animal imaging applications.

It is useful at this stage to write the CW component of the net fluorescence signal, which is obtained by simply integrating (1) over all times. It is readily found using (2) and (3) that

(8)

where ^B are the standard CW Born-sensitivity functions for an absorption perturbation [67]. Comparison of the time-domain asymptotic problem given by (5) with the CW forward problem given by (8), clearly illustrates the advantage of the TD case. Although both the asymptotic TD and the CW measurements are related to a CW weight matrix, the TD measurement effectively separates into multiple CW forward problems for each lifetime. A CW measurement, however, is related to the sum of the yield distributions for all the lifetime components present in the imaging medium.

III. Materials and Methods

A. TD Apparatus

All measurements were carried out using a small animal molecular imaging system, schematically depicted in Fig. 1. The primary excitation source was a Ti:Sapphire laser (Mai Tai, Newport-Spectra Physics, Mountain View, CA ~ 150-fs pulse width, 80-MHz repetition rate, 750–850-nm tuning range). To broaden the range of excitation wavelengths, a super-continuum source setup was also constructed using poly-crystalline fiber (PCF) technology. The PCF (Femtowhite 800, Crystal Fibre, Denmark) had a 1.8-μm core diameter and a zero dispersion wavelength of 750 nm. For optimal performance, the input to the PCF was near 800 nm, which is in the “anomalous” dispersion regime [63]. The inset in Fig. 1 shows a typical output power spectrum from the PCF measured for ≈ 500-mW input at 800 nm. The output of the PCF had a > 500-mm spectral range, with pulse widths of a few tenths of picoseconds and average power of a few mW/10 nm. For applications requiring wavelengths in the range 750–850 nm, as in the experimental results reported in this paper, the Ti:Sapph was used directly as the excitation source.

Fig. 1

Schematic of the free space TD fluorescence tomography system. The source consists of either a Ti:Sapphire laser (750–850 nm) or a photonic crystal fiber (PCF) driven by the Ti:Sapph output. The light is launched into a 200 μm step-index (more ...)

The light output from the Ti:Sapph was attenuated to 5–20 mW depending on the sample thickness and absorption, and launched into a 200-μm core SMA connectorized fiber using a fiber collimation package (FC). (All opto-mechanical components were from Thorlabs Inc., except for components in the PCF source setup, which were from Newport Corporation). The light at the output end of the fiber was focused, using another collimation package, to a fine spot on the imaging surface, resulting in an approximately point excitation source. The collimation package system was mounted on a computer controlled, micrometer precision XY translation stage system (Velmex Inc., Bloomfeld, NY), which allowed for scanning the source over the imaging surface. The source fiber with the collimation system was placed below the phantom plate for transmission measurements. Reflectance measurements could also be performed, if needed, by mounting the fiber above the phantom plate (not shown in Fig. 1).

The detection was performed in free space using a CCD camera (Picostar HR-12 CAM 2, LaVision GmbH, Goettingen, Germany) with a quantum efficiency of 65%. A time gated image intensifier (Picostar HR-12, LaVision GmbH, Goettingen, Germany), with a 200-ps minimum gate width, was used to provide nanosecond time resolution. A fast HRI delay unit (1-ms switching, 25-ps minimum steps) was used to delay the trigger output (80 MHz) from the laser to the intensifier unit, across one duty cycle of 12.5 ns. A 60-mm focal length camera lens (AF Nikkor, f2.8, Nikon) mounted on the CCD camera allowed a well-focused image of the entire surface of the mouse head to be formed on the CCD image plane. The maximum resolution of the CCD camera was 1376 × 1040 pixels. For better noise statistics, the camera images were 4 × 4 binned (averaged) at the hardware level to a reduced size of 344 × 260 pixels. Filters could be attached optionally to the camera lens to collect multispectral, time-resolved data. Since the scanning speed of the stepper motor was much slower than the delay switching time, the full temporal signal was collected for one source position at a time. All control software was written within the software package (DAVIS 6.1) accompanying the LaVision camera.

B. Measuring the Impulse Response, Time Origin, and Detector Positions

In order to match the model prediction with the measured surface tomographic data, it is important to obtain the system impulse response function (IRF) as well as an accurate estimate of the initial time t₀. The time t₀ refers to the time when the excitation pulse is incident on the surface of the imaging medium. A correct estimate of t₀ ensures that the relative amplitudes of the multiple lifetime components are correct, and is crucial for estimating the optical properties using TD data, as well as for fluorescence reconstructions using the early time gates [33]. An advantage of noncontact detection is that t₀ can be estimated from the IRF, which can in-turn be measured directly from the source illumination. This is possible since the position of the CCD camera (i.e., the detectors) is unchanged before and after the sample is placed on the imaging plate, unlike in contact geometries with fiber based detection.

The excitation light (at 755 nm) was focused onto a piece of paper placed on the imaging plate (without the imaging subject). The full temporal response, which is the IRF, was then measured using a 10-nm band pass filter at 790 nm (to minimize damage to the intensifier unit), for each position of the source. Fig. 2(a) shows the integrated (CW) IRFs for a grid of 21 source positions spaced 2 mm along the X and Y translation directions. The source grid image in Fig. 2(a) also serves to determine the coordinates of the sources for forward modeling. These were obtained as the points of maximum intensity, after the image for each source position was median filtered. Fig. 2(b) shows the corresponding IRF’s for all source positions. The width of the IRF was nearly equal to the gate width setting of 300 ps. The asymmetry in the response could be attributed to a difference in the charging-discharging speeds of the intensifier, and could also be caused by a weak autofluorescence from the paper. Since the variation between the IRF was minimal across the sources (< 10 ps), the initial time t₀ was assumed to be identical for all source positions and estimated from the mean of the impulse response as the time at 1% of the peak. This value can be assumed to be the initial excitation time at the surface of the imaging medium, provided the boundary of the imaging volume is approximately flat. For the experiments reported here, the range of source positions was limited to a relatively flat region in the bottom surface of the mouse phantom.

Fig. 2

(a) Direct CCD image of the source illumination, obtained by focusing the fiber output on a thin sheet of paper placed on the imaging plate (without the imaging subject present). The image shown is the integrated TD signal for 21 source positions, used (more ...)

Note that the time t₀ measured from the IRF includes the time for free space propagation of light from the source to the camera. However, the propagation time for the fluorescence emitted from the imaging surface to the detector will be slightly shorter than t_0, due to the finite thickness of the mouse. Before analyzing the fluorescence data, t₀ was therefore offset by ≈ 50 ps, corresponding to a mouse thickness of 2 cm.

The measured IRF can be directly forward convolved into the model before tomographic inversion. This procedure is superior to a deconvolution of the IRF from the raw fluorescence data, which is a highly ill-posed problem. From Fig. 2(b), it is clear that the IRF can reasonably well be approximated as a square function. For the asymptotic reconstructions, which is the main focus of this paper, the effect of the IRF on the forward model can then be readily calculated. Using (4), the signal for a square IRF of width T is expressed as

(9)

Thus, the effect of the system IRF on the measured decay amplitudes is a scaling factor that depends on the lifetime and the width of the IRF.

C. Mouse Phantom Preparation

Mouse shaped phantoms were prepared using the negative mold of a sacrificed mouse. The mold for the phantom consisted of a two-part epoxy resin and hardener mixture (Douglas and Sturgess Inc., San Francisco, CA). A combination of TiO₂ (paint) and ink was added to the resin, to achieve approximate background optical properties of equation M28

and μ_a = 0.1 cm⁻¹ [68]. The optical properties were also independently verified in the reconstruction stage from the excitation and emission data, as described in Section III-F. Two inclusions (cavities) were created inside the phantom by positioning low-melting-point agar beads (5 mm diameter) that were originally cast in a spherical negative mold, using pairs of syringe needles. The epoxy mixture was then poured into the negative mold and allowed to cure for 12 h. Once cured, the negative mold was removed along with the needles. Fresh needles were then attached along with polypropylene tubing, enabling dynamic injection of fluorophores into the inclusions. Fig. 3 shows a photograph of the mouse phantom along with the injection tubes.

Fig. 3

Photograph of epoxy mouse phantom used in the experiments. The fluorophore injection needles (two for each inclusion) are also indicated. (A → top inclusion; B → bottom inclusion.)

To determine the location of the inclusions, CT images of the mouse phantom were performed using a ultra-high resolution, flat panel CT scanner (Siemens Medical Solutions, Forchheim, Germany) at MGH. The inclusions were filled with an aqueous solution of Iodine. Fig. 8 displays the CT reconstructions (as grayscale images) along two sagittal planes that roughly contained the centroid of the inclusions. Note that the inclusions (visible as white spots on the gray mouse background) were not spherical and were smaller than the original size of the spherical beads. Moreover, the inclusions were much less clearly defined and had a larger spatial extent along the normal to the sagittal plane, due to the presence of the tubing and injection needles. The centroids of the inclusions were estimated from the CT image to be separated by 6 mm along the vertical dimension, and located at depths of 4 and 10 mm from the top surface of the mouse head. The separation of the centroids along the length of the mouse was near 1.6 mm. The separation of the inclusions along the width of the mouse (normal to the sagittal) was not clearly defined in the CT reconstruction, due to the presence of the injection tubing and needles.

Fig. 8

Normalized depth profiles for the asymptotic and CW reconstructions shown in Fig. 6. The profiles for each case show the depth dependence along Z, at the X–Y location of the peak of the corresponding yield reconstruction. The depth profiles recovered (more ...)

D. Three-Dimensional Camera for Surface Coregistration

For accurate forward modeling of the light propagation, the boundary information of the phantom or animal was obtained using a photogrammetric camera system (3-D Facecam 100, Genex Technologies Inc., Kensington, MD). The 3-D camera captures a surface image of the mouse as triangulated vertex data with a resolution of 100 μm. This high density data is first reduced to a set of unique surface mesh points, which were then interpolated using the Matlab (The Mathworks, Natick, MA) function “ndgrid,” to provide a surface on a uniform mesh. Before being employed in modeling, the 3-D surface image needs to be coregistered with the planar image of the CCD camera with which the fluorescence data are acquired. The coregistration procedure used here is largely similar to previous work [9], and will not be discussed in detail. Briefly, the transformation was effected as a 2-D affine transformation between the CCD image and the 3-D surface image, using four fiducial points in the form of metal screws placed on the subject plate at points within the filed-of-view of the 3-D camera and the CCD camera. The mouse phantom was placed on the imaging plate and its boundary secured using additional screws for repeatability in positioning. The 3-D image was then acquired with the optional mirror in place (Fig. 1). A single 3-D image thus acquired along with the fiducials can be used as the reference image for any optical measurement using the phantom, provided that the fiducial points are visible on the CCD image and the phantom is located at the same relative position with respect to the fiducials.

The coregistration results for the phantom reported in this work are shown in Fig. 4. Fig. 4(a) shows the CCD image with the mouse phantom, the fiducials, and a grid of 93 detectors (used for the tomographic measurements reported here) arranged on the camera image. The detectors had a 2 mm sepa along the length of the mouse head and 1 mm separation along the width. Each detector was assigned a rectangular grid of 9 × 5 pixels (2 mm × 1 mm) and the signal was calculated as the mean within this area. Fig. 4(b) shows the planar image of the 3-D surface used for the coregistration, along with source and detector coordinates mapped on the surface image using the affine transformation.

Fig. 4

Coregistration of CCD camera image with 3-D camera surface image. (a) Room light illuminated image of the mouse phantom placed on the imaging plate. The + symbols indicate the fiducial points used for co-registration. The sources (x) determined from the (more ...)

E. Multiexponential Fits Using Global Lifetime Analysis

An important step in the asymptotic approach is the extraction of the decay amplitudes for all the lifetimes present in the medium. The recovery of lifetimes and decay amplitudes from multiexponential fits to decay data is a nonlinear problem, which is further complicated by the presence of noise. However, the asymptotic separation of the spatial and temporal components of the DFTR, as expressed in (4), allows for a significant simplification of the fitting process. Since according to the asymptotic model in (4), the set of discretized lifetime components τ_n are independent of the measurements locations (i.e., the S–D coordinates), the multiexponential analysis can be performed in two stages. First, the total measured DFTR is calculated as a sum over all S–D pairs, to obtain a high signal-to-noise ratio (SNR) temporal data set. This “global” signal is composed of all the lifetime components present in the system, and allows a more robust determination of the lifetimes through a nonlinear analysis. In the second step, the lifetimes determined from the surface-integrated DFTR are used in a linear fit of the DFTR for each individual S–D pair. Besides improving the robustness of the fitting procedure for the lifetimes and decay amplitudes, the global analysis is computationally much less cumbersome than performing a nonlinear fit for every S–D measurement. The nonlinear fits for the global DFTR were numerically implemented using the Matlab function, “fminsearch,” which employs an unconstrained nonlinear optimization using the Nelder–Mead simplex approach.

F. Computation of the Forward Problem

The computation of the forward problem was carried out using a transport based approach, via the Monte-Carlo (MC) method [70], [72]. The first step is the determination of the background optical properties at the excitation and emission wavelengths, which in general includes the fluorophore absorption. Typically, these are determined from separate measurements of the incident fields at the excitation and emission wavelengths U^x,x and U^m,m for all the S–D pairs and time gates. For the reconstructions reported here, it was sufficient to estimate the optical properties using the corresponding CW components Ū^x,x and Ū^m,m of the TD data as follows. It is known that due to nonuniqueness [69], absorption, and scattering properties cannot be uniquely determined from CW data alone. For the mouse phantom used here, we assume a fixed value for the reduced scattering at 10 cm⁻¹, known to within 5% from the phantom preparation, and given the high robustness of fluorescence yield reconstructions to errors in the estimation of background scattering [16]. It was also assumed that equation M29

given the small separation (40 nm) between the wavelengths λ_x and λ_m.

With a fixed scattering coefficient, the MC approach allows dramatic flexibility for estimating the background absorptions equation M30

and

. Each MC simulation stores the photon path histories between each S–D pair. Since the photon trajectory in a transport medium is independent of absorption, the photon weights for each S–D pair can be quickly re-estimated for any given absorption through a Beer-law type relationship [71]. If the MC generated CW Greens’s functions for a pair of optical absorptions equation M32

and

are

and

, the incident signal at λ_x and λ_m can be written as

(10)

(11)

where Θ(r_s, r_d) is an S–D dependent coupling strength and An external file that holds a picture, illustration, etc.
Object name is nihms221199ig1.jpg

An external file that holds a picture, illustration, etc.
Object name is nihms221199ig2.jpg

are scaling factors that account for other parameters including the filter attenuation coefficient, intensifier gain, camera aperture and laser power. Using the experimentally measured excitation and emission data Ū^x,x and Ū^m,m (10) and (11) can be used to fit simultaneously for both the S–D coefficients and the absorption at the emission and excitation wavelengths [73]. Note that the background medium need not be homogenous, since the MC simulation stores the path histories for a large number of finite subregions within the simulation volume. Thus, the distribution of the optical absorption can be rapidly fit for any desired segmentation of the imaging medium. This is particularly advantageous in the common scenario of imperfect uptake or nonspecific binding of the disease markers [6], [23], where it is necessary to incorporate the highly absorbing fluorophore distribution present in the background medium into the forward model. Equations (10) and (11) thus offer a powerful and efficient nonlinear approach for a rapid determination of the background absorption (including fluorophore) and S–D coefficients, which are required to calculate the fluorescence yield distribution.

For the MC computation, the 3-D surface was digitized into a volume image by straight line descent along the Z-direction, for each X–Y coordinate. The mouse phantom volume was divided into 0.5 mm × 0.5 mm × 0.5 mm voxels filling a field-of–view (FOV) of 4 cm (sagittal) × 5 cm (coronal) × 2.5 cm (axial), resulting in a total of 413 100 voxels. The initial MC simulations to calculate the path lengths used 10⁸ photons per source and took 1–3 h for each source on a dual opteron processor. A parallel computing cluster with multiple dual opteron nodes was employed for the computation of the photon path histories for 21 sources within a few hours. Once the history files were computed and stored in memory, the computation time for a recalculation for any absorption was near 15 s for ≈ 2000 S–D pairs. Assuming a homogeneous background and with the reduced scattering coefficient fixed at 10 cm⁻¹, (10) and (11) were used to fit for the S–D coefficients simultaneously with the background absorption of the mouse phantom. The absorption at excitation and emission wavelengths were found to be nearly equal, with values of equation M38

and

Once the background optical properties are determined, the two-point Greens functions from each source and detector to all the medium voxels, viz., equation M40

and

, are calculated using the MC simulation, using the optical properties at λ_x and λ_m, respectively. The CW weight functions for the inversion equation M42

were then calculated in the adjoint formulation by a direct multiplication of the source and detector Green’s functions as in (7).

G. Inverse Problem

The inversion of the fluorophore yield was carried out using Tikhonov regularization, in conjunction with singular value decomposition (SVD). Due to the dramatic variation of the sensitivity profile with depth, a spatially variant preconditioning of the weight matrix was employed [62]. For an under-determined forward problem of the form y = Wx, with W as a general sensitivity matrix, the pseudo-inverse equation M43

defined by equation M44

is given as

(12)

where = WL⁻¹, with L as a diagonal matrix whose diagonal elements are equation M46

, λ is the regularization parameter and α = (max(diag(^T))). Note that since (W^TW)_ij = Σ_l W_liW_lj, we have

(13)

Thus, the diagonal elements of L are the column norms of the sensitivity function W, or the net sensitivity at each voxel summed over all measurement pairs. The normalization by L therefore has the effect of annulling the strong spatial variations in the sensitivity function. (In practice, a small threshold, typically max(L)/100 was added to L to avoid infinities in the normalization.) In order to reduce the computational time for evaluating the inverse in (12) for multiple regularization parameters, a SVD analysis was employed. Writing = USV^T, it can be readily shown that

(14)

For the reconstructions reported here, the data vector y consisted of either the asymptotic amplitudes a_n(r_s, r_d) or the CW components Ū^x,m both of which were renormalized by Θ(r_s, r_d).

IV. Results: Tomography of the Mouse Phantom

The asymptotic lifetime-based algorithm was validated on the TD system using the mouse phantom described in Section III-C. The two inclusions in the phantom were filled with a NIR dye (3,3′-Diethylthiatricarbocyanine) with absorption and emission maxima near 755 and 790 nm. To achieve lifetime contrast, two solutions were prepared with 0.1 mg of the dye dissolved in 50-mL distilled water and glycerol solutions. The lifetimes of the dye in aqueous and glycerol solvents were measured to be near 0.5 and 1 ns, respectively. The increase in fluorescence lifetime with viscosity reflects the reduced ability of the dye molecules in their excited state to relax through rotational modes, which in turn reduces the nonradiative decay rate [74].

The longer lived dye solution in glycerol was injected in inclusion B and the shorter lived aqueous dye solution of the dye in inclusion A (see Fig. 3). The mouse phantom was placed in the imaging setup shown in Fig. 1 and the data were collected in the transillumination geometry for 21 source positions and 93 assigned detectors on the camera image (see Figs. 2 and and4).4). Before dye injection into the two inclusions, the incident fields U^x,x and U^m,m were measured for excitation at 755 and 790 nm, respectively. The fluorescence data was collected after dye injection using a 790–10 nm width band pass filter (CVI – laser), for excitation at 755 nm. Typical values for the other experimental parameters were intensifier gain: 500 V, gatewidth: 300 ps, and CCD integration time: 100 ms. In Fig. 5(a), a representative sample of the DFTR for a single S–D pair is plotted, along with the IRF and the temporal response of the excitation light. The initial time t₀ determined from the IRF is also indicated (see Section III-B).

Fig. 5

(a) Representative TD data of the impulse response (dotted line), diffuse excitation signal, U^x,x at 755 nm (dashed line) and diffuse fluorescence signal U^x,^m at 790 nm (dots) for a single S–D pair. The vertical arrow indicates t₀ determined from (more ...)

The global DFTR (see Section III-E) for a single S–D measurement with the mouse phantom is shown in Fig. 5(b) (inset). The dramatic improvement in the noise statistics of the global DFTR over that for a single SD pair is evident from Fig. 5(b). Although the global DFTR can recover two independent lifetimes, further improvement in the robustness of the fit can be obtained if one of the lifetimes is known in advance, as in a controlled phantom experiment as reported here. This is also possible in practical in vivo applications, if the lifetime of the dye within a control animal can be estimated separately. Considering this scenario, the global fit was performed with the shorter lifetime fixed at 0.5 ns, and the second lifetime was recovered as τ₂ = 0.95 ns, as expected. Fig. 5(b) also shows the raw data and the linear fit for a representative DFTR for a single S–D pair, along with the individual contributions from the shorter and longer lived components.

The S–D coupling coefficients and the background absorption were fit using (10) and (11), as described in Section III-F. A homogenous background absorption was assumed as a first-pass approach, although the formalism presented here can readily incorporate the effect of background absorption due to fluorophores themselves. For the mouse phantom used here, where the top inclusion (A) is located further along the source to detector axis than the bottom inclusion (B) (Fig. 3), the homogenous (or linear) approximation primarily amounts to a neglect of preabsorption of the excitation light to A by any fluorophores present in inclusion B, and the reabsorption of the fluorescence emitted from B by the fluorophores present in A. Since the inclusions are small and the fluorophores are effectively nonscattering, these two effects are expected to only alter the relative intensities of the fluorophore strengths in the inclusions and not their localization (see below), which is the main concern of this paper. Moreover, the actual volume of dye present within the inclusion was not known precisely due to imperfections in the phantom.

The decay amplitudes for the two lifetime components were inverted to recover the fluorescence yield localizations as prescribed in (5). The set of amplitude coefficients of each life-time component for all S–D pairs formed the data vector for the inversion of the yield distributions η₁ (r) and η₂ (r) using (14). Although the Monte-Carlo weight functions were calculated for the full set of 413 000 medium voxels (0.125 mm³ volume) of the original mouse segmentation (see Section III-F), the inversion was performed with a reduced weight matrix of 50 000 voxels (1 mm³) for computational ease. The reconstruction results are shown in Fig. 6 as 3-D images overlayed with the surface boundary of the phantom. The images are iso-surface maps at 95% peak intensity. The regularization parameter λ was near unity for the reconstructions shown and was chosen based on an empirical assessment of image quality. The yield reconstruction for the shorter lifetime of 0.5 ns is seen to have a shallower depth than that for the longer lifetime of 0.95 ns, as expected. Fig. 6(d)–(e) shows the reconstructions performed using the CW component of the entire TD data [(8)]. It is clear that the CW reconstruction fails to resolve the two targets, and is biased towards the inclusion containing the longer lifetime. This bias could be attributed to the fact that the longer lived fluorescence of the dye in the glycerol solution leads to a higher CW intensity than the shorter lived one, since it loses less energy through non-radiative channels [57], assuming identical concentration and excitation conditions. This fact is also evident from the expression for the CW forward model in (8), where the contribution of yield for each lifetime component is scaled by the lifetime. Also, note from the raw time resolved signal in Fig. 5(b) that the shorter lived component has a much weaker contribution to the net fluorescence decay than the longer lived component. We note here that this bias could not possibly arise from any sensitivity variations between the regions occupied by the two inclusions, since they were located approximately symmetrically with respect to the sources and detectors.

Fig. 6

3-D reconstructions of TD fluorescence data using the asymptotic and CW components. (a) Asymptotic lifetime based reconstructions of the DFTR from the mouse phantom injected with the 0.5 ns dye solution in the top inclusion (A in Fig. 3) and the 1 ns (more ...)

The asymptotic TD reconstructions were compared with CT images of the mouse phantom (Section III-C), to verify accuracy. Coregistration was performed between the high-resolution (0.125 mm³) image of the 3-D camera and the CT image, using an affine-transformation with fiducial points chosen from identifiable features on the outline of the images. Once the transformation matrix was determined, the optical images were rotated to the CT-image plane after being first interpolated to revert to the original resolution of the 3-D camera image. The overlayed images in Fig. 7 show that the peak of the asymptotic yield reconstructions match the true locations of the inclusions as recovered by the CT image reasonably well. In future work with phantoms, we will use more well defined fiducials common to the optical and CT images to improve the accuracy of the coregistration and to better quantify the accuracy of the asymptotic reconstructions.

Fig. 7

Co-registration of the optical reconstructions with the CT reconstructed images of the mouse phantom. Two sagittal planes containing the centroid of each inclusion are shown. The CT image is shown in gray scale, with the mouse background visible as dark (more ...)

For comparison, separate measurements were performed with fluorophores occupying only one inclusion at a time, viz., with the glycerol solution occupying B with A empty, and with the aqueous solution occupying A with B empty. Fig. 8 shows a vertical slice of the depth dependent yield for all cases, at the X, Y location of the maxima of the corresponding reconstructed yields. Also shown in the figure are the CT recovered centroids of the inclusions, along the Z-direction. The optical reconstructions are seen to match the locations recovered by CT. It is also seen that the reconstructions with either inclusion A or B filled with the short or the long lifetime fluorophores nearly match the corresponding yield localizations recovered from the combined data set. This indicates that the neglect of background fluorophore absorption has not affected the yield localization in a significant way. However, it is also noted from Fig. 8 (legend) that the strength of the reconstructed yield for inclusion A was approximately a factor of 3 smaller than that for inclusion B, although equal concentrations of dye were used in both. This could be attributed to several factors. The homogeneous assumption could underestimate the yield distribution of inclusion A more than that of B. This is due to the fact that preabsorption of the excitation light by fluorophores in B occurs near the peak of the fluorophore absorption, whereas the reabsorption of the fluorescence from B by the fluorophores in A occurs much further from the absorption maximum. In addition, the presence of glycerol increases the lifetime (as discussed earlier in this section), and could also alter the extinction coefficient, which was not taken into account in this work. Finally, there was also an experimental uncertainty in that the volume of fluorophore solution in each inclusion was not known accurately, due to a slow leakage of the dye from the inclusions back to the injection tubes. A more controlled experiment is needed to verify the quantitation accuracy of the reconstructions and will be carried out in future work.

Fig. 8 also shows that the widths of the asymptotic yield reconstructions (> 1 cm at 80% of the peak) are much larger than the center-to-center separation between inclusions (≈ 6 mm). This illustrates the key advantage of the asymptotic approach: it can resolve lifetime-sensitive targets that are well within the resolution limits imposed by diffuse optics (which determines the widths of the individual distributions). This is made possible since the yield distribution in the asymptotic approach is a distinct image for each lifetime, unlike the CW method which recovers the sum of the two yield distributions. It is also noteworthy that the lateral (i.e., along the length of the mouse, or Y-axis in Figs. 6 and and8)8) center-to-center separation of the inclusions was near 1.6 mm, which is more than a factor of 2 smaller than what was previously achieved (4.3 mm) using free space CW measurements [8].

The asymptotic TD yield reconstructions will be affected by errors in lifetime estimation. An error in estimating lifetimes affects the yield reconstructions indirectly through the corresponding error in the decay amplitudes. The problem of fitting noisy decay data with multiple exponentials is complex and has been dealt with extensively in the literature related to microscopic FLIM (see, for example, [77]–[80]). The global analysis employed in this work offers a robust approach to this problem, since the determination of lifetime is not subject to the highly noisy data sets from the individual measurement pairs. The stability of the global analysis was verified as follows. The biexponential fit for the global analysis was performed with one lifetime having a range of 0.4–0.6 ns, which is ±20% with respect to the true value of the short lifetime component of 0.5 ns. This mimics the scenario where the a priori lifetime is not known accurately. For each case, the second lifetime component was determined from the global fit, and found to be within the range, 0.9–1 ns, which is a ≈ ±6% error with respect to the true value of 0.95 ns. Subsequently, the amplitude coefficients for both lifetimes were determined for each S–D pair using a linear fit as described in Section III-E. Tomography was then performed for each case, and the maximum error in depth localization of the short lived component (inclusion A) was estimated to be 1 mm. The error in the localization of the long lived component (inclusion B) was well within the 1-mm voxel resolution used, and was estimated from the relative variation in the lifetime of this component (±6%) compared to that of the shorter component (±20%) to be near 0.3 mm. These results indicate the robustness of lifetime based tomography using a global lifetime estimation approach, given that a ±20% error bar in lifetime estimation is well above what we can typically achieve experimentally. In future work, we will present a more detailed study of the influence on the lifetime based reconstructions, of various other factors such as the noise statistics of the data collected, the values of the lifetimes involved, their separation, the point at which the asymptotic data is truncated, and the fluorophore distribution.

V. Discussion

We have presented theory and experimental methods for the inversion of fluorescence lifetime and yield distributions within complex turbid media, from TD fluorescence measurements. The application of the imaging system described here for small animal imaging was demonstrated using mouse shaped phantoms, embedded with inclusions for dynamic fluorophore injection. The ability to recover fluorescence lifetime in a direct way rather than indirectly using tomographic inversion, as is necessary in FD approaches, is a unique feature of TD imaging as exploited in this work. This feature is further enhanced by the tomographic algorithms that help separately localize multiple lifetime components within tissue from noninvasive surface measurements. Although the application in the paper demonstrated recovery of ≈ 1-cm-deep targets, the asymptotic model equally well applies to tissue of any thickness, provided the fluorescence lifetime is longer than the tissue absorption timescale [33]. Fluorescence lifetime imaging (FLIM) [65], [77]–[80] is already a well-established microscopy technique that is used to probe lifetime contrast in thin tissue sections. The approach for TD fluorescence tomography from diffuse media as presented in this work serves as a natural extension of time resolved FLIM methods to turbid imaging in thick tissue.

We showed that the asymptotic form of the multiexponential model for diffuse fluorescence, characterized by spatially invariant lifetimes, allows for a global lifetime analysis. This approach is similar to global fitting methods employed for FLIM with thin tissues sections [79], and offers a significant reduction in the computation time while providing a robust way to estimate lifetimes and amplitudes from surface TD measurements on diffuse media.

The approach for TD fluorescence utilized here assumes discrete lifetimes. This assumption is valid in light of the direct observation of these lifetimes in the raw data and their separability using multiexponential fits. The actual constitution of lifetimes within the tissue medium will depend on the specific application of interest. A discrete lifetime model can also accommodate a general scenario where the lifetimes form a continuous distribution. Assuming the long lifetime condition is satisfied for all lifetimes present, the distribution of lifetimes present could be first recovered from the asymptotic (nonexponential) decay using numerical regularization techniques for Laplace inversion, such as the maximum entropy method [75]. The surface amplitude in this case is a distribution for each S–D pair, and could be inverted to obtain the yield distribution for any desired lifetime component.

It should be stressed that an important advantage of fluorescence reconstructions from turbid media is their relative in-sensitivity to errors in estimating the background tissue optical properties [16]. This observation may be readily understood from the second-order nature of the fluorescence excitation-emission process: In the Born approximation, the background Greens functions appear in the fluorescence expression twice, once each for the excitation and the emission processes. Thus, an error in each of the Greens functions due to incorrect background optical properties is propagated to second order (i.e., squared) in the fluorescence sensitivity functions. (Note that direct imaging of intrinsic absorption contrast is much more susceptible to background optical property errors since the measured signal also has first order terms corresponding to the incident field.) On the other hand, it has also been shown that the estimation of fluorescence lifetime as a distribution is more sensitive to optical property estimation errors than the yield distribution [16]. These findings directly support the approach presented here, wherein the lifetime is estimated directly from the raw surface decay data without the need for assumptions or estimations of background properties. The only requirement is that the fluorescence lifetimes be longer than the intrinsic absorption timescale present in the medium, a condition widely satisfied for biological tissue [31]–[33].

We have not made an attempt in this work to compare or quantify the absolute resolution achievable by the TD method against FD and CW approaches. The tomographic analysis presented here was focussed on the decay portion of the TD data, which provides superior localization but suffers from poorer resolution due to the lower SNR compared to the peak of the DFTR. It is noteworthy that even so, the asymptotic portion shows superior performance over the CW reconstructions. A more complete analysis of the TD data that incorporates the early rise and peak portions of the DFTR is likely to improve the resolution of the asymptotic reconstructions [33], and will be dealt with in future work.

Several modifications to the TD system presented here can be envisaged to improve the performance of the system. For instance, in order to fully utilize the tomographic potential offered by the photogrammetric method, multiple view angles could be captured to provide a more complete 3-D surface image of the phantom or animal. The horizontal placement of the imaging subject as used here is readily amenable to this modification. Application of the TD system to living mice is bound to pose additional challenges, such as the influence of breathing and other motion artifacts on the surface boundary acquisition. However, the methodology and preliminary experimental phantom results presented in this work suggest the practical feasibility and advantages of TD fluorescence tomography for small animal imaging. We hope that this work will provide useful guidelines for lifetime based imaging for drug discovery and development in small animals, given the huge variety of biomedically relevant parameters that can alter fluorescence lifetime in tissue, such as pH [57], viscosity [74], oxygen concentration [76], and tissue autofluorescence, as well as molecular interactions such as Forster resonance energy transfer (FRET) [64], [66]. It is also hoped that this work would further motivate the engineering of molecular probes that shift lifetime upon binding. In future work, we will consider in vivo applications of diffuse fluorescence lifetime imaging as presented in this work, in living mouse models of disease.

Acknowledgments

The authors would like to thank A. H. Slocum, Dr. M. K. Kalra, and Dr. R. Gupta (MGH) for help with the CT measurement and reconstructions. The authors would also like to thank Advanced Research Technologies (ART) for loaning the mouse phantom that was used to construct the negative mold.

This work was supported in part by the National Institutes of Health under Grant EB000768, Grant AG026240, and Grant P41-RR14075 and in part by the Advanced Research Technologies.

Footnotes

¹A derivation and discussion of the conditions for the validity of (1)–(3) can be found in [33].

Contributor Information

Anand T. N. Kumar, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Charlestown, MA 02129 USA.

Scott B. Raymond, Harvard-MIT Division of Health Sciences and Technology, Cambridge, MA 02139 USA.

Andrew K. Dunn, Department of Biomedical Engineering, University of Texas–Austin, Austin, TX 78712 USA.

Brian J. Bacskai, Department of Neurology and Harvard Medical School, Massachusetts General Hospital, Charlestown, MA 02129 USA.

David A. Boas, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Charlestown, MA 02129 USA.

References

1. Ntziachristos V, Weissleder R. Experimental three-dimensional fluorescence reconstruction of diffuse media by use of the normalized born approximation. Opt Lett. 2001 Jun;26(12):893–895. [PubMed]

2. Klose AD, Hielscher AH. Fluorescence tomography with simulated data based on the equation of radiative transfer. Opt Lett. 2003;28:1019–1021. [PubMed]

3. Zacharakis G, Ripoll J, Weissleder R, Ntziachristos V. Fluorescent protein tomography scanner for small animal imaging. IEEE Trans Med Imag. 2005 Jul;24(7):878–885. [PubMed]

4. Patwardhan SV, Bloch SR, Achilefu S, Culver JP. Time-dependent whole-body fluorescence tomography of probe bio-distributions in mice. Opt Express. 2005;13:2564–2577. [PubMed]

5. Soubret A, Ripoll J, Ntziachristos V. Accuracy of fluorescent tomography in the presence of heterogeneities: Study of the normalized born ratio. IEEE Trans Med Imag. 2005 Oct;24(10):1377–1386. [PubMed]

6. Soubret A, Ntziachristos V. Fluorescence molecular tomography in the presence of background fluorescence. Phys Med Biol. 2006;51:3983–4001. [PubMed]

7. Deliolanis N, Lasser T, Hyde DD, Soubret A, Ripoll J, Ntziachristos V. Free-space fluorescence molecular tomography utilizing 360 degrees geometry projections. Opt Lett. 2007;32:382–384. [PubMed]

8. Graves EE, Ripoll J, Weissleder R, Ntziachristos V. A sub-millimeter resolution fluorescence molecular imaging system for small animal imaging. Med Phys. 2003 May;30(5):901–911. [PubMed]

9. Schulz RB, Ripoll J, Ntziachristos V. Experimental fluorescence tomography of tissues with noncontact measurements. IEEE Trans Med Imag. 2004 Apr;23(4):492–500. [PubMed]

10. Schulz RB, Peter J, Semmler W, D’Andrea C, Valentini G, Cubeddu R. Comparison of noncontact and fiber-based fluorescence-mediated tomography. Opt Lett. 2006;31:769–771. [PubMed]

11. Montet X, Figueiredo JL, Alencar H, Ntziachristos V, Mahmood UU, Weissleder RR. Tomographic fluorescence imaging of tumor vascular volume in mice. Radiology. 2007;242:751–758. [PubMed]

12. O’Leary MA, Boas DA, Li XD, Chance B, Yodh A. Fluorescence lifetime imaging in turbid media. Opt Lett. 1996 Jan;21(2):158–160. [PubMed]

13. Paithanker D, Chen A, Pogue B, Patterson M, Sevick-Muraca EM. Imaging of fluorescent yield and lifetime from multiply scattered light re-emitted from tissues and other random media. Appl Opt. 1997;36:2260–2272. [PubMed]

14. Eppstein MJ, Hawrysz DJ, Godavarty A, Sevick-Muraca EM. Three-dimensional, bayesian image reconstruction from sparse and noisy data sets: Near-Infrared fluorescence tomography. Proc Nat Acad Sci. 2002 Jul;99(15):9619–9624. [PubMed]

15. Lee J, Sevick-Muraca EM. Three-dimensional fluorescence enhanced optical tomography using referenced frequency-domain photon migration measurements at emission and excitation wavelengths. J Opt Soc Amer. 2002 Apr;19(4):759–771. [PubMed]

16. Milstein AB, Oh S, Webb KJ, Bouman CA, Zhang Q, Boas DA, Millane RP. Fluorescence optical diffusion tomography. Appl Opt. 2003 Jun;42:3081–3094. [PubMed]

17. Godavarty A, Eppstein MJ, Zhang CY, Theru S, Thompson AB, Gurfinkel M, Sevick-Muraca EM. Fluorescence-enhanced optical imaging in large tissue volumes using a gain modulated ICCD camera. Phys Med Biol. 2003;48:1701–1720. [PubMed]

18. Thompson AB, Hawrysz DJ, Sevick-Muraca EM. Near-infrared fluorescence contrast-enhanced imaging with area illumination and area detection: The forward imaging problem. Appl Opt. 2003;42(19):4125–4136. [PubMed]

19. Joshi A, Bangerth W, Sevick-Muraca EM. Non-contact fluorescence optical tomography with scanning patterned illumination. Opt Express. 2006;14:6516–6534. [PubMed]

20. Eppstein MJ, Fedele F, Laible J, Zhang C, Godavarty A, Sevick-Muraca EM. A comparison of exact and approximate adjoint sensitivities influorescence tomography. IEEE Trans Med Imag. 2003 Oct;22(10):1215–1223. [PubMed]

21. Godavarty A, Sevick-Muraca EM, Eppstein MJ. Three-dimensional fluorescence lifetime tomography. Med Phys. 2005;32:992–1000. [PubMed]

22. Roy R, Thompson AB, Godavarty A, Sevick-Muraca EM. Tomographic fluorescence imaging in phantoms: A novel reconstruction algorithm imaging geometry. IEEE Trans Med Imag. 2005 Feb;24(2):137–154. [PubMed]

23. Sahu AK, Joshi A, Kupinski MA, Sevick-Muraca EM. Assessment of a fluorescence-enhanced optical imaging system using the hotelling observer. Opt Express. 2006;14(17):7642–7660. [PMC free article] [PubMed]

24. Roy R, Godavarty A, Sevick-Muraca EM. Fluorescence-enhanced optical tomography of a large tissue phantom using point illumination geometries. J Biomed Opt. 2006;11:044007. [PubMed]

25. Rasmussen JC, Joshi A, Pan T, Wareing T, McGhee J, Sevick-Muraca EMEM. Radiative transport in fluorescence-enhanced frequency domain photon migration source. Med Phys. 2006;33:4685–4700. [PubMed]

26. Roy R, Godavarty A, Sevick-Muraca EM. Fluorescence-enhanced three-dimensional lifetime imaging: A phantom study. Phys Med Biol. 2007;52:4155–4170. [PubMed]

27. 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–292.

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

29. Hattery D, Chernomordik V, Loew M, Gannot I, Gandjbakhche A. Analytical solutions for time-resolved fluorescence lifetime imaging in a turbid medium such as tissue. J Opt Soc Amer. 2001;18:1523–1530. [PubMed]

30. 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–2260. [PubMed]

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

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

33. Kumar ATN, Raymond SB, Boverman G, Boas DA, Bacskai BJ. Time resolved fluorescence tomography based on lifetime contrast. Opt Express. 2006;14:12255–12270. [PMC free article] [PubMed]

34. Hassan M, Riley J, Chernomordk V, Smith P, Pursley R, Lee SB, Capala J, Gandjbakhche AH. Fluorescence lifetime imaging system for in vivo studies. Molec Imag. 2007;6:229–236. [PMC free article] [PubMed]

35. Soloviev VY, McGinty J, Tahir KB, Neil MAA, Sardini A, Hajnal JV, Arridge SR, French PMW. Fluorescence lifetime tomography of live cells expressing enhanced green fluorescent protein embedded in a scattering medium exhibiting background autofluorescence. Opt Lett. 2007;32:2034–2036. [PubMed]

36. Folli S, Westermann P, Braichotte D, Pelegrin A, Wagnieres G, van de Bergh H, Mach JP. Antibody-indocyanin conjugates for immunophotodetection of human squamous cell carcinoma in nude mice. Cancer Res. 1994;54:2643–2649. [PubMed]

37. Jain RK. Barriers to drug delivery in solid tumors. Sci Am. 1994;271:58–65. [PubMed]

38. Ballou B, Fisher GW, Waggoner AS, Farkas DL, Reiland JM, Jaffe R, Mujumdar RB, Hakala TR. Tumor labeling in vivo using cyanine-conjugated monoclonal antibodies. Cancer Immunol Immunother. 1995;41:257–263. [PubMed]

39. Neri D, et al. Targeting by affinity-matured recombinant antibody fragments of an angiogenesis associated fibronectin isoform. Nat Biotechnol. 1997;15:1271–1275. [PubMed]

40. Goldsmith SJ. Receptor imaging: Competitive or complementary to antibody imaging. Seminard Nucl Med. 1997;27:85–93. [PubMed]

41. Tung CH, Bredow S, Mahmood U, Weissleder R. Preparation of a cathepsin D sensitive near-infrared fluorescence probe for imaging. Bioconjugate Chem. 1999 Sep 20;10(5):892–896. [PubMed]

42. Weissleder R, Tung CH, Mahmood U, Bogdanov A., Jr In vivo imaging of tumors with protease-activated near-infrared fluorescent probes. Nature Biotechnol. 1999;17:375–378. [PubMed]

43. Mahmood U, Tung CH, Bogdanov A, Weisslede R. Near-infrared optical imaging of protease activity for tumor detection. Radiology. 1999;213:866–870. [PubMed]

44. Becker A, Riefke B, Ebert B, Sukowski U, Rinneberg H, Semmler W, Licha K. Macromolecular contrast agents for optical imaging of tumors: Comparison of indotricarboyanine-labeled human serum albumin and transferrin. Photchem Photbio. 2000;72:234–241. [PubMed]

45. Achilefu S, Dorshow RB, Bugaj JE, Rajagopalan R. Novel receptor targeted fluorescent contrast agents for in vivo tumor imaging. Invest Radiol. 2000 Aug 1;35(8):479–485. [PubMed]

46. Soukos NS, Hamblin MR, Keel S, Fabian RL, Deutsch TFTF, Hasan T. Epidermal growth factor receptor-targeted immunophotodiagnosis and photoimmunotherapy of oral precancer in vivo. Cancer Res. 2001;61:4490–4496. [PMC free article] [PubMed]

47. Zaheer A, Lenkinski RE, Mahmood A, Jones AG, Cantley LC, Frangioni JV. In vivo near-infrared fluorescence imaging of osteoblastic activity. Nat Biotechnol. 2001;19:1148–1154. [PubMed]

48. Bremer C, Tung CH, Weissleder R. In vivo molecular target assessment of matrix metalloproteinase inhibition. Nat Med. 2001;7:743–748. [PubMed]

49. Becker A, Hessenius C, Licha K, Ebert B, Sukowski U, Semmler W, Wiedenmann B, Grotzinger C. Receptor-targeted optical imaging of tumors with near-infrared fluorescent ligands. Nat Biotech. 2001;19:327–331. [PubMed]

50. Bugaj JE, Achilefu S, Dorshow RB, Rajagopalan R. Novel fluorescent contrast agents for optical imaging of in vivo tumors based on a receptor-targeted dye-peptide conjugate platform. J Biomed Opt. 2001 Apr;6(2):122–133. [PubMed]

51. Achilefu S, Jimenez HN, Dorshow RB, Bugaj JE, Webb EG, Wilhelm RR, Rajagopalan R, Erion JJJL. Synthesis, in vitro receptor binding, and in vivo evaluation of fluorescein and carbocyanin peptide-based optical contrast agents. J Med Chem. 2002;45:2003–2015. [PubMed]

52. Hintersteiner M, et al. In-vivo detection of amyloid-beta deposits by near-infrared imaging using an oxazine-derivative probe. Nat Biotechnol. 2005;23:577–583. [PubMed]

53. Ntziachristos V, Ripoll J, Wang LV, Weissleder R. Looking and listening to light: The evolution of whole-body photonic imaging. Nature Biotechnol. 2005;23:313–320. [PubMed]

54. Ye YP, Bloch S, Xu BG, Achilefu S. Design, synthesis, and evaluation of near infrared fluorescent multimeric RGD peptides for targeting tumors. J Med Chem. 2006;49:2268–2275. [PMC free article] [PubMed]

55. Wang W, Ke S, Kwon S, Yallampalli S, Cameron AG, Adams KE, Mawad ME, Sevick-Muraca EM. A new optical and nuclear dual-labeled imaging agent targeting interleukin receptor alpha-chain. Bioconj Chem. 2007;18:397–402. [PubMed]

56. Wang W, Davis EA, Ke S, Kwon S, Davis AR, Lazard ZW, Cameron A, Sevick-Muraca EM. Multiple target-specific peptides imaging heterotopic bone re-generation process. Biopolymers (Biospectroscopy) 2007;88:527–1527.

57. Lakowicz JR. Principles of Fluorescence Spectroscopy. 2. Norwell, MA: Kluwer; 1999.

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

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

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

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

62. Culver JP, Durduran T, Furuya D, Cheung C, Greenberg JH, Yodh AG. Diffuse optical tomography of cerebral blood flow, oxygenation and metabolism in rat during focal ischemia. J Cereb Blood Flow Metabol. 2003;23:911–924. [PubMed]

63. Agrawal GP. Nonlinear Fiber Optics. 4. San Diego, CA: Academic; 2007.

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

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

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

67. Arridge SR. Optical tomographyin medical imaging. Inverse Prob. 1999;15:R41–R93.

68. Boas DA. PhD dissertation. Univ. Pennsylvania; Philadelphia, PA: 1997. Diffuse photon probes of structural and dynamical properties of turbid media: Theory and biomedical applications.

69. Arridge SR, Lionheart WRB. Nonuniqueness in diffusion-based optical tomography. Opt Lett. 1998;23:882–884. [PubMed]

70. Wang L, Jacques SL, Zheng L. MCML-Monte carlo modeling of light transport in multi-layered tissues. Comput Meth Programs Biomed. 1995;47:131–146. [PubMed]

71. Hayakawa CK, Spanier J, Bevilacqua F, Dunn AK, You JS, Tromberg BJ, Venugopalan V. Perturbation monte carlo methods to solve inverse photon migration problems in heterogeneous tissues. Opt Lett. 2001;26(17):1335–1337. [PubMed]

72. Boas DA, Culver JP, Stott JJ, Dunn AK. Three dimensional monte carlo code for photon migration through complex heterogeneous media including the adult human head. Opt Express. 2002;10:159–170. [PubMed]

73. Boas DA, Gaudette TJ, Arridge SR. Simultaneous imaging and optode calibration with diffuse optical tomography. Opt Express. 2001;8:263–270. [PubMed]

74. Taylor JR, Adams MC, Sibbett W. Investigation of viscosity dependent fluorescence lifetime using a synchronously operated picosecond streak camera. Appl Phys. 1980;21:13–17.

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

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

77. James DR, Ware WR. A fallacy in the interpretation of fluorescence decay parameters. Chem Phys Lett. 1985;120:455–459.

78. Siegel J, Elson DS, Webb SED, Gambaruto GL, Lévêque-Fort S, Stamp GWH, Wallace AL, Alvarez AF, French PMW. Studying biological tissue with fluorescence lifetime imaging: Microscopy, endoscopy, and complex decay profiles. Appl Opt. 2003;42:2995–3004. [PubMed]

79. Pelet S, Previte MJR, Lalho LH, So PTC. A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation. Biophys J. 2004;87:2807–2817. [PubMed]

80. Jo JA, Fang Q, Papaioannou T, Baker JD, Dorafshar AH, Reil T, Qiao J, Fishbein MC, Freischlag JA, Marcu L. La-guerre-based method for analysis of time-resolved fluorescence data: Application to in-vivo characterization and diagnosis of atherosclerotic lesions. J Biomed Opt. 2006;11:021004–1. [PMC free article] [PubMed]