Abstract
Bioluminescence tomography (BLT) is a molecular imaging modality which can three-dimensionally resolve the molecular processes in small animals in vivo. The ill-posedness nature of BLT problem makes its reconstruction bears nonunique solution and is sensitive to noise. In this paper, we proposed a sparse BLT reconstruction algorithm based on semigreedy method. To reduce the ill-posedness and computational cost, the optimal permissible source region was automatically chosen by using an iterative search tree. The proposed method obtained fast and stable source reconstruction from the whole body and imposed constraint without using a regularization penalty term. Numerical simulations on a mouse atlas, and in vivo mouse experiments were conducted to validate the effectiveness and potential of the method.
1. Introduction
Due to its ability of monitoring physiological and pathological activities at the molecular level, small animal optical molecular imaging has become an important method for biomedical research. Bioluminescence imaging, as one of optical molecular imaging modalities, has attracted attention for its advantages in sensitivity, specificity, and cost effectiveness in cancer research and drug development [1–3]. Compared with planar bioluminescent imaging, BLT can three-dimensionally reconstruct the bioluminescent probes in small animals [4].
The generic BLT model is ill-posed. However, it has been theoretically proven that the solution uniqueness can be established under practical constraints using a priori knowledge [5]. In most existing reconstructions, multi spectral measurement [6–10], permissible source region (PSR) [9–12], and sparse reconstruction [13–15] are three common strategies to reduce ill posedness of BLT. Although multi spectral techniques improve reconstruction qualities to a certain degree by increasing the measurable data, they in turn impose some limitations in practical applications such as increased signal acquisition time and a high computational cost [16]. Besides, the PSR strategy can significantly improve the location accuracy of reconstructed source and reduce the computational cost by limiting the reconstruction region into a small area. However, in practical applications, both the size and position of the permissible region have significant impact on imaging results [5]. Additionally, since the bioluminescent source distribution is usually sparse in practical applications and only insufficient boundary measurements are available, the compressed sensing can bring benefits in spatial resolution and algorithm stability to BLT reconstruction. Recently, many sparse reconstruction methods have emerged in BLT [17–19]. The majority of them reformulate the BLT inverse problem into minimizing an objective functions that integrate a sparse regularization term with a quadratic error term and solve it via computationally tractable convex optimization methods, such as linear programming and gradient methods, However, the efficacy of the sparse regularization methods strongly depends on the choice of regularization parameter in practical applications [19].
The existing works have demonstrated that PSR can improve reconstruction qualities by reducing the number of unknown variables. Although the meaningful results can be obtained by using the PSR that is based on the bioluminescent signals and a priori knowledge available from a specific biomedical application [4, 11, 12], in most cases, it is rather difficultly to manually select such a small and appropriate region. Recently, some optimal permissible source region methods have emerged in BLT. Feng et al. presented a reconstruction algorithm for a spectrally resolved BLT based on an adaptive rough estimate of an optimal permissible source region and multilevel finite element method approach (FEM), where Tikhonov regularization was used to solve the constrained BLT inverse problem [9]. Naser and Patterson proposed a two-step reconstruction algorithm of bioluminescence, in which the permissible regions were shrunk by using an iterative minimization solution based on the L1 norm [10]. These works demonstrated the feasibility and potential of the optimal PSR techniques with numerical simulation. However, both of two previous reconstructions also needed the regularization methods to reconstruct the sources in the allowed region, which made the reconstruction results also depend on the choice of regularization parameter. Furthermore, they were demonstrated with only regular phantoms simulations and presented no in vivo experiment validation.
In this work, a novel BLT reconstruction algorithm based on the semi greedy method was proposed. The optimal PSR problem was cast into a search for the correct support of source distribution among a number of dynamically evolving candidate subsets, and the optimal PSR was chosen automatically by using an iterative search tree. Therefore, the columns of the system matrix were treated as the nodes for building up the search tree where each path from the root to a leaf node denoted a candidate. The search tree was initialized with some unspecific nodes. At each iteration, new nodes were appended to the most promising path, which were selected to minimize the cost function based on the residue. The permissible source region was expended by adding nodes with high a probability to contribute to the source. Among the system matrix, the columns that were corresponding to the nodes contained in the most promising path were selected to obtain the source distribution. By automatically choosing an optimal PSR, the method reduced the ill posedness of the problem and imposed constraint without using a regularization penalty term.
This paper is organized as follows. In Section 2, the forward photon propagation model, the inverse problem for BLT with FEM, and the proposed algorithm are introduced. In Section 3, the numerical simulations in a mouse atlas demonstrate the performance of the proposed method. In Section 4, an in vivo mouse experiment is conducted to further evaluate its reliability. Finally, we discuss the results and conclude this paper.
2. Method
Light propagation in biological media is essentially the transport of radiant energy. The radiative transfer equation (RTE) can rigorously describe light transport in turbid media [20]. Compared with the other approximations of RTE such as simple spherical harmonics, spherical harmonics and discrete ordinates, the following steady-state diffusion approximate equation (DA) in (1) is the most popular one as a result of its moderate computational efficiency and explicit physical meaning [4, 5, 11, 15, 21]. Assuming that the bioluminescence imaging experiment is performed in a totally dark environment and no photon travels into through the boundary , the equation is subject to the Robin boundary condition in (2) as follows: where is the domain of the problem, donates the source energy distribution, represents the photo fluence rate, is the absorption coefficient, is the reduced scattering coefficient, indicates the optical diffusion coefficient, donates the boundary of the problem, and represents the mismatch coefficient between and its surrounding medium. The measured quantity on the boundary is given by the outgoing radiation as follows:
FEM is a powerful tool for solving the DA equation [4, 6–10]. By using FEM to discretely approach the solving domain and making a series of transformations and rearrangements, the linear relationship links the source distribution inside the heterogeneous medium, and the photon fluence rate on the surface is established as follows: where is the source distribution of the interior nodes, is the measurable photon flux photon on the boundary nodes, is the positive definite matrix, and is the source weight matrix. The nonmeasurable entries in and corresponding rows in can be removed. Then a new linear relationship can be obtained as follows [4, 22]:
For BLT, the domains of the bioluminescent sources are usually very small and sparse compared with the entire reconstruction domain. That means that there are only fewer nonzero components in . Therefore, the system matrix can be seen as the dictionary, and has a less-term representation over the dictionary. As a consequence, the aim of the proposed method is imposing the constraint on the source space by choosing only the part that contributes to the source distribution. In the language of sparse approximation, greed pursuit algorithms are useful methods for solving this problem [23, 24]. For instance, Orthogonal Matching Pursuit (OMP) is to pick columns in a greedy fashion [25]. According to the introduction of OMP, the reconstruction starts with an empty index set. At each iteration, we choose the single column of that is most strongly correlated with the remaining part of . Then, we subtract its contribution to and iterate on the residual. The reconstruction is stopped after a number of iterations. Unfortunately, from experimental results shown in Figure 2 and Table 2, we found that OMP, as the single-path algorithm, could not achieve the desired expectations for reconstructing the bioluminescent source. In the experiments, when the computation of a single path selects a wrong column, the correct one is still in the set of candidate representations. Therefore, incorporation of a multipath search strategy is motivated to improve reconstruction. In this section, the semigreedy method was used to search for the correct support of among a number of dynamically evolving candidate subsets.
A general best first (GBF) is a search algorithm which constructs a tree by expanding the most promising node chosen according to a specified rule. Search algorithm , as one of the most studied versions of GBF, can find path in combinatorial search. It selects an optimal path by minimizing an additive evaluation function , where is the cost of the currently evaluated path from start node to , and is a heuristic estimate of the cost of the path remaining between and some goal nodes [26–29]. In our problem, the search tree was iteratively built up by nodes which represent the dictionary atoms. Each path from the root to a leaf node denoted a subset of dictionary atoms which was a candidate support for . Let us define the notation. and denote the atoms contained in path and the vector of corresponding coefficients obtained after orthogonal projection of the residue onto the set of selected columns [25]. Similarly, and represent the selected atom at the th node on path and corresponding coefficient.
The search tree starts with less unspecified nodes. A simple way is selecting the nodes that have the highest absolute inner product with . In order to find the fewest possible nodes, the search must constantly make an evaluation to decide which available paths should be expended next. Therefore, the evaluation function is defined as follows:
Beside the evaluation function, the auxiliary function is also needed to assess the cost brought by adding a preferred goal node to the path. Generally, according to the expectation on average equal contribution of unopened nodes, the auxiliary function can be built as follows: where coefficient is defined as . is a ratio between the number of the nonzero entries and the zero entries in the solution . It is well known that the sparse solution has only less nonzero entries. Therefore, in most practical application, is very smaller than 1. Here, we selected . If the source distribution could be seen as the K-sparse signal, could be computed by . The cost function can be written as follows: where is a constant. A lot of experiments for different reconstruction models including 2D and 3D experiments were performed to evaluate the impact of on the source reconstruction. We found that better results could be obtained when it varied in the interval [1.00, 1.25]. Therefore, in our experiments, its value was taken from the interval and set to be1.05.
In practice, if all of the children of the most promising partial path are added to the search tree at each iteration, the tree might have too many search paths. Therefore, the following pruning strategies are employed a guide on how the tree grows.
The first one is about extensions per path. At each step, it is not necessary to have all of the unopened atoms added to the current optimal path. We can expand the search tree only by the children which have the highest absolute inner-product with the residue to the selected path. This pruning strategy decreases the upper bound to on the number of paths. Practically, and are selected to be much smaller than , which can drastically decrease the paths involved in the search. Although the number of extensions per path is limited to , it is also necessary to control the size of path. That is because that adding new paths at each iteration continues increasing required memory. Therefore, to reduce memory requirements, we adopted the “beam search” strategy [30] and limit the maximum number of paths in the tree by the beam width . When this limit is exceeded, the paths with maximum cost are seen as the worst paths and are removed from the tree until paths remain. Here, and were set to be 4 and 200, respectively. Moreover, since order of nodes along a path is unimportant, amalgamating the equivalent path is also important to improve the search performance. For this purpose, we define a path equivalency notion; and are two paths with different length . If all atoms of can be found in and these composed the continuous subset of , we define the above two paths as being equivalent. Consequently, the insertion of into the tree is unnecessary.
After the growing of the search tree is finished, the linear relationship between the observation and the selected PSR can be established as follows:
Since the nodes contained in the optimal path can be much smaller than the number of all nodes , (9) is an overdetermined linear equation. Therefore, a limited memory variation of the Broyden Fletcher Goldfarb Shanno (LBFGS) [31] was used to directly solve (9).
3. Simulation Studies in the Mouse Atlas
In this subsection, heterogeneous simulations were presented to demonstrate the performance of the proposed method for mouse applications. The experimental data were acquired by a dual-modality BLT/micro-CT system developed in our lab [32, 33]. By using image processing and interactive segmentation technology, heterogeneous model including heart, lungs, liver, bone, spleen, and muscle was built. The optical coefficients for each organ are listed in Table 1 [34]. Here, the torso section with a height of 25 mm was selected as the reconstruction region.
3.1. Reconstruction in a Single Source Case
In the experiment, a spherical source with a 0.6 mm diameter was placed in the liver with the center at (18.24 mm, 31.58 mm, 47.29 mm) as shown in Figure 1(a). The source was modeled as isotropic point sources whose strength was set to be 2 nW/mm3. As for the forward problem, the FEM was used to obtain the synthetic measurements on the boundary. The atlas model was discretized into a tetrahedral-element mesh with 30892 nodes and 167841 elements. The generated simulated photon distribution on the boundary is presented in Figure 1(b). Then the forward solutions were projected onto a single coarse mesh consisting of 20068 elements and 3098 nodes, which was used for reconstructing the source.
(a)
(b)
(a)
(b)
(c)
(d)
To better illustrate the performance of the proposed method, we compared the proposed method with OMP [25] and FIST-L1 [35, 36]. The former is a typical greedy pursuit method for sparse signal recovery. The latter, as a sparse regularization method, can be viewed as a standard approach to ill-posed linear inverse problems and has been used in fluorescence molecular tomography (FMT) and BLT. Here, the step size in FIST-L1 was computed by using the estimation algorithm introduced in [35]. Since the regularization parameter plays an important role in the regularization methods, we performed two experiments with different regularization parameters that were set to be 4e-11 and 4e-10, respectively. All of the reconstructions were carried out on a personal computer with 3.2 GHz Intel Core2 duo CPU and 2 GB RAM.
The qualities of the reconstruction were quantitatively assessed in terms of location error and the maximum reconstructed intensity. The location error was defined as Euclidean distance between and , where and were the real locations of the source center and the location of the node with the maximum reconstructed value, respectively. The visual effects of the reconstruction results are presented in the form of slice images and iso-surfaces, as shown in Figure 2. Additionally, the detailed information about parameters and the final quantitative reconstruction results are summarized in Table 2.
We found that the reconstructed positions by L1 regularization with an optimal regularization parameter and the proposed method were identical. Specifically, the reconstructed center was (18.26 mm, 31.97 mm, 47.28 mm) with a location error of 0.3995 mm from the actual source, whereas the location error by OMP was 2.9605 mm. The performance of OMP was inferior to the other two methods. L1 regularization performed well and obtained satisfactory source localizations and maximum reconstruction value. However, the selection of the regularization parameter had a great impact on the reconstruction results. As for the proposed method, it performed slightly better than L1 regularization with manually optimized regularization parameter in terms of maximum reconstruction value. Moreover, it was also an efficient reconstruction method.
The above experiments were performed without noise. In order to evaluate the sensitivity of the proposed method to various noise levels, six cases were carried out where the measurements were added to 5%, 10%, 15%, 20%, 25% and 30% Gaussian noise, respectively. We also made a comparison between the proposed method and FIST-L1 with the regularization parameter that was set to be 4e-11. The reconstruction results under different noise levels are compiled in Table 3 and Figure 3 showing that the proposed method was robust to measurement noise.
(a)
(b)
(c)
(d)
(e)
(f)
3.2. Double Source Case and MultiSource Case
Dual source setting was also considered in order to evaluate the proposed method. Two sources have the same size as one used in the single source case. Their strength and position were set to be 2 nW/mm3 and (18.24 mm, 31.58 mm, 47.29 mm), 1 nW/mm3 and (18.74 mm, 39.15 mm, 47.06 mm), respectively. The reconstruction results are shown in Figure 4 and Table 4.
Multiple sources setting simulation experiment was also presented to further demonstrate the ability of the proposed method. Based on the setting in double source case, the third source with the same size and shape was added. Its strength and position were set to be 1.5 nW/mm3 and (23.60 mm, 37.94 mm, 47.45 mm). The final reconstruction results are presented in Figure 5 and Table 5. The result of two group experiments indicated that the sources can be accurately distinguished by using the proposed method.
4. In Vivo Experiment Validation
Besides the numerical simulations with mouse atlas, an in vivo experiment was carried out on a mouse to further test the proposed method. The experiment was also performed with a dual-modality BLT/micro-CT system developed in our lab [32, 33, 37]. A nude hairless mouse (Nu/Nu, Laboratory Animal Center, Peking University, China) was used in this experiment. To simulate a known bioluminescent source, a home-made cylindrical light source about 3 mm long and 1.5 mm in diameter was implanted into the abdomen of the mouse in this experiment. The source was made of a catheter filled with luminescent liquid and emitted a red luminescent light that had a similar emission spectrum with a firefly luciferase-based source.
Before the beginning of the experiment, the CCD (VersArray, Princeton Instruments, Trenton, NJ, USA) was cooled to 110° by using liquid nitrogen. The mouse was anesthetized and placed in a mouse holder. The mouse holder was set to rotate to 0°, 90°, 180° and 270°. At each of four positions, the mouse was photographed by the CCD camera. After the optical data were acquired, the mouse was scanned by using the micro-CT to obtain the anatomical maps which could provide structural information for the source reconstruction. Then the CT data were segmented into five regions represent muscle, lungs, heart, liver, and kidneys, respectively, as shown in Figure 6(). The heterogeneous model including five tissues was discretized into the mesh containing 11917 tetrahedral elements and 2557 nodes. The optical parameters for different tissues were calculated based on the literature as listed in Table 6 [34, 38]. The optical data was registered with the volumetric mesh, and measured data were mapped onto the surface of the mesh. The result of mapped photon distribution is shown in Figure 6().
It took about 8 seconds to complete the reconstruction using the proposed method. The final results are presented in Figure 7, where the reconstruction source center is (37.17 mm, 38.82 mm, and 20.92 mm) with a deviation about 2 mm to the actual center. As can be seen in the reconstruction results, the proposed methods could obtain satisfactory bioluminescent source localizations.
5. Discussion and Conclusion
In this paper, we have proposed a new method based on the semigreedy for bioluminescence tomography. The reconstruction results of the simulations on a mouse atlas demonstrate that the proposed reconstruction method is able to accurately and stably localize bioluminescent source from whole body, even with noisy measurements. The in vivo experiment further shows its performance.
The PSR strategy can significantly improve reconstruction qualities. However, in most cases, empirically selecting such small and appropriate region is unconvenient, even available. In this study, The optimal PSR problem is cast into a search for the correct support of source distribution among a number of dynamically evolving candidate subsets. In view of the characteristics in BLT sparse distribution, only the columns that contribute to the source reconstruction are chosen automatically by using semi-greedy method. The constraint imposed on the source space reduces the ill posedness of the problem and computational cost.
It is noted that in vivo experiment is not as accurate as simulations. Some reasons can be explained for this phenomenon. First of all, the error was generated, when the energy distribution was mapped from 2D images to a 3D mouse surface. Secondly, there were only five main segmented tissues used to build a heterogeneous model while others simply were regarded as the muscle, which also led to errors. Finally, the accuracy of the photon propagation model was very important for source reconstruction. The diffusion approximation was used due to its moderate computational efficiency and explicit physical meaning. However, it has some limitations in certain regions, such as void or more absorptive regions. Therefore, the error brought on by the DA model is inevitable. As discussed above, our future work will focus on studying more accurate forward models to describe photon propagation in biological tissues and improving the experimental procedures and imaging system to further promote the reconstruction quality.
Acknowledgments
This paper is supported by the National Basic Research Program of China (973 Program) under Grants no. 2011CB707700, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant no. KGCX2-YW-907, the National Natural Science Foundation of China under Grant Nos. 30970780, 81000624, 81027002, 81071205, 30970778, 81071129, and 81101095, the Beijing Natural Science Foundation No. 4111004, the Science and Technology Key Project of Beijing Municipal Education Commission under Grant No. KZ200910005005, the Doctoral Fund of the Ministry of Education of China under Grant No. 20091103110005, the Fellowship for Young International Scientists of the Chinese Academy of Sciences under Grant No. 2010Y2GA03, and the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists under Grant No. 2010T2G36.