## Abstract

In this paper, we propose a dual-excitation-mode methodology for three-dimensional (3D) fluorescence molecular tomography (FMT). For this modality, an effective reconstruction algorithm is developed to reconstruct fluorescent yield and lifetime using finite element techniques. In the steady state mode, a direct linear relationship is established between measured optical data on the body surface of a small animal and the unknown fluorescent yield inside the animal, and the reconstruction of fluorescent yield is formulated as a linear least square minimization problem. In the frequency domain mode, based on localization results of the fluorescent probe obtained using the first mode, the reconstruction of fluorescent lifetime is transformed into a relatively simple optimization problem. This algorithm helps overcome the ill-posedness with FMT. The effectiveness of the proposed method is numerically demonstrated using a heterogeneous mouse chest phantom, showing good accuracy, stability, noise characteristics and computational efficiency.

©2005 Optical Society of America

## 1. Introduction

Near infrared light (NIR) can travel a few centimeters in the tissue due to the low photon absorption in spectral window of 700–900nm, and can be captured on the tissue surface. Hence, the fluorescence enhanced NIR optical tomography was developed to probe the biological processes at anatomical, functional, cellular and molecular levels by tomographically determining the fluorescent yield and lifetime[1]. Fluorescence molecular tomography has great potential to become a promising imaging modality for small animal imaging *in vivo* [2, 3].

Rapid advancement has been made in the research of fluorescent imaging recently. Several reconstruction approaches for FMT have been proposed, including Born type approximation techniques based on the diffusion model [4, 5], reconstruction algorithms of the fluorescent yield with the steady-state equation of radiative transfer (ERT) model [6, 7], finite element based methods for reconstruction of both lifetime and fluorophore yield using frequency domain data [8, 9], an adaptive finite-element-based method [10], time resolved techniques for fluorescence diffuse optical tomography [11], and a reconstruction algorithm of lifetime and yield using a two stage Bayesian approximate extended Kalman filter [12].

In fluorescence imaging, the continuous wave (CW) and modulated light are two popular excitation modes. The CW mode uses light of constant intensity, and offers excellent detection characteristics, high quantum efficiency and optimum signal-to-noise performance[3]. The major disadvantage of the steady-state method is the inability to detect fluorescence lifetime information. The frequency domain technique uses light of a modulated intensity at frequency ranging from 70MHz–150MHz, and the measurement of the light intensity and phase shift of the photon wavefront reveals the fluorescence lifetime information [12].

In this paper, we combine the advantages of both modes, and propose a dual-excitation-mode method to reconstruct fluorophore yield and lifetime in the steady state and frequency domain, respectively. We excite the fluorophore targets and measure the photon density on the boundary of an imaged object in both continuous and modulated light modes. By utilizing the dual-excitation-mode, we develop a finite-element, diffusion model based method for FMT. A direct linear relation between fluorescent yield and the photon density measured on the boundary in the steady state is established by the finite element discretization. The problem is then converted to that of a standard linear least square optimization. We incorporate a *priori* knowledge of feasible fluorophore concentration regions as constraints in our reconstruction algorithm, which effectively reduce the ill-posedness of the problem [13, 14, 15]. We then reconstruct the lifetime information from the frequency domain measurement by the finite element method. The reconstructed fluorophore’s spatial distribution provides a strict constraint that makes the optimization process better posed, ensuring accurate and stable reconstruction of the lifetime.

In the next section, we present our finite-element-based diffusion equation method. In the third section, we demonstrate our method by reconstructing the fluorescent yield and lifetime in numerical simulation based on a heterogeneous phantom. In the last section we discuss the relevant issues and present the conclusions.

## 2. Methodology

#### 2.1. Diffusion models

In a dual-excitation-mode fluorescence imaging experiment, a continuous wave (CW) external excitation light source at wavelength *λ*_{x}
was first applied at a point on the surface of tissue, followed by the modulated excitation light of the same wavelength with modulation frequency *ω*. The excitation photons penetrate the tissue, where the scattering predominates over the absorption, and the excitation photon propagation can be quite accurately described by the diffusion approximation to the radiative transport equation [16, 17]

