1. Introduction

Fluorescence molecular tomography (FMT) can volumetrically image and three-dimensionally resolve molecular function by reconstructing the distributions of the fluorescence molecular probes [1–6]. In FMT, the data can be acquired as a steady-state intensity with continuous wave (CW), or as a complexity of measurable amplitude and phase (frequency-domain, FD), or as a time-varying intensity the system responses to an ultra-short input pulse (time-domain, TD). The image reconstruction approaches in steady-state model have been fairly demonstrated [7–11]. For example, Ntziachristos et al used the CW to image the tumors targeted by special fluorescence molecular probes in mice. The fluorescence lifetime image reconstructed by frequency-domain and time-domain models can provide further information about the surroundings of the biological fluorophore, such as pH, enzyme, oxygen, etc. [12–16]. For lifetime imaging, most of the reported works employed frequency-domain lifetime model, and a number of simulative and experimental methods have been demonstrated [7–15]. Boas et al proposed Bayesian image reconstruction for 3D fluorescence lifetime tomography using experimental measurement [8]. Yuan and Godavarty et al separated the imaging procedure into two steps to reconstruct the fluorescent yield and lifetime respectively [11–12]. In addition, multiple-frequency data were used to improve image quality by A. B. Milstein et al [13].

Comparing with CW and FD, time-domain has some advantages: the fluorescence yield and lifetime can recover simultaneously and the analysis of multiple components can be analyzed in a direct way [16]. In addition, because typical fluorescence lifetime is very short, to increase phase shift needs higher modulation frequency, but the frequency up to 1GHz and higher is not technically achievable yet as well as not practical [17], so time-domain system is another option.

Time-domain system employs an ultra-short pulsed light for the excitation and a time-resolved device, such as a streak camera or a time-correlated single photon counting system, to detect the time-dependent fluorescent light. The reconstruction based on analytical solution has been studied with simulative data by Lam et al [18], while numerical methods are more practical and primarily based on the full time-resolved scheme or feature-data one [19–23]. Although the time-resolved scheme provides the richest information on photon migration in tissue, it is generally impractical to the realistic three dimensional cases due to a large computation cost. Instead, some featured data types are usually employed for characterizing the temporal profile and saving the computation time. Hillman et al have proposed a Mellin-transform method. However, many practices have revealed that high-order Mellin’s transforms are rather sensitive to noise [21]. Therefore, both the Laplace and the Fourier transforms are more preferably used to analyze the time-resolved signals [19, 23]. It has been pointed out that the Laplace transform is more practically suitable since the periodic nature of the Fourier transform usually results in strong oscillation in the ratios [24].

Because photon migration in tissue is inherently three dimensional (3D) and the fluorescence lifetime usually a meaningful indicator of biochemical environment in tissue, developing the three-dimensional time-domain FMT is necessarily important. In this paper, we briefly describe a finite element method (FEM) solution to the Laplace-transformed time-domain coupled diffusion equations at first. Then a linear inversion scheme for simultaneous reconstruction of both the fluorescent yield and lifetime images are proposed within the normalized Born-ratio formulation. A variety of 3D numerical phantoms are devised to assess the algorithm performance. With assumption of the Gaussian noise model, we finally investigate the influence of the transform-factors on the image quality and find that an appropriate choice of the transform-factors might considerably increase the noise-robustness of the reconstruction.

2. Methodology

2.1 Forward model

In time-domain FMT, when an ultra-short pulsed laser illuminates tissue surface, the excitation wave propagating in scattering-absorbing medium is highly absorbed by both exogenous and endogenous chromophores and fluorophores. Then the excitation and fluorescent emission lights propagating through tissue and can be measured by a time-resolved system at the boundary. The propagation of excitation and emission lights can be mathematically modeled by the well-known time-domain coupled diffusion equations shown as Eq. (1)

