## Abstract

As a visualizing and quantitative method, Fluorescence Molecular Tomography (FMT) has many potential applications in biomedical field and its three-dimensional (3D) implementation is needed in both theory and practice. In this paper, we propose a 3D scheme for time-domain FMT within the normalized Born-ratio formulation. A finite element method solution to the Laplace transformed time-domain coupled diffusion equations is employed as the forward model, and the resultant linear inversions at two distinct transform-factors are solved with an algebraic reconstruction technique to separate fluorescent yield and lifetime images. The algorithm is validated using simulated data for 3D cylinder phantoms, and the spatial resolution and quantitativeness of the reconstruction assessed. We demonstrate that the proposed approach can accurately retrieve the positions and shapes of the targets with high spatial resolution and quantitative accuracy, and tolerate a signal-to-noise ratio down to 25dB by appropriately choosing the transform factors.

© 2008 Optical Society of America

## 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}

_{0}

*c*Φ

_{x}(

**r**,

*t*′)Γ(

**r**)e

^{[-(t-t′)Γ(r)]}d

*t*′(Γ(

**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*)/[2(1+

_{f}*R*)], where

_{f}*R*≈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].

_{f}#### 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**

*, and an integral equation is obtained according to Eq. (2)*

_{s}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**

*, and complies with the reciprocity theorem.*

_{d}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

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\approx \frac{1}{{I}_{x}({\mathbf{r}}_{d},{\mathbf{r}}_{s},\beta )}\sum _{e=1}^{E}c{\overline{G}}_{m}({\mathbf{r}}_{d},\mathbf{r},\beta ){\overline{\Phi}}_{x}(\mathbf{r},{\mathbf{r}}_{s},\beta )\overline{x}(\mathbf{r},\beta ){V}_{e}$$

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*. Eq. (4) is a set of linear equations that can be expressed in matrix form

_{e}where *I ^{nb}*(

*β*)=[

*I*(

^{nb}**r**

^{1}

*,*

_{s}**r**

^{1}

*,*

_{d}*β*),

*I*(

^{nb}**r**

^{2}

*,*

^{s}**r**

^{1}

*,*

_{d}*β*),…,

*I*(

^{nb}**r**

*,*

^{S}_{s}**r**

*,*

^{D}_{d}*β*)]

*,*

^{T}**r**

*and*

^{i}_{s}**r**

*denote the*

^{j}_{d}*i*th source and the

*j*th detector respectively. The unknown quantity

*x*(

**r**,

*β*)≈

**x**

*(*

^{T}*β*)

**u**(

**r**), where

**x**(

*β*)=[

*x*

_{1}(

*β*),

*x*

_{2}(

*β*), …

*x*(

_{N}*β*)]

*and*

^{T}**u**(

**r**)=[

*u*

_{1}(

**r**),

*u*

_{2}(

**r**),…,

*u*(

_{N}**r**)]

*is the shape function.*

^{T}**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**

*,*

^{i}_{d}**r**,

*β*) and Φ

_{x}(

**r**,

**r**

*,*

^{i}_{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{-

*μ*, -1/

_{aν}c*τ*}, 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 *β*
_{1}, *β*
_{2} and the corresponding reconstructed *x*(**r**,*β*
_{1}) and *x*(**r**,*β*
_{2}), 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.

#### 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}, *τ*
_{1}=200 ps) and (ημ_{af2}=0.002mm^{-1}, *τ*
_{2}=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.

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

#### Case 4: different positions in Z-axis

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.

#### Case 5 Noise-robustness

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

*,*

^{i}_{d}**r**

*,*

^{j}_{ss}*β*)=

*I*(

_{υ}**r**

*,*

^{i}_{d}**r**

*,*

^{j}_{s}*β*)(1+10

^{-χ/20}

*R*), where

_{N}*χ*is the SNR and

*R*is the standard normally distributed random number. In our study we empirically selected a pair of factors

_{N}*β*=∓1/[1/(

_{L}*μ*

^{(B)}

*)+1/(*

_{ax}c*μ*

^{(B)}

*)+*

_{am}c*τ*

^{(B)}], where the superscript

*B*denotes background. A variety of transform-factors are employed, these values change with a step of ∓0.1

*β*. 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/(

_{L}*μ*

^{(B)}

*)+1/(*

_{ax}c*μ*

^{(B)}

*)+*

_{am}c*τ*

^{(B)}] are extremum, since reconstructions with a transform-factor larger than 1/[1/(

*μ*

^{(B)}

*)+1/(*

_{ax}c*μ*

^{(B)}

*)+*

_{am}c*τ*

^{(B)}] or smaller than -1/[1/(

*μ*

^{(B)}

*)+1/(*

_{ax}c*μ*

^{(B)}

*)+*

_{am}c*τ*

^{(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)}
* _{ax}c*)+1/(

*μ*

^{(B)}

*)+*

_{am}c*τ*