where **r** is the position vector, *Φ*
_{x} (**r**, *ω*) the excitation photon density, Ω the imaged region, *c* the velocity of light. *D*_{x}
(**r**) and *μ _{ax}* (

**r**) are the diffusion and absorption coefficients at excitation wavelength

*λ*

_{x}. The diffusion coefficient is given by

*D*

_{x}(

**r**) = (3(

*μ*(

_{ax}**r**) + (1 -

*g*)

*μ*(

_{sx}**r**)))

^{-1}, where

*μ*(

_{sx}**r**) is the scattering coefficient, and

*g*the anisotrophy parameter. The diffusion equation is complemented by the Robin-type boundary conditions on the boundary

*∂*Ω [18, 19]

where **n** denotes the outward normal of the tissue surface, *S*(**r**, *ω*) the excitation light source on the boundary, and *ρ* a constant depending upon the optical reflective index mismatch at the boundary [19].

The fluorophore probe inside the tissue absorbs the incident photons and decays to its ground state, emitting photons at wavelength *λ*_{m}
, where *λ*_{m}
> *λ*_{x}
. The emission photons escape from the tissue and are measured on the surface *∂*Ω. In the modulated excitation light mode, the light has a constant modulation frequency *ω* and can be mesured with a gain modulated image intensified CCD camera system [20]. In the CW excitation mode (*ω* = 0), a steady-state dataset are collected using the CCD camera. The emission photon propagation in the tissue is modeled by the following diffusion equation [1,8]

where **Φ**
_{m} (**r**, *ω*) is the emission photon density, *D*_{m}
(**r**) and *μ _{am}* (

**r**) are the diffusion and absorption coefficients of tissue at emission wavelength

*λ*

_{m}. The diffusion coefficient is given by

*D*

_{m}(

**r**) = (3(

*μ*(

_{am}**r**) + (1 -

*g*)

*μ*(

_{sm}**r**)))

^{-1}, where

*μ*(

_{sm}**r**) is the scattering coefficient at

*λ*

_{m}.

*η*(

**r**)

*μ*(

_{af}**r**) is the fluorescent yield.

*η*is the fluorophore’s quantum efficiency and

*μ*(

_{af}**r**) is the fluorescent absorption coefficient, which is directly proportional to the fluorophore concentration. τ (

**r**) is the lifetime, which relates to the local chemical enviroment. The emission light is subjected to the Robin-type boundary condition

#### 2.2. Reconstruction method

To reconstruct the fluorescent yield and lifetime, the region Ω is first discretized into *N*_{e}
elements (Ω_{1},Ω_{2}, ⋯, Ω_{Ne}) and *N*_{v}
vertex nodes (*V*
_{1},*V*
_{2}, ⋯,*V _{Nv}*), such that Ω = ∪

^{Ne}

_{k=1}Ω

_{k}. Based on the discretization of the region, the photon density

**Φ**

_{x}(

**r**,

*ω*) is approximated with piecewise polynomial interpolation functions[16]

Using the finite element method, the excitation light diffusion equation Eq. (1) is formulated into the following matrix equation

The components *a*_{x}*, c*_{x}*, p*_{x}*, q*_{x}
and *h*_{x}
corresponding to matrices **A**
_{x}, **C**
_{x}, **P**
_{x}, **Q**
_{x} and **H**
_{x} are given as follows

**Φ**
_{m} (**r**, *ω*) and *U*(**r**, *ω*, *η*, *μ _{af}*,

*τ*) are approximated with piecewise polynomial basis functions [19]

The emission light diffusion equation Eq. (3) is also transformed into matrix equation

The components *a*_{m}*, c*_{m}*, p*_{m}*, q*_{m}
and *f*_{m}
corresponding to matrices **A**
_{m}, **C**
_{m}, **P**
_{m}, **Q**
_{m} and **F**
_{m} are given as follows

### 2.2.1. Reconstruction of fluorescent yield

