Applying time-dependent data for fluorescence tomography

Can time-resolved, high-resolution data as acquired by an intensified gated CCD camera (ICCD) aid in the tomographic reconstruction of fluorescence concentration? Usually it is argued that fluorescence is a linear process and thus does not require non-linear, time-dependent reconstructions algorithms, unless absorption and scattering coefficients need to be determined as well. Furthermore, the acquisition of a number of time frames is usually prohibitive for fluorescence measurements, at least in small animals, due to the increased total measurement time. On the other hand, it is obvious that diffusion is less pronounced in images at early gates, due to selective imaging of photons of lower scatter order. This will be the case also for photons emitted by fluorescent sources. Early-gated imaging might increase the contrast in acquired images and could possibly improve fluorescence localization. Herein, we present early gated fluorescence images obtained from phantoms and compare them to continuously acquired data. Increased contrast between background and signal maximum can be observed in time-gated images as compared to continuous data. To make use of the properties exhibited by early gated frames, it is necessary to use a modified reconstruction algorithm. We propose a variant of the well-known Born approximation to the diffusion equation that allows to take into account single time frames. The system matrix for the time-dependent Born approach is more complex to calculate, however the complexity of the actual inverse problem (and the acquisition times) of single-frame reconstructions remains the same as compared to continuous mode.


INTRODUCTION
In recent years, innovations in photomultiplier design have led to the development of gated intensified CCD camera systems with a gate width of only a few hundred picoseconds. These systems have been long awaited by the biomedical optics community, as they allow the acquisition of time-resolved, high-resolution datasets for low-light applications. Furthermore, they enable high-resolution acquisitions of early photons, which have undergone only a few scattering events and thus can be used for direct reconstruction based on line integrals (i.e., the Radon transform) instead of diffusion 1 . While it has been shown in the past that time-dependent information (may it be time-domain or frequencydomain data) is necessary to resolve both absorption and scattering coefficients in turbid media 2 , it is an open question whether or not this type of information could help in improving reconstructions of fluorescence data. Fluorescence emission of lowly concentrated dyes is considered a linear process, i.e., it is usually modeled using a linear differential equation technique, based on the Green functions of L, for the estimation of q, making reconstruction a linear inverse problem with no specific need for time-resolved information. For this reason, it is questionable whether much additional information can be gained using time-gated imaging.
On the other hand, it can be expected that in early gates, photons of a lower scatter order will be more pronounced and that fluorescent inclusion will therefore appear less blurred, increasing the dynamic content of acquired images and hopefully increasing spatial resolution and / or sensitivity.
Herein we present experimentally acquired datasets indicating that structural differences can be observed between fluorescence data of early time gates and time-independent (CW) fluorescence measurements, leading to a higher signalto-background ratio. These differences can be mainly attributed to reduced diffusion of fluorescence photons within the observed early gate. To enable 3D reconstructions of single-frame time-dependent measurements, we propose a linear, diffusion-based algorithm derived from the well-known Born approach, based on the time-dependent diffusion equation. While this algorithm is of higher computational complexity than the continuous approach, total measurement times are comparable to CW measurements. The reconstruction algorithm is subsequently applied to simulated data sets.