^{(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).

## References

**1. **R. Weissleder, C.H. Tung, U. Mahmood, and A. Bogdanov, “In vivo imaging with protease-activated near-infrared fluorescent probes,” Nat. Biotechnol. **17**, 375–378 (1999). [CrossRef] [PubMed]

**2. **V. Ntziachristos, C-H Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. **8**, 757–60 (2002). [CrossRef] [PubMed]

**3. **Achilefu, R. Dorshow, J. Bugaj, and R. Rajapopalan, “Novel receptor-targeted fluorescent contrast agents for in vivo tumor imaging,” Invest. Radiol. **35**, 479–485 (2000). [CrossRef] [PubMed]

**4. **T.F. Massoud and S.S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. **17**, 545–580 (2003). [CrossRef] [PubMed]

**5. **V. Ntziachristos, C. Bremer, E.E. Graves, J. Ripoll, and R. Weissleder, “In vivo tomographic imaging of near-infrared fluorescent probes,” Molecular Imaging. **1**, 82–88 (2002). [CrossRef]

**6. **S.R. Cherry, “In vivo molecular and genomic imaging: new challenges for imaging physics,” Phys. Med. Biol. **49**, R13–48 (2004). [CrossRef] [PubMed]

**7. **H. Jiang, “Frequency-domain fluorescent diffusion tomography: a finite-element-based algorithm and simulation,” Appl. Opt. **37**, 5337–5343 (1998). [CrossRef]

**8. **A.B. Milstein, S. Oh, K.J. Webb, C.A. Bouman, Q Zhang, D.A. Boas, and R.P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. **42**, 3081–94 (2003). [CrossRef] [PubMed]

**9. **A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3D fluorescence tomography,” Opt. Express **13**, 9847–9857 (2005). [CrossRef] [PubMed]

**10. **M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three-dimensional, Bayesian image reconstruction from sparse and noisy data sets: Near-infrared fluorescence tomography,” Pans. **99**, -9624 (2002) [CrossRef]

**11. **A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein, “Three-dimensional fluorescence lifetime tomography,” Med. Phys. **32**, 992–1000 (2005) [CrossRef] [PubMed]

**12. **B. Yuan and Q. Zhu, “Separately reconstructing the structural and functional parameters of a fluorescent inclusion embedded in a turbid medium,” Opt. Express **14**, 7172–7187 (2006). [CrossRef]

**13. **A. B Milstein, J.J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffusion tomography using multiple-frequency data”, Opt. Soc. Am. A **21**,1035–1049 (2004) [CrossRef]

**14. **A. Soubret and Vasillis Ntziachristos, “Fluorescence molecular tomography in the presence of background fluorescence,” Phys. Med. Biol. **51**, 3983–4001 (2006). [CrossRef] [PubMed]

**15. **A.T.N. Kumar, J. Skoch, B.J. Bacskai, D.A. Boas, and A.K. Dunn, “Fluorescent-lifetime-based tomography for turbid media,” Opt. Lett. **30**, 3347–3349 (2005). [CrossRef]

**16. **F. Gao, H.J. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomogrphy,” Opt. Express **14**, 7109–7124(2006). [CrossRef] [PubMed]

**17. **V. Y. Soloviev, K. B. Tahir, J. McGinty, D. S. Elson, M. A. A. Neil, A. Sardini, P. M. W. French, and S. R. Arridge “Fluorescence lifetime tomography by using time-gated data acquisition,” Appl. Opt. **46**, 7384–7391(2007). [CrossRef] [PubMed]

**18. **S. Lam, F. Lesage, and X. Intes, “Time domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express **13**, 2263–2275 (2005). [CrossRef] [PubMed]

**19. **F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. **41**, 778–791 (2002). [CrossRef] [PubMed]

**20. **H.J. Zhao. Zhao, F. Gao, Y. Tanikawa, K. Homma, and Y. Yamada, “Time-resolved optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,” Appl. Opt. **43**, 1905–1916 (2005). [CrossRef]

**21. **E.M.C. Hillman, J.C. Hebden, M. Schweiger, H. Dehghani, F.E.W. Schmidt, D.T. Delpy, and S.R. Arridge, “Time resolved optical tomography of the human forearm,” Phys. Med. Biol. **46**, 1117–1130 (2002). [CrossRef]

**22. **F. Gao, Y. Tanikawa, H.J. Zhao, and Y. Yamada, “Semi-three-dimensional algorithm for time-resolved diffuse optical tomography by use of the generalized pulse spectrum technique,” Appl. Opt. **41**, 7346–7358(2002). [CrossRef] [PubMed]

**23. **V. Y. Soloviev, J. McGinty, K. B. Tahir, M. A. A. Neil, A. Sardini, J. V. Hajnal, S. R. Arridge, and P. M. W. French, “Fluorescence lifetime tomography of live cells expressing enhanced green fluorescent protein embedded in a scattering medium exhibiting background autofluorescence,” Opt. Lett. **32**, 2034–2036 (2007). [CrossRef] [PubMed]

**24. **J. Wu, “Convolution picture of the boundary conditions in photon migration and its implications in time-resolved optical imaging of biological tissues,” J. Opt. Soc. Am. A **14**, 280–287(1997). [CrossRef]

**25. **A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging **24**, 1377–1386 (2005). [CrossRef] [PubMed]

**26. **A.C Kak and M. Slaney, *Principle of Computerized Tomography Imaging*, (IEEE Press, New York, 1988).

**27. **J Wu, L. Perelman, R. R. Dasari, and M. S. Feld, “Fluorescence tomography imaging in turbid media using early-arriving photons and Laplace transforms,” Proc. Natl. Acad. Sci. USA. **94**, 8783–8788 (1997). [CrossRef] [PubMed]

**28. **F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg, and C. Depeursinge, “IN vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. **38**, 4939–4950(1999). [CrossRef]

**29. **J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, “Full adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue,” Opt. Express **15**, 6955–6975 (2007). [CrossRef] [PubMed]

**30. **D. Wang, X. Song, and J. Bai, “Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution,” Opt. Express **15**, 9722–9730 (2007). [CrossRef] [PubMed]

**31. **S. C. Davis, H. Dehghani, J. Wang, S. D. Jiang, B. W. Pogue, and K. D. Paulsen, “Image-guided diffuse optical fluorescence tomography implemented with Laplacian-type regularization,” Opt. Express **15**, 4066–4082(2007) [CrossRef] [PubMed]