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 Ne elements (Ω1,Ω2, ⋯, ΩNe) and Nv vertex nodes (V 1,V 2, ⋯,VNv), 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 ₀. U_p corresponds to the permissible region Ω_p, in which the fluorescent probe is located. U ₀ corresponds to the forbidden region Ω₀ with no fluorophore concentration, that is U ₀ = 0. Let B ^y = ${\mathbf{K}}_m^{-1}$ F _m, we remove the columns of matrix B ^y corresponding to the vector U ₀, 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.

 

Fig. 1. Mouse chest phantom. (a) finite element mouse chest geometrical model (each region is labeled by a letter, H for heart, L for lung, B for bone, and M for muscle), (b) a cross section of the model.

Download Full Size | PDF

Table 1. Optical parameters

View Table | View all tables in this article

 

Fig. 2. Imaging geometry. (a) excitation light source positions, (b) feasible region.

Download Full Size | PDF

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.05mm ^-1 and the lifetime τ 0.8ns. 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^TX, 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.

 

Fig. 3. Target location reconstruction. (a) actual target position (pointed by the arrow), (b) reconstructed target position(pointed by the arrow).

Download Full Size | PDF

 

Fig. 4. Fluorescent yield reconstruction. (a) actual fluorescent yield distribution, (b) fluorescent yield reconstructed with measurement corrupted by 10% guassian noise.

Download Full Size | PDF

Table 2. Reconstruction result of a single fluorescent target.

View Table | View all tables in this article

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.05mm ^-1 and 0.03mm ^-1 and the lifetimes τ 0.8ns and 1.2ns 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.

 

Fig. 5. Fluorescent lifetime reconstruction. (a) actual fluorescent lifetime distribution, (b) fluorescent lifetime reconstructed with measurement corrupted by 10% guassian noise.

Download Full Size | PDF

 

Fig. 6. Targets’ locations reconstruction. (a) actual target position(pointed by the arrow), (b) reconstructed target position(pointed by the arrow).

Download Full Size | PDF

Table 3. Reconstruction results of two fluorescent targets.

View Table | View all tables in this article

 

Fig. 7. Fluorescent yield reconstruction. (a) actual fluorescent yield distribution, (b) fluorescent yield reconstructed with measurement corrupted by 10% guassian noise.

Download Full Size | PDF

 

Fig. 8. Fluorescent lifetime reconstruction. (a) actual fluorescent lifetime distribution, (b) fluorescent lifetime reconstructed with measurement corrupted by 10% guassian noise.

Download Full Size | PDF

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.

 

Fig. 9. Comparison between the actual measurement data at different detection positions in the frequency domain mode and the corresponding profiles of surface photon density computed based on the reconstructed fluorescent yield and lifetime. (a) The actual measurement data and the corresponding fitting profiles with the first excitation light source, and (b) the counterparts with the second excitation light source.

Download Full Size | PDF

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.