Abstract
Fluorescence molecular tomography (FMT) is a promising technique for in vivo small animal imaging. In this paper, a two-stage reconstruction method based on normalized Born approximation is developed for FMT, which includes two steps for quantitative reconstruction. First, the localization of fluorescent fluorophore is determined by -norm regularization method. Then, in the location region of fluorophore, which is provided by the first stage, algebraic reconstruction technique (ART) is utilized for the fluorophore concentration reconstruction. The validity of the two-stage quantitative reconstruction algorithm is testified by simulation experiments on a 3D digital mouse atlas and physical experiments on a phantom. The results suggest that we are able to recover the fluorophore location and concentration.
1. Introduction
Fluorescence molecular tomography (FMT) is an emerging imaging modality that allows noninvasive, quantitative, three-dimensional (3D) imaging of fluorescence probes associated with molecular and cellular functions [1–3]. It can be applied to drug discovery and therapeutic response in small animals [4, 5]. In this imaging modality, an external excitation light illuminates the tissue surface of the small animal, which is labeled with fluorescence probes beforehand. Then, the fluorescence probes are excited to emit photons, which are detected as fluorescence on the surface. With the measured data on the surface, the spatial distribution and concentration of fluorophore can be reconstructed based on certain inverse mathematical models [6, 7].
In recent years, great efforts have been made in developing new reconstruction algorithms [7–10] and imaging system [11–13]. However, challenges remain in quantitative fluorescence tomography. In [7], Song et al. developed a novel adaptive finite element algorithm based on the longest refinement method for free-space fluorescent tomography. This algorithm recovered the spatial distribution of the fluorescence target, but the reconstructed fluorescent yield only reached a half of the true fluorescence yield. In order to reduce the influence of optical properties, the normalized Born ratio has been developed for FMT [8]. To employ a priori information of the sparseness of fluorescence target, Han et al. proposed a fast reconstruction algorithm based on sparsity regularization for FMT [9]. The fluorescence source can be localized, but the fluorescence intensity was recovered with nearly 38% error. An integrated trimodality system that combines fluorescence, diffuse optical, and X-ray tomography (FT/DOT/XCT) was developed for quantitative fluorescence tomography [10]. The diffuse optical tomography (DOT) provided the functional optical background and X-ray tomography (XCT) provided the structural information, which included the location of the fluorophore [10, 13]. The fluorophore could be localized accurately and the fluorophore concentration was recovered within 8% error only when both functional and structural a priori information was utilized for fluorescence tomography (FT) reconstruction algorithm. This paper inspired us to propose the two-stage method.
In this paper, we proposed a two-stage reconstruction algorithm based on normalized Born approximation, which can largely reduce the effects of incorrect assumption of optical properties [8]. First, we localized the fluorophore by -norm regularization method to fully take advantage of the sparsity of fluorophore. Then, in the region where the fluorophore is located, algebraic reconstruction technique (ART) is utilized to recover the fluorophore concentration. The first stage not only localizes the fluorophore but also provides the reconstruction region for the second stage. So, the first stage is very important for the whole quantitative reconstruction. The developed method was validated by both simulation experiments on a 3D digital mouse and physical experiments on a cubic phantom. The results have shown the feasibility and potential of the algorithm for quantitative fluorescence tomography.
2. Method
2.1. Forward and Inverse Problems
In the near-infrared (NIR) bandwidth (700–900 nm), the light propagation in the biological tissue can be modeled using diffusion equation (DE) with Robin-type boundary condition [14, 15]. The Green’s function due to a continuous wave (CW) point excitation source describes the light propagation. It can be obtained as follows: where is the diffusion coefficient at position . is the reconstruction domain with a boundary of . is the absorption coefficient and is a constant depending on the optical reflective index mismatch at the boundary. represents the outward normal vector to the boundary. The diffusion equation can be solved using finite element method (FEM) [15, 16]. In order to reduce the influence of tissue heterogeneity, normalized Born approximation is used [8]. The normalized ratio of the measured fluorescence signal and corresponding excitation signal at detector is given as follows [8, 17]: where is a calibration factor which depends on the genetic system characteristics. is the point inside the volume . denotes the fluorescence yield to be reconstructed, which is directly related to the fluorophore concentration by the formula , where is the molar extinction coefficient, is the quantum efficiency, and is the concentration of the fluorophore [18]. describes the light propagation at the excitation wavelength due to a point source at , and corresponds to the light propagation at the emission wavelength due to a point source at . By discretizing the volume into voxels, a linear system can be generated for multiple source-detector pairs: where is weight matrix and each element represents the contribution of th voxel of the discrete field of interest to the measurement due to the th source-detector pair. The unknown fluorescence yield in each voxel is recovered by solving (3) using certain regularization method. Detailed descriptions can be found in [17].
2.2. Two-Stage Reconstruction Algorithm
Due to the ill-posed and the ill-conditioned nature of FMT inverse problem, many kinds of regularization methods are employed to obtain a stable solution. Here, we developed a two-stage reconstruction algorithm for quantitative fluorescence tomography. Considering the fact that the volume of the fluorophore is very small compared with the whole reconstruction region, the fluorophore can be seen as sparsely distributed in the region of the interested. First, -norm regularization is utilized to fully take advantage of the sparseness of the fluorophore. As a result, (3) can be solved by the following minimization problem: where and is the regularization parameter. From the results in the first stage, we determined the permissible fluorophore region for the second stage by selecting the nodes with 30% of the maximum reconstructed value. Since the reconstruction region is small, the -norm regularization method is not suitable for the reconstruction. Then, algebraic reconstruction technique (ART) is utilized to recover the fluorescence yield. Its form is as follows [19, 20]: where is the relaxation parameter. and are the “”th and ()th estimation of the th element of the unknown variables . The initial value of the ART is provided by the first stage.
Algorithm. The flow of the proposed method is depicted in Figure 1, which mainly includes.
Stage 1. Reconstruct the distribution of fluorophore on the whole body of the mouse by minimizing (4), that is, localize the fluorophore by -norm regularization method.
Stage 2. Recover the fluorescence yield using ART on the permissible fluorophore region which is provided by the first stage, that is, quantify the fluorescence concentration.
3. Experiments and Results
In this section, two groups of verification experiments were conducted to validate the potential and feasibility of the proposed two-stage reconstruction algorithm for quantitative fluorescence tomography. Reconstructions were carried out on a personal computer with a 2.66 GHz Intel(R) Xeon(R) CPU E5430 and 8.00 GB RAM. To analyze the results quantitatively, location error (LE), which is the Euclidean distance between the centers of the reconstructed and the actual fluorophore, is defined by , where is the reconstructed coordinate with the maximum reconstructed value and is the actual coordinate. The relative fluorescence yield error is defined to be , where and are the recovered maximum fluorescence yield and the actual one, respectively.
3.1. Digital Phantom Experiment
3.1.1. Setup for Simulation Studies
In the simulation experiments, a digital mouse atlas of CT and cryosection data was utilized to provide anatomical information [21]. We only considered the torso section of the mouse with a height of 35 mm. A cylindrical fluorescent target with a 0.8 mm radius and 1.6 mm height was placed in the liver with center at 11.9 mm, 6.4 mm, and 16.4 mm as shown in Figure 2(a), and the fluorescence yield of the fluorescent target was set to be 0.05. The optical parameters for different organs are shown in Table 1 [6, 22]. The excitation and emission data were synthesized using finite element method (FEM) based on (1). In this experiment, we considered 36 point sources at different positions in sequence as shown in Figure 2(b). For each point source, 39 detectors were used for collecting fluorescence signal inside 120° field of view (FOV), as shown in Figure 2(c). So there are total 1404 source-detector pairs.
(a)
(b)
(c)
3.1.2. Reconstruction of the Single Fluorophore
The reconstruction region was discretized into 4157 nodes and 21811 tetrahedral elements. In the first stage, we employed the incomplete variables truncated conjugate gradient method (IVTCG) to solve (4) [23]. Since determining a good regularization parameter is a very challenging task, we did not represent this process here in detail. Instead, we selected the range between and which is sufficient for our experiments. In the second stage, the ART iterations were terminated when the iteration number reached the maximum ART iteration number 20 or the relative change between neighboring iteration was less than 0.5%. To clearly show the performance of the proposed two-stage reconstruction for FMT, the recovered results are showed in Figure 3. And Table 2 gives the LE and relative error of fluorescence yield at each stage. It is clear that the fluorescence yield on the selected slice for the second stage is brighter than that of the first one in Figure 3.
(a)
(b)
(c)
(d)
We also used algebraic reconstruction technique (ART) and -norm regularization method to solve this problem, respectively, as a comparison. Figure 4 shows the reconstructed results. It is clear that even though both the ART and -norm regularization methods can provide a good location of the fluorescent target, yet they both fail to recover the fluorescence yields accurately. Compared with our proposed method, they recovered the fluorescence yields with a large error.
(a)
(b)
3.2. Physical Phantom Experiment
3.2.1. Setup for Experiment Studies
A full angle noncontact FMT-imaging system is employed for collecting fluorescence signals [6]. The schematic of this system is shown in Figure 5 and detailed descriptions can be found in [6]. A 672 nm continuous wave (CW) laser source illuminated the target, which was placed at the rotational stage. A highly sensitive charge-coupled device (CCD) camera (Princeton Instruments PIXIS 2048B, Roper scientific, Trenton, NJ), which was cooled to −70°C, was utilized for recording the transmitted light at both the excitation and emission wavelengths using appropriate filters. Here, custom-made 40 nm band-pass filter (HZXD, Beijing, China) centered at 660 nm and 35 nm band-pass filter (HZXD, Beijing, China) centered at 720 nm were used to allow light transmission at the excitation and emission wavelength, respectively.
A cubic phantom with a side length of 20 mm was placed on the rotational stage, shown in Figures 6(a)–6(c). The phantom was made from polyoxymethylene. The optical parameters for both excitation and emission wavelengths are mm−1, mm−1, mm−1, and mm−1 respectively, which were determined by diffuse optical tomography. There is a small hole with a 1mm radius and 2 mm height, which was used to emplace the 3804.4 nM Cy5.5 solution (with the quantum efficiency of 0.23 at the peak excitation wavelength of 671 nm and the extinction coefficient of about 0.019 mm−1 [24]). The fluorophore was excited by 8 point sources at mm plane, and there are 160 detectors being located on the boundary inside 90° FOV corresponding to each point source, as shown in Figures 6(d) and 6(e).
(a)
(b)
(c)
(d)
(e)
3.2.2. Reconstruction for Experimental Data
The cubic in Figure 7 was discretized into 3191 nodes and 15870 tetrahedral elements for the inverse problem. The reconstruction process was the same as in Section 3.1.2. Figure 7 shows the recovered results by our method. The reconstructed center coordinate of Cy5.5 is 14.0 mm, 9.1 mm, and 7.7 mm while the real center coordinate is 16.0 mm, 8.0 mm, and 8.5 mm with location error of 2.42 mm. The recovered concentration of Cy5.5 dye is about 3478.3 nM with a relative error of 8.6%. Figure 8 shows the corresponding results by algebraic reconstruction technique (ART) and -norm regularization method. It is shown that their fluorescence yields were smaller than the actual one.
(a)
(b)
(c)
(a)
(b)
4. Discussion and Conclusion
The quantitative reconstruction of fluorescence molecular tomography (FMT) is meaningful for biological applications. In this paper, we developed a two-stage reconstruction method based on normalized Born approximation for quantitative fluorescence tomography. This algorithm includes two steps. The first step is to localize the permissible region of fluorophore. In this stage, the fluorescent target is very small compared with the whole reconstruction region, which can be considered as a priori information in the form of -norm penalty. In the second step, algebraic reconstruction technique (ART) is used to recover the fluorescence yield in the permissible region of fluorophore, which was provided by the first step.
Both the numerical experiment and physical experiment were conducted to validate the feasibility and potential of the proposed algorithm. In all experiments, our method was compared with ART and -norm regularization method, and the developed algorithm showed a better performance. In the physical phantom experiment, the differences in the positions of reconstructed fluorescent target and actual one for all methods were obvious. There are many possible reasons for this phenomenon. First, the strong tissue scattering makes the optical imaging in thick tissue with low-resolution. Second, the ill-posed and ill-conditioned nature of the inverse problem makes the solution unstable and nonunique. Therefore, the recovered results are easily affected by the measurement noise, which is inevitable during the signal detection. Last but not least, the system errors caused by geometry mismatch and optical modeling approximation might also lead to large errors in the reconstruction algorithm [25].
Generally speaking, compared with conventional methods, our method is easy to use with high accuracy. The outstanding advantage of this method is that it only utilizes the optical information to provide the permissible region of target. In [10, 13], the authors used X-ray computed tomography (XCT) to provide an accurate exterior boundary of object and anatomical images of the animal as a priori information, on which the fluorescent target can be localized. In some clinical problems, however, the targets cannot be obtained by XCT, like tumors or anomaly tissues. In this case, our method showed its advantages, as we employ the -norm regularization to localize the object instead of XCT. The weakness of the proposed approach is that it would not perform well when localization of object cannot be obtained accurately in the first stage. If the target is localized with a large error in the first stage, the conduction of the second stage would become of no meaning.
In conclusion, we have proposed an effective method for quantitative fluorescence tomography. Numerical experiments and physical experiments were both conducted to validate the feasibility and potential of our method. This algorithm would be applied in in vivo small animal tumor experiments in the future.
Acknowledgments
This work was supported by the Program of the National Basic Research and Development Program of China (973) under Grant no. 2011CB707702, the National Natural Science Foundation of China under Grant nos. 81090272, 81101083, 81101084, 81101100, 81000632, and 30900334, the National Key Technology Support Program under Grant no. 2012BAI23B06, and the Fundamental Research Funds for the Central Universities.