EXPERIMENTAL SETUP AND DATA ANALYSIS
The experimental setup is similar to a non-contact fluorescence tomography setup described previously 3 and is depicted in Fig. 1. The output beam of a pulsed laser diode (80MHz repetition rate, emission at 650nm, PDL800 by PicoQuant, Berlin, Germany) was collimated onto a cylindrical phantom (4cm diameter, 10cm height) having tissue-like absorption and scattering properties (µa = 0.1cm -1 , µs' = 10cm -1 ). The phantom contained a number of cylindrical inclusions (with 2mm, 4mm and 6mm diameter) along its main axis, see Fig. 2. According to the actual specifications given below the inclusions were filled with an intralipid / ink solution to mimic the background properties, with some of the inclusions additionally containing different concentrations of Cy5.5 fluorochrome. For non-contact detection, an intensified gated CCD (ICCD) camera (PicoStar HR, LaVision, Göttingen, Germany) was coupled to a wide-angle objective (CNG, Schneider, Kreuznach, Germany) with high numerical aperture (NA = 1.4). The camera system was placed at a distance of 6.2cm from the phantom's main axis. The ICCD was triggered by the diode driver using a proprietary physical delay line. Optical filtering of the fluorescence signal was performed using a combination of glass filters (OG695 and OG715, Schott, Germany). To yield different positions of the incident laser source, the phantom was rotated around its main axis in 10° steps, resulting in 36 different projections. A schematic drawing of this setup can be found in fig. 1.

Pulse Generator
Laser Diode

Delay Generator
Gated ICCD Camera Phantom Figure 1: Experimental Setup. The collimated output beam of a pulsed laser diode excites fluorochromes in a turbid tissue phantom. Fluorescence signal is acquired using a gated intensified CCD camera (details see text).
For each experiment, six different time gates (300ps width) were acquired per projection, triggered at 500ps, 700ps, 900ps, 1100ps, 1500ps, and 2500ps after arrival of the optical pulse at the camera. Additionally, the image intensifier was also used in continuous mode to yield CW data. Exposure times were chosen between 10s and 30s for gated fluorescence acquisition, 2s for continuous acquisition of fluorescence, and 0.1s for the acquisition of images of the excitation light (i.e., not using any filters). Depending on the gate delay, the gain voltage of the image intensifier was chosen to yield optimal signal amplitude. The gain varied between different delays but was kept constant for all projections. Dark noise was eliminated by subtracting background images.
The influence of filter leakage was calculated by using the earliest gate acquisition which was expected to contain nearly no fluorescence due to fluorescence lifetime effects. The excitation image was scaled to optimally fit the early-gated fluorescence images in all projections. A scaling factor of q = 2·10 -³ was determined, such that the actual fluorescence image f i used for data analysis was obtained as i i i f m qx , where m and x denote the acquired images of the emission and excitation signal, respectively, at projection i. To quantify the difference between images of time-gated and continuously acquired images we analyzed image contrast of time-gated data with respect to gate delay and in comparison to the contrast visible in CW data. Image contrast c i of a diffuse projection image i was defined as max min / max min where min i and max i denote minimum and maximum pixel values in the image. As diffuse projection images are expected to change pixel values only very smoothly, the contrast between neighboring pixels instead of overall image contrast could be chosen as an estimator for pixel noise.

PROPOSED RECONSTRUCTION MODEL
for a source at r s and a point r m inside the medium.
Due to the time convolution, eq. (2) is of higher computational complexity; however, the matrix required to determine n is of equal size if only one time step t d is to be used. Computational complexity increases linearly with the number of time-frames that have to be computed. Using a time-stepping of 20ps, for a gate at 700ps, at least 35 steps have to be computed, yielding an increase in computational load of ~35 as well.
The normalized Born method for continuous wave mode was introduced to abolish source coupling and detector efficiency , s d by division of the fluorescence signal with the according signal obtained without a filter 4 . For the model of eq. (2), as the efficiency factor is contained in the equation independent of the time frame, one could normalize with the according frame of the transmitted light, with a different time frame, or simply with a continuous wave, timeindependent image, depending on whether or not , s d changes in between measurements. If only source coupling is of importance, it is probably advisable to normalize with the time-independent transmitted data, as these images have the best signal-to-noise ratio. However, we did not make use of normalization in the experiments and simulations conducted herein and just used the data as is.
The Green functions of eq. (2) were calculated using a time dependent finite element model, implemented in C++ using the deal.II library 5 and using the grid depicted in fig. 2. This grid was created using the CuBIT meshing system (Sandia National Laboratories, New Mexico, USA). The mesh was optimized to yield voxels of 1.5mm side length, limiting the resolution but decreasing computation times. The convolution was performed by transforming the time series of the Green functions for each voxel into Fourier space and multiplying the three functions there. Our implementation used the freely available FFTW3 library 6 . Inversion of the resulting sensitivity matrix was performed using the well-known algebraic reconstruction method 3,4 with 100 iterations and a relaxation parameter of = 0.1.

RESULTS
Changes in image contrast and signal shape are visible in all experimental setting performed. Especially large is the increase in the setup with two 4mm tubes (fig 2b; the resulting images for a number of time frames along with appropriate profiles are depicted in fig. 3). As expected, especially the early frames show a significantly different signal as compared to the continuously acquired images. However, at 1500ps this effect diminishes. Subsequently, we performed a contrast analysis as described in section 2. To compare the model we used with experimental data, we chose to simulate the signal by analytically modeling the cylindrical inclusions and performing a diffusion-based forward computation. The contrast curves obtained from simulated data are depicted in figure 4, while the experimentally acquired data is given in figure 5. It can be noted that for the experimentally acquired data as well as for the simulated data, a contrast increase can be noted at early gates which is safely above noise level. The actual increase in contrast depends on the situation. For a setting where a big inclusion is located in the center of the phantom (see fig. 2c, fig. 4c, and fig. 5c), a strong fluorescence signal is visible in all projection, increasing the contrast in all projections for early gates, while a setting where the fluorescence signal is sometimes at the opposite of the camera, leading to a strongly diffused signal, the contrast increase is not so high but reaches a level of about 1.2-1.5 above the contrast in the continuous data, which is still significantly above the noise level. However, the simulated and the observed data show some differences in their contrast behavior, which might be due to fluorescence lifetime effects or model inconsistencies.
To examine possible changes in the quality of reconstructed images, we reconstructed the simulated data using the approach described in section 3 for a weight matrix considering a gate at 700ps after source pulse arrival, 36 projections at 10° stepping, and a non-contact detector grid spanning an area of 2cm along the axis and 120° around the cylinder, having 3x19 detectors, as described previously 3 . Data for all three settings (a)-(c) of figure 2 were calculated and inverted using algebraic reconstruction (100 iterations, = 0.1). Central transversal slices of the reconstructed datasets are depicted in figure 5 for the time-independent, and in figure 6 for the time-gated data, respectively.
While only in the setting with two 4mm tubes the inclusions can be separated well, differences are visible between timegated and continuous approach in the other settings as well: the six inclusions appear in a better defined triangular shape, while the 6mm inclusion surrounded by a 2mm tube becomes more eccentric in the direction of the small inclusion. The difference in image quality is visible already at earlier iterations of the matrix inversion (data not shown) and thus indicates that the time-gated system offers a less ill-posed problem.

DISCUSSION
The data presented in here shows systematical differences between gated images of fluorescence, acquired at an early time and continuously acquired data, with gated images exhibiting higher contrast (while having slightly worse signal-tonoise characteristics). These differences indicate that early fluorescence images are less diffuse and thus might increase the achievable spatial resolution in the reconstruction process. While the standard Born approach does not support timedependent information, we present a similar method allowing a single, early time frame to be used. The forward model is more complex to establish, but the resulting inverse problem has the same complexity (matrix size) as the continuous model. Furthermore, the time needed to acquire experimental data does not increase as compared to continuous acquisitions, as only a single time gate has to be measured.
In reconstructions of simulated data sets, it could be seen that the resulting sensitivity matrix of the time-gates approach yields results in better and faster convergence of the data. However, the overall resolution was in both cases far below 2mm, due to the coarse mesh used to minimize computation times. In a next step, we will try to reconstruct the experimentally acquired data with that model. This will require some adaptation of the model to compensate for source pulse width, and gate width of the intensifier.
Furthermore, we proposed normalization methods similar to the normalized Born approach known from continuous wave imaging to eliminate source-detector coupling issues. These approaches need further examination, in terms of the effects that a normalization with continuously acquired data might have on the images.