We take *ω* = 0 in Eq. (6) to describe the excitation photon propagation in the steady state. The excitation photon density can be solved from Eq. (6)

where **K**
_{x} = (**A**
_{x} + **C**
_{x} + **Q**
_{x}). Set *ω* = 0 in Eq. (9), the emission photon density in steady state is formulated as

where **K**
_{m} = (**A**
_{m} + **C**
_{m} + **Q**
_{m}) and *U*(*ημ _{af}*) =

*ημ*.

_{af}The reconstruction of the fluorescent yield is to recover the vector *U* from measured photon density *Φ*
_{m} on *∂*Ω in the CW mode, so **Φ**
_{m} can be partitioned into measurable boundary values ${\mathbf{\Phi}}_{m}^{\mathit{\text{meas}}}$
and unmeasurable interior values ${\Phi}_{m}^{\mathit{\text{igno}}}$
. Since the reconstruction of fluorescent yield is underdetermined and ill-posed, we incorporate the *a priori* knowledge obtained from measured data as well as physiological and anatomical information. From the *a priori* knowledge, the vector *U* is divided into two parts: *U*_{p}
and *U*
_{0}. *U*_{p}
corresponds to the permissible region Ω_{p}, in which the fluorescent probe is located. *U*
_{0} corresponds to the forbidden region Ω_{0} with no fluorophore concentration, that is *U*
_{0} = 0. Let **B**
^{y} = ${\mathbf{K}}_{m}^{-1}$
**F**
_{m}, we remove the columns of matrix **B**
^{y} corresponding to the vector *U*
_{0}, and remove the rows of **B**
^{y} corresponding to ${\Phi}_{m}^{\mathit{\text{igno}}}$
to obtain ${\mathbf{B}}_{p}^{y}$
. Therefore, we obtain a linear relationship between ${\Phi}_{m}^{\mathit{\text{meas}}}$
and *ημ _{af}*

The measured data in fluorescent imaging are usually noisy, hence it is impractical to solve for ${U}_{p}^{y}$ from the linear system Eq. (13) directly. An optimization approach is employed to find a regularized solution by minimizing the following objective function

where *α* is the upper bound that is chosen to be physically meaningful, and *u*_{k}
represents the component values ${U}_{p}^{y}$
. Λ is the weight matrix and ∥**V**∥_{Λ} = **V**
^{T}
**ΛV**, *β* is a stabilizing function, and σ the regularization parameter. This minimization is a standard linear least square problem.

### 2.2.2. Reconstruction of fluorescent lifetime

We set *ω* to the modulation frequency of the excitation light in Eq. (6) for expressing the excitation photon propagation through tissue in the frequency domain. The photon density can be obtained from

where **K**
_{x} = (**A**
_{x} + **C**
_{x} + *iω*
**P**
_{x} + **Q**
_{x}). In the frequency domain the emission photon density of Ω is formulated as

where **K**
_{m} = (**A**
_{m} + **C**
_{m} + *iω*
**P**
_{m} + **Q**
_{m}) and $U\left(\omega ,\tau \right)={\mathit{\eta \mu}}_{{a}_{f}}\frac{1+\mathit{i\omega \tau}}{1+{\omega}^{2}{\tau}^{2}}.$

Reconstructed fluorescent yields determines the fluorescent target positions, which allows the reconstruction of lifetime to be limited to specific locations, so that the number of unknown lifetime variables is minimized. To determine actual locations of fluorescent yields, the reconstructed fluorescent yields should be filtered using a threshold to eliminate undesirable components due to noise. Similar to what we did in the fluorescent yield reconstruction, by taking advantage of the reconstructed fluorescent yield distribution, we reduce **B**
^{l} to ${\mathbf{B}}_{p}^{l}$
where **B**
^{l} = ${\mathbf{K}}_{m}^{-1}$
**F**
_{m}. Therefore, we have a nonlinear equation about lifetime *τ*