where the subscripts x and m denote the excitation and emission wavelengths respectively; c is the speed of light in medium; the optical properties involved are the absorption coefficient μ_aν(r)(ν∈[x,m]), the reduced scattering coefficient μ′_sν(r)(mm^-1) and the diffusion coefficient κ_ν(r,t)=c/[3μ′_sν(r)]; the fluorescent yield ημ_af(r) (product of the quantum efficiency and absorption coefficient of the fluorophore) and lifetime τ(r) are the fluorescence properties. In Eq. (1), the source term of the excitation equation, i.e., S_x(r,r _s,t)=δ(r-r _s,t), is an impulse function at excitation position r _s=1/μ′_sx and t=0, while the source term of the emission equation, i.e., s_m(r,r _s,t)=ημ_af(r)∫^t ₀ cΦ_x(r,t′)Γ(r)e^{[-(t-t′)Γ(r)]}dt′(Γ(r)=1/τ(r)), is originated from the propagation of the excitation light in fluorescent target, and might be distributed throughout the medium.

By Laplace-transforming the above time-domain coupled diffusion equations with a complex transform-factor β, the following steady-state coupled equations are obtained [16]

where Φ _ν(r,r _s,β) is the Laplace transform of the time-dependent photon density Φ_ν(r,r _s,t). The measurable flux I _ν(r _d,r _s,β) at the boundary can be simply calculated by the Fick’s law under Robin boundary condition, i.e., I_ν(r _d,r _s,β)=cΦ_ν(r _d,r _s,β)(1-R_f)/[2(1+R_f)], where R_f≈0.53 represents the internal reflection coefficient at the air-tissue boundary. In this paper, the coupled diffusion equation is solved using a finite element method with the support of the extrapolated boundary condition [16].

2.2 Inverse issue

The image reconstruction in FMT, i.e., the inverse issue, is to recover the three-dimensional map of the fluorescent property from the boundary measurements at emission wavelength. In this work the normalized Born-ratio formulation is employed to improve image quality for its advantages of independence of source intensity and less sensitivity to systematic errors [25]. With this approach, the emission flux is normalized with the excitation one measured at boundary site r _d with the excitation source at r _s, and an integral equation is obtained according to Eq. (2)

where x(r,β)=ημ_af(r)/[1+βτ(r)]; G_m(r _d,r,β) is the Green function of emission diffusion equation that describes the propagation of the emission wave emitted from a fluorochrome at r to a detector at r _d, and complies with the reciprocity theorem.

To solve the inverse problem, the volume Ω is discretized into E elements joined at N vertex nodes by finite element method. Accordingly, the volume integral in Eq. (3) reduces to the following

where V_e corresponds to a particular voxel in the discretized geometry. Ḡ_m(r _d,r,β), $\overline{\Phi }$ _x(r,r _s,β) and x̄(r,β) denote mean values in V_e. Eq. (4) is a set of linear equations that can be expressed in matrix form

where I^nb(β)=[I^nb(r ¹ _s,r ¹ _d,β), I^nb(r ^{2 s},r ¹ _d,β),…, I^nb(r ^S_s,r ^D_d,β)]^T, r ⁱ_s and r ^j_d denote the ith source and the jth detector respectively. The unknown quantity x(r,β)≈x ^T(β)u(r), where x(β)=[x ₁(β), x ₂(β), …x_N(β)]^T and u(r)=[u ₁(r),u ₂(r),…,u_N(r)]^T is the shape function. W(β)is a matrix of S×D rows and N columns. Any element in the matrix can be given as

