|Home | About | Journals | Submit | Contact Us | Français|
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.
Fluorescence diffuse optical tomography – 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 –.
Fluorescence tomographic imaging can be broadly classified based on the excitation mechanism employed: continuous-wave (CW) methods use steady state excitation sources –, frequency-domain (FD) methods use frequency modulated (~ MHz) sources –, and time-domain (TD) methods – 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 –, , . 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 . 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 , . 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 . We note that FD techniques – and derived data types from TD measurements such as the Fourier , Laplace , or Mellin transforms (i.e., moments)  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 –, 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 –, , . 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.
A theory of TD fluorescence in turbid media was recently presented . 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 rd and time t, for an impulsive excitation at wavelength λx, position rs and t = 0, can be expressed as a multiexponential sum
It should be noted that the decay amplitudes An depend on time, in addition to the source and detector locations. Within the Born approximation, which assumes a single fluorescence emission event, the An’s define a linear inverse problem for the yield distributions of each lifetime component1
The Green’s functions used in (3) are simply the time-dependent background Green’s functions for light propagation (described by the transport or diffusion quations ) at the excitation and emission wavelengths, Gx and Gm, but evaluated with reduced absorptions, 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 and on the respective total background absorptions and . The background absorptions can in general include the unknown fluorophore absorption in a nonlinear treatment. Also, 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)−1 , , ), the TD expression in (1) smoothly goes over to an asymptotic limit where the amplitudes An become time independent . In this limit, the fluorescence decay factorizes into spatial and temporal components
The amplitudes, an, are the long time values of An, obtained by extending the upper limit of the time-integral on the right-hand side of (2) to ∞ with the understanding that the has completed its evolution and is constant in the asymptotic region 
where the CW weight functions are given by
where and are the CW counterparts of and . The weight functions 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 an 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.
where B are the standard CW Born-sensitivity functions for an absorption perturbation . 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.
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 . 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.
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.
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 t0. The time t0 refers to the time when the excitation pulse is incident on the surface of the imaging medium. A correct estimate of t0 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 . An advantage of noncontact detection is that t0 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 t0 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.
Note that the time t0 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 t0, due to the finite thickness of the mouse. Before analyzing the fluorescence data, t0 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
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.
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 TiO2 (paint) and ink was added to the resin, to achieve approximate background optical properties of and μa = 0.1 cm−1 . 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.
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.
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 , 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.
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.
The computation of the forward problem was carried out using a transport based approach, via the Monte-Carlo (MC) method , . 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 Ux,x and Um,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 , 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−1, 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 . It was also assumed that 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 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 . If the MC generated CW Greens’s functions for a pair of optical absorptions and are and , the incident signal at λx and λm can be written as
where Θ(rs, rd) is an S–D dependent coupling strength and , 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 . 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 , , 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 108 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−1, (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 and .
Once the background optical properties are determined, the two-point Greens functions from each source and detector to all the medium voxels, viz., and , are calculated using the MC simulation, using the optical properties at λx and λm, respectively. The CW weight functions for the inversion were then calculated in the adjoint formulation by a direct multiplication of the source and detector Green’s functions as in (7).
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 . For an under-determined forward problem of the form y = Wx, with W as a general sensitivity matrix, the pseudo-inverse defined by is given as
where = WL−1, with L as a diagonal matrix whose diagonal elements are , λ is the regularization parameter and α = (max(diag(T))). Note that since (WTW)ij = Σl WliWlj, we have
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 = USVT, it can be readily shown that
For the reconstructions reported here, the data vector y consisted of either the asymptotic amplitudes an(rs, rd) or the CW components Ūx,m both of which were renormalized by Θ(rs, rd).
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 .
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 Ux,x and Um,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 t0 determined from the IRF is also indicated (see Section III-B).
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 τ2 = 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 η1 (r) and η2 (r) using (14). Although the Monte-Carlo weight functions were calculated for the full set of 413 000 medium voxels (0.125 mm3 volume) of the original mouse segmentation (see Section III-F), the inversion was performed with a reduced weight matrix of 50 000 voxels (1 mm3) 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 , 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.
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 mm3) 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.
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 .
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, –). 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.
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 . Fluorescence lifetime imaging (FLIM) , – 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 , 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 . 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 . 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 . 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 –.
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 , 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 , viscosity , oxygen concentration , and tissue autofluorescence, as well as molecular interactions such as Forster resonance energy transfer (FRET) , . 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.
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.
1A derivation and discussion of the conditions for the validity of (1)–(3) can be found in .
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.