where ${U}_{p}^{l}$
(*ω, τ*) is a vector function of lifetime *τ*, supported only on the reconstructed fluorescent target positions. We perform an optimization procedure to reconstruct the lifetime *τ* as follows

## 3. Simulations and results

In this section, we presents numerical results to demonstrate the effectiveness of our reconstruction algorithm, which is implemented in MATLAB.

#### 3.1. Numerical experiments setup

We perform the numerical experiments based on the heterogeneous mouse chest geometrical phantom, as in Fig. 1. The phantom is a cylinder, 30mm in diameter and 26mm in height, containing regions that resemble lung, heart, bone and muscle. The optical parameters of each region for both excitation and emission wavelength are listed in Table. 1 [9, 14]. The numerical phantom is discretized into 6576 nodes and 11340 finite elements, which are prisms. The emission photon density, which is generated by the finite-element forward model, is collected on the surface of the phantom. The measurement data consists of 1024 data points corresponding to the nodes on the phantom’s side surface.We took two measurements in both the steady state and frequency domain with two distinct positions of the excitation light sources on the phantom’s side surface. The two excitation light sources positions are shown in Fig. 2(a). The power of excitation light sources are 30mW. For frequency domain measurement, the excitation light with a modulation frequency of 100MHz is used. The measured data for both modes are corrupted with 10% random guassian noise.

#### 3.2. Single fluorescent target

First, we performed the numerical experiment with a single fluorescent target. The target is embedded in the right lung region, as shown in Fig. 3(a). The fluorophore’s quantum yield *η* was set to 0.16, the fluorophore absorption coefficient *μ _{af}* 0.05

*mm*

^{-1}and the lifetime

*τ*0.8

*ns*. We assumed no fluorophore concentration in the background region. The feasible region was selected to be the middle section of the right lung area, which consists of 300 finite elements, as shown in Fig. 2(b). In fact, it is typical to have bright and dark regions in views of surface measurment, and the regions can be constructed around the bright spots. In the numerical implementation, we set the stabilizing function

*β*(

*X*) =

*X*

^{T}

*X*, and the regularization parameter σ = 3.0 × 10

^{-9}. Then, the proposed method described in Section 2 was applied to reconstruct the fluorescent yield and lifetime of the target. The 3D localization of reconstructed target is shown in Fig. 3(b), which matched the actual target position in Fig. 3(a). The reconstructed fluorescent yield and lifetime, are shown in Fig. 4 and Fig. 5. The reconstruction results and relative errors are listed in Table. 2.

#### 3.3. Two fluorescent targets

We also performed the numerical experiment with two fluorescent targets to further demonstrate the performance of the algorithm. Both targets were embedded in the left lung region, as shown in Fig. 6(a). The two targets had separation distance of 6.8mm. The fluorophore’s quantum yield *η* was set to 0.16, the fluorophore absorption coefficients *μ _{af}* were 0.05

*mm*

^{-1}and 0.03

*mm*

^{-1}and the lifetimes

*τ*0.8

*ns*and 1.2

*ns*in target one and two, respectively. The feasible region remained the same as in the single target experiment, as in Fig. 2.

The reconstructed results are shown in Fig. 6(b), Fig. 7 and Fig. 8. These results further show that the maximum values of the reconstructed fluorescent yield and lifetime were indeed at the centers of the actual fluoresent yield and lifetime distributions, and some small reconstructed values, which are caused by noise, were scattered around the actual targets. The reconstructed targets well matched the actual counterparts. The reconstructed values of the fluorescent yield and lifetime, as well as their relative errors, are summarized in Table 3. The relative error of fluorescent yield (lifetime) for a target is herein defined as the ratio between the absolute difference between the reconstructed fluorescent yield (lifetime) and the true fluorescent yield (lifetime) over the true fluorescent yield (lifetime). The computed surface photon density distributions based on the reconstructed fluorescent yield and lifetime were in good agreement with the measured counterparts, as shown in Fig. 9. The reconstruction computation took 256 seconds CPU time on a 3.2GHz Pentium 4 PC with 1Gb memory.

## 4. Discussion and conclusion

The use of the dual modes for FMT has been shown synergistic. In the steady state mode, a linear relationship is established between measured surface optical data and the unknown fluorescent yield using the finite element analysis. The reconstruction of fluorescent yield becomes a linear inversion problem. By incorporating *a priori* knowledge of feasible fluorochrome concentration regions, the ill-posedness of the FMT problem has been significantly reduced. In the frequency domain mode, based on localization results of the fluorescent probe in the first stage, we have reduced the number of unknowns on fluorescent lifetime, and put the reconstruction of fluorescent lifetime in a simple minimization framework, which helps further overcome the ill-posedness for FMT reconstruction. Our numerical results have indicated that the proposed method has good numerical accuracy, stability, robustness, and efficiency.

Because FMT is to recover a 3D source distribution from 2D surface-based optical data, it is a typical underdetermine and ill-posed problem. To regularize the reconstruction we have imposed a permissible region constraint to reduce the number of variables. This type of *a priori* knowledge can be derived from luminescent views of the object, since it is common to have bright and dark clusters in these views, and permissible regions can be constructed around the bright spots. According to the error between the measured and computed data on the surface of the object, the permissible regions can be iteratively refined.

The emphasis of this study is on the feasibility of the proposed reconstruction method for FMT. It is recognized that there are major differences between the real small animal and the numerical model we have used as well as between the radiative transport equation and its diffusion approximation. These important factors surely affect the performance of the reconstruction method. Although we have included strong Gaussian noise in our numerical simulation to show the robustness of our algorithm, it is underlined that there may be larger errors in the reconstruction results from an *in vivo* small animal study. Although the geometrical model of our numerical phantom herein is rough relative to structures of a true small animal, there are little problems in theory and technology to handle a rather complex heterogeneous geometrical model like a mouse using finite element techniques. In practical applications, multiple excitation points and more sampling locations can be employed to obtain more information to enhance the stability of the FMT reconstruction. Also, it is noted the autofluorescence of the biological tissue leads to emission of lower-energy light when is is excited by higher-energy light. This signal may compromise the contrast in data and images. Hence, the measurement data should be filtered to suppress the autofluorescence effects [21].

In conclusion, we have developed a dual mode strategy for FMT. In doing so, the complex nonlinear inverse problem of FMT has been simplified into a linear inverse problem for fluorescent yield and an optimization problem for lifetime. Our work has provided a promising method for FMT, which is important for small animal imaging. Currently, we are planning to test our FMT method using physical phantoms and living small animals. Relevant results will be presented later.

## Acknowledgments

This work is supported by a NIH/NIBIB grant EB001685. The authors are grateful for the editorial help from Deepak K. Bharkhada.

## References and links

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

**
2
. **
G.
Choy
,
P.
Choyke
, and
S. K.
Libutti
, “
Current advances in molecular imaging: noninvasive in vivo bioluminescent and fluorescent optical Imaging in Cancer Research
,”
Mol. Imag.
**
2
**
,
303
–
312
(
2003
). [CrossRef]

**
3
. **
V.
Ntziachristos
,
J.
Ripoll
,
L. V.
Wang
, and
R.
Weissleder
, “
Looking and listening to light: the evolution of whole-body photonic imaging
,”
Nature Biotech.
**
23
**
,
313
–
320
(
2005
). [CrossRef]

**
4
. **
V.
Ntziachristos
and
R.
Weissleder
, “
Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation
,”
Opt. Lett.
**
26
**
,
893
–
895
(
2001
). [CrossRef]

**
5
. **
M. A.
O’Leary
,
D. A.
Boas
,
X. D.
Li
,
B.
Chance
, and
A. G.
Yodh
, “
Fluorescence lifetime imaging in turbid media
,”
Opt. Lett.
**
21
**
,
158
–
160
(
1996
). [CrossRef] [PubMed]

**
6
. **
A. D.
Klose
and
A. H.
Hielscher
, “
Fluorescence tomography with simulated data based on the equation of radiative transfer
,”
Opt. Lett.
**
28
**
,
1019
–
1021
(
2003
). [CrossRef] [PubMed]