where V_e numerates all the elements joined at the nth node; ${\overline{G}}_m^{\left(V_e\right)}({\mathbf{r}}_d^i,\beta )$ and ${\overline{\Phi }}_x^{\left(V_e\right)}({\mathbf{r}}_s^j,\beta )$ are the mean values of G_m(r ⁱ_d,r,β) and Φ_x(r,r ⁱ_s,β) at the nodes of the element respectively; Eq. (5) is usually of large scale and ill-posed. The algebraic reconstruction technique (ART) is used to solve the inverse problem. This method generates reconstructions via an iterative process which begins with an initial estimate of the object being reconstructed and then improves on this initial estimate via a sequence of estimates that presumable converge to ‘optimum’ reconstruction after some large number of iterations[25]. The primary advantage of this implicit linear method over other solutions to a linear system, such as the truncated singular-value-decomposition (SVD), and conjugate-gradient (CG) algorithm, is its near independence of memory occupation since only one row in the weight matrix is needed at one time. In our study, the number of ART iterations is 20 and a relaxation parameter is 0.5. The relaxation parameter has a significant influence on the results and has been proved in a range of [0, 2] [26]. The initial point is set to be homogeneous to the background.

Generally speaking, different Laplace transform-factor values emphasize different parts of the time-resolved curve. Thus analyses with different values of transform-factor β possess different capabilities for extracting spatial information. The choice of β value is a trade-off between resolution (better at earlier times) and the signal-to-noise ratio (SNR) (better at later times) [23, 27]. The Laplace transformed time-dependent diffusion model with a real transform-factor β could be mathematically interpreted as the time-independent coupled diffusion equations with absorption coefficients of (μ_aν(r)c+β≥0(ν∈{x,m}). In addition, we also know that S_m(r,r _s,β)=cΦ_x(r,r _s,β)ημ_af(r)/[1+βτ(r)], so there should be [1+βτ(r)]≥0. Considered the two conditions together, the choice of the transform-factor β should meet the condition that β≥max{-μ_aνc, -1/τ}, so that the time-independent diffusion model is physically meaningful as well mathematically stable. To investigate the effect of transform-factors on reconstructed images, a variety of transform-factors are empirically employed.

In biologic tissue, one kind of biochemical molecule of interest is targeted by the specific probe and a set of two unknown distribution i.e. the fluorescent yield and lifetime of the probe are to be recovered. The reconstructed quantity x(r,β) is a function of the fluorescence yield, lifetime and transform-factor β, so employing a pair of factors β ₁, β ₂ and the corresponding reconstructed x(r,β ₁) and x(r,β ₂), we can get the profiles of yield and lifetime simultaneously.

3. Validations

In this section, we present numerical results to demonstrate the effectiveness of the proposed method. The algorithm is validated using simulated data generated from the forward diffusion models with the corresponding optical and fluorescent properties. The phantoms shown in Fig. 1 are cylinders with a radius of 15mm and a height of 40mm, containing background and fluorescent targets that are spheres with a radius of 3mm. 64 coaxial source-detector optodes are assumed at four different layers (Z=6mm, Z=16mm, Z=26mm and Z=36mm). On every layer 16 source-detector pairs are distributed at equal spacing. When 64 sources illuminate the surface in serial, 64 detectors will collect the photons in parallel. This leads to a total of 64×64 time-resolved measurements. A finite element method with 52377 nodes and 97200 elements (prisms) is used to solve the diffusion equations. The process is implemented in MATLAB.

In our work, the optical properties of the background and targets are both set to μ_ax,m=0.035mm^-1 and μ′_sx,m=1.0mm^-1 for the excitation and emission wavelengths, and the fluorescent properties of background are ημ_af=0.001 mm^-1 and τ=100 ps. These values are in the range of the optical properties for in vivo muscle [28]. Five cases are performed to comprehensively validate the efficiency of the algorithm. The data-types are noiseless sets expect for the fifth case.

 

Fig. 1. phantom1 (left) and phantom2 (right)

Download Full Size | PDF

 

Fig. 2. Original (I) and reconstructed (II) images of different fluorescent yield and lifetime for target 1 and target 2 at Z=22mm.

Download Full Size | PDF

Case1: different yield and lifetime

Firstly, we perform the numerical experiments with two fluorescent targets of different yield and lifetime at the same positions along Z-axis. Target 1 and 2 are symmetric about the X-axis with their centers located at (X=-5mm, Y=0mm, Z=22mm) and (X=5mm, Y=0mm, Z=22mm), respectively, as shown in Fig. 1(a). Different fluorescence properties of (ημ_af1=0.003mm^-1, τ ₁=200 ps) and (ημ_af2=0.002mm^-1, τ ₂=300 ps) are assumed for the two targets to investigate the ability of the algorithm to discern the difference of the fluorescent properties within the targets. Fig. 2 shows the original and reconstructed 2-D slice images for the fluorescent yield and lifetime at Z=22mm. The reconstructed images show a reasonable fidelity in the position and size. More precisely, the reconstructed target sizes in the yield images are approximately the same as the original one while reconstructed target sizes in lifetime images are slightly larger than the true ones. However, the grayscale reconstruction is under-estimated for both the parameters.

Case2: quantitativeness

The quantitativeness of the reconstructed yield and lifetime is investigated in this case. The positions of target 1 and target 2 are the same as case 1, but the contrast of both yield and lifetime to background varies within 1.5:1(ημ^1,2 _af=0.0015mm^-1, τ^1,2=150 ps), 2:1(ημ^1,2 _af=0.002mm^-1, τ^1,2=200 ps), 3:1(ημ^1,2 _af=0.003mm^-1, τ^1,2=300 ps), 4:1(ημ^1,2 _af=0.004mm^-1, τ^1,2=400 ps) and 5:1(ημ^1,2 _af=0.005mm^-1, τ^1,2=500 ps). A measure, referred to as quantitativeness ratio, is introduced for the evaluation, which defined as the ratio of the peak value of the reconstructed yield (lifetime) in the X-profile to original value. Fig. 3(a) presents the original and reconstructed cross-sectional images of 3:1 contrast for target 1 and target 2 at Z=22mm, Fig. 3(b) demonstrates their profiles along X-axis. The quantitativeness ratio of reconstructed yield and lifetime as a function of the contrast between the target and the background are given in Fig. 3(c). It is shown that the quantitativeness ratio of the reconstructed yield consistently decreases with the increase of the target contrast, while a decrease at the lower contrast is observed for the lifetime reconstruction. For both the yield and lifetime reconstructions, the quantitativeness ratios are over 65% for all the cases. In addition, it is qualitatively seen from Fig. 3 (a) that the spatial resolution of yield is superior to that of lifetime.

 

Fig. 3. (a) Original (top) and reconstructed (bottom) images of fluorescent yield and lifetime with the contrast of 3:1 for target 1 and target 2 at Z=22mm, (b) original (blue) and reconstructed (red) profiles along X-axis for 3:1 contrast, and (c) the quantitativeness ratio of reconstructed yield (dot line) and lifetime (square line) as a function of the target contrast.

Download Full Size | PDF

Case3: the spatial- resolution

Our reconstruction algorithm is also performed with a variable center-to-center separation (CCS) of targets to evaluate its performance in spatial resolution. For simplicity, the fluorescent yield and lifetime of Target 1 and 2 are both set to ημ_af=0.002mm^-1 and τ=300 ps. Fig.4 (a) illustrates the reconstructed images for the CCS of 10mm and 8mm. It can be seen from Fig.4 (b) that the valley value between two targets in the reconstructed X-profile υ(x) increases with the CCS deceasing from 12mm to 6mm, exhibiting an increasing degradation in the spatial-resolution. To quantify the spatial-resolution performance, we define R_υ=(υ _max-υ(0))/(υ _max-υ _min), υ∈(ημ_af,τ), where υ _max=max[υ(x)], υ _min=min[υ(x)], as the measure. With the definition, R_υ=1 and R_υ=0 represent the maximal and minimum spatial-resolution, respectively. For the cases here, R_υ at the CCS of 12mm, 10mm, 8mm and 6mm are calculated as 0.894, 0.590, 0.298, 0.132 and 0.0757, 0.390, 0.123, 0.0734 for the yield and lifetime, respectively. The results show that even the edge-to-edge separation of the two targets decreases to zero, i.e. CCS=6mm, a feeble separation is still observable. Overall, the spatial-resolution of the reconstructed yield is better than that of the reconstructed lifetime.

 

Fig. 4. (a) The reconstructed yield and lifetime images for CCS=10mm (I) and 8mm (II) at Z=2mm, and (b) profiles of reconstructed images along the X-axis for a variety CCS of 12mm (blue), 10mm (red), 8mm (green) and 6mm (black).

Download Full Size | PDF

Case 4: different positions in Z-axis

 

Fig. 5. (a) The reconstructed images of fluorescent yield and lifetime on the different cross sections at Z=24mm (I), Z=23mm (II) and Z=22mm (III) respectively, and (b) their images for longitudinal section at Y=0mm(center).

Download Full Size | PDF

In this case, the ability of the algorithm to reconstruct different positions in Z-axis is investigated for Target 2 and 3. The phantom is shown in Fig.1 (b). The center of target 3 is placed aced at Z=25mm, while the X and Y coordinates are the same as those of Target 1 in Fig.1 (a). The reconstructed cross-sectional images of the yield and lifetime at Z=24mm, Z=23mm and Z=22mm, respectively, are shown in Fig.5 (a). The figure shows that with Z increasing, the radius of target 2 increases but that of target 3 decreases. The variations in the target sizes correctly reflect the real situation. The Fig.5 (b) shows the longitudinal sectional images of the yield and lifetime at Y=0mm. To test the ability of algorithm to reconstruct target’s locations, we restrict two regions of interest with a radius of 5mm and the same centers as Target 2 and 3, and then define $X_c=\frac{\left(\sum _i{\upsilon }_{\left(i\right)}{x}_{\left(i\right)}\right)}{\sum _i{\upsilon }_{\left(i\right)}}$ , $Y_c=\frac{\left(\sum _i{\upsilon }_{\left(i\right)}{y}_{\left(i\right)}\right)}{\sum _i{\upsilon }_{\left(i\right)}}$ , $Z_c=\frac{\left(\sum _i{\upsilon }_{\left(i\right)}{z}_{\left(i\right)}\right)}{\sum _i{\upsilon }_{\left(i\right)}}$ and υ _(i)∈(μη _af(i),τ _(i)) to calculate the reconstructed target centers for quantifying the positional change. Table 1 gives the original and reconstructed “centers of gravity” of target 2 and target 3. The results show a very faithful localization performance of the algorithm.

Table 1. Original and reconstructed centers of Target 2 and 3

View Table

Case 5 Noise-robustness

 

Fig. 6. (a) The reconstructed images of fluorescent yield and lifetime at Z=22mm with the SNR of 35dB (top), 25dB (bottom) and a pair of transform-factors ∓0.5β_L, and (b) the reconstructed images at Z=22mm with SNR of 35dB (top), 25dB (bottom) and a pair of transform-factors ∓β_L

Download Full Size | PDF

Finally, we evaluate the noise-robustness of the algorithm. The fluorescent yield and lifetime of target 1 and target 2 are the same as case 3 and the CCS=10mm. Originally, the noise is embedded in the measured time-resolved flux and is composed of many types. But here we model the noise directly in each of the featured data-types, i.e., the Laplace transforms of the time-resolved data, as an additive Gaussian random variable with a standard deviation proportional to the data-type: σ_υ(r ⁱ_d,r ^j_ss,β)=I_υ(r ⁱ_d,r ^j_s,β)(1+10^-χ/20 R_N), where χ is the SNR and R_N is the standard normally distributed random number. In our study we empirically selected a pair of factors β_L=∓1/[1/(μ ^(B) _axc)+1/(μ ^(B) _amc)+τ ^(B)], where the superscript B denotes background. A variety of transform-factors are employed, these values change with a step of ∓0.1β_L. Here we just show the reconstructed images based on two pairs of transform-factors and two kinds of SNR for the limitation in paper length. Fig.6 (a) illustrates the reconstructed cross-sectional images with SNRs of 35dB and 25dB, and a pair of transform-factors ∓0.5β_L. Fig.6 (b) shows the cross-sectional images with SNRs of 35dB and 25dB and a pair of transform-factors ∓β_L. The results reveal that the image quality with ∓β_L is improved evidently and the reconstruction of the fluorescent yield is more robust to the noise. More detailed investigation also shows that the quality of reconstructed images improves with increase in the difference between the two transform-factors. It seems that β_L=∓1/[1/(μ ^(B) _axc)+1/(μ ^(B) _amc)+τ ^(B)] are extremum, since reconstructions with a transform-factor larger than 1/[1/(μ ^(B) _axc)+1/(μ ^(B) _amc)+τ ^(B)] or smaller than -1/[1/(μ ^(B) _axc)+1/(μ ^(B) _amc)+τ ^(B)] always fails to produce the correct images.

4. Discussion and conclusions

It is clearly seen from the reconstructed images that the shape recovery is approximately correct, while the centers of the reconstructed targets appear to be darker and the reconstructed sizes are larger than the original ones, especially for lifetime. The latter adversities can be attributed to ill-posedness of the inverse issue and to the use of the row-fashioned ART inversion scheme in part. It is believed that large-size targets, adaptive mesh, more measurable data and a priori knowledge incorporation will be beneficial to high image quality [29–31]. In addition, the case 2 reveals that a high-contrast always goes against reconstruction for the yield, which might be interpreted as the result of breakdown in the first-order diffusion theory at the interface between drastically distinct optical media. It is worthy to note the reconstructed yields almost are superior to the lifetime in any aspect, but there is no reasonable explanation so far. In addition, the reconstructed yields are more obvious than the reconstructed lifetimes in the normally prevalent edge artifacts. Incorporating structural information and prior regularization map to reduce the artifacts around the source and detector boundaries has been discussed in Ref. 31.

In this work, we employ the ART due to its modest memory requirements for large inversion problems as well its performance stability. Working with the ART, a small relaxation parameter value has been proved to make the inversion more robust but slow down the convergence. In practice, however, its selection is still empirical and situation-dependent.

It is also worthy to point out that, different transform-factor pairs have no noticeable effect on image quality on the condition of ideal data (these results are not shown in this paper), but with noise data taken into consideration, their effects become obvious as shown in Fig.6. Empirically, the reconstruction with the extremum values for the transform-factor, β=∓1/[1/(μ ^(B) _axc)+1/(μ ^(B) _amc)+τ ^(B)], has shown an optimal noise-robustness.

In summary, we have developed a 3D image reconstruction scheme for time-domain FMT. The scheme employs on a finite element method solution to the Laplace transformed coupled diffusion equations and solves ART-based linear inversion with a pair of noise-robust transform-factors. The emphasis of this study is on the feasibility of the proposed algorithm to reconstruct the fluorescent yield and lifetime simultaneously through simulated data. We have demonstrated the ability of the algorithm to reconstruct fluorescent yield and lifetime, evaluated the spatial-resolution and quantitativeness performances, and investigated the feasibility of improving the algorithm noise-robustness by selecting appropriate transform-factors, which is crucial to the FMT applications. The results show that the method accurately retrieves the position and shape of the target. For the contrast of 1.5:1, the quantitativeness ratios of reconstructed yield and lifetime are up to 94.3% and 80% respectively, and are over 65% for both the parameters at other contrasts. The spatial-resolution achieved is promisingly high. With the proposed noise-robust factors, the algorithm can endure a SNR down to 25dB.

The authors acknowledge the funding supports from the National Natural Science Foundation of China (60578008, 60678049), the National Basic Research Program of China (2006CB705700) and Tianjin Municipal Government of China (07JCYBJC06600).