**
7
. **
A. D.
Klose
,
V.
Ntziachristos
, and
A. H.
Hielscher
, “
The inverse source problem based on the radiative transfer equation in optical molecular imaging
,”
J. Comput. Phys.
**
202
**
,
323
–
345
(
2004
). [CrossRef]

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

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

**
10
. **
A.
Joshi
,
W.
Bangerth
, and
E. M.
Sevick-Muraca
, “
Adaptive finite element based tomography for fluorescence optical imaging in tissue
,”
Opt. Express
**
12
**
,
5402
–
5417
(
2004
)
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-12-22-5402
[CrossRef] [PubMed]

**
11
. **
S.
Lam
,
F.
Lesage
, and
X.
Intes
, “
Time domain fluorescent diffuse optical tomography: analytical expressions
,”
Opt. Express
**
13
**
,
2263
–
2275
(
2005
)
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-7-2263
. [CrossRef] [PubMed]

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

**
13
. **
W.
Cong
,
D.
Kumar
,
Y.
Liu
,
A.
Cong
, and
G.
Wang
, “
A practical method to determine the light source distribution in bioluminescent imaging
,”
Proc. SPIE
**
5535
**
,
679
–
686
(
2004
). [CrossRef]

**
14
. **
W.
Cong
,
G.
Wang
, and
D.
Kumar
,
*
et al.
*
, “
Practical reconstruction method for bioluminescence tomography
,”
Opt. Express
**
13
**
,
6756
–
6771
(
2005
)
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-13-18-6756
. [CrossRef] [PubMed]

**
15
. **
G.
Wang
,
Y.
Li
, and
M.
Jiang
, “
Uniqueness theorems in bioluminescence tomography
,”
Med. Phys.
**
31
**
,
2289
–
2299
(
2004
). [CrossRef] [PubMed]

**
16
. **
S. R.
Arridge
,
M.
Schweiger
,
M.
Hiraoka
, and
D. T.
Delpy
, “
A finite element approach for modeling photon transport in tissue
,”
Med. Phys.
**
20
**
,
299
–
309
(
1993
). [CrossRef] [PubMed]

**
17
. **
E.
Shives
,
Y.
Xu
, and
H.
Jiang
, “
Fluorescence lifetime tomography of turbid media based on an oxygen-sensitive dye
,”
Opt. Express
**
10
**
,
1557
–
1562
(
2002
)
http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-26-1557
. [PubMed]

**
18
. **
R. C.
Haskell
,
L. O.
Svaasand
,
T. T.
Tsay
,
T. C.
Feng
,
M. S.
McAdams
, and
B. J.
Tromberg
, “
Boundary conditions for the diffusion equation in radiative transfer
,”
J. Opt. Soc. Am. A
**
11
**
,
2727
–
2741
(
1994
). [CrossRef]

**
19
. **
M.
Schweiger
,
S. R.
Arridge
,
M.
Hiraoka
, and
D. T.
Delpy
, “
The finite element method for the propagation of light in scattering media: Boundary and source conditons
,”,
Med. Phys.
**
22
**
,
1779
–
1792
(
1995
). [CrossRef] [PubMed]

**
20
. **
A.
Godavarty
,
M. J.
Eppstein
,
C.
Zhang
,
S.
Theru
,
A. B.
Thompson
,
M.
Gurfinkel
, and
E. M.
Sevick-Muraca
, “
Fluorescence-enhanced optical imaging in large tissue volumes using a gain-modulated ICCD camera
,”
Phys. Med. Biol.
**
48
**
,
1701
–
1720
(
2003
). [CrossRef] [PubMed]

**
21
. **
J. R.
Mansfield
,
K. W.
Gossage
,
C. C.
Hoyt
, and
R. M.
Levenson
, ”
Autofluorescence removal, multiplexing, and automated analysis methods for in-vivo fluorescence imaging
,”
J. Biomed. Opt.
**
10
**
,
041207
(
2005
). [CrossRef]