Open Access
Add to BibSonomyAbstract
Bioluminescence tomography (BLT) has become a powerful tool for whole-body small animal imaging. In this contribution, an adaptive improved element free Galerkin method (IEFGM) is presented to perform a quantitative reconstruction of the internal light source using quasi- or multi-spectral information, which not only can avoid the time-consuming mesh generation but also can reduce the ill-posedness of BLT effectively. In the algorithm, the reconstruction can be largely enhanced by an adaptive technology based on a posteriori error estimation. Finally, the numerical and physical phantom experiment results show that the bioluminescent source can be recovered accurately.
©2009 Optical Society of America
1. Introduction
As a promising molecular imaging modality, bioluminescence tomography (BLT) has been widely studied for its excellent performance and high cost-effectiveness in recent years [1, 2]. The goal of BLT is to reconstruct the internal light source with the transmitted and scattered bioluminescent signal detected on the external surface of a small animal, which can promote tumorigenesis study, drug development, gene therapy, etc. [1, 3, 4].
At present, permissible source region strategy or multi-spectral information is commonly used to reduce the ill-posedness of BLT [5, 6, 7, 8, 9]. However, the determination of permissible source region requires the surface light power distribution, the anatomical structure and the peak intensity attenuation of photon transportation [5, 10]. Furthermore, the permissible source region can not always be inferred, especially for the deeper and/or multiple sources cases. In addition, the size of permissible source region has an important effect on numerical stability and efficiency [5]. On the other hand, multi-spectral data acquisition becomes easier with novel multiple luciferase-tagged technique and advanced spectral separation methods [11, 12]. With multi-spectral information, the ill-posed inverse problem of BLT will become well-posedness, and the solution existence and uniqueness of multi-spectral BLT have also been demonstrated in the reference [13]. Thus, multi-spectral BLT has attracted increasing attention and a consensus that spectrally resolved method can enhance the uniqueness and stability of BLT solution has been achieved [7, 8, 9]. In previous studies, multi-spectral information is commonly captured by the band-pass filter [7, 8], whose transmissivity is too low for the weak boundary bioluminescent signal. From the spectrogram of cutoff filter, it is conceivable that the photons, whose wavelength is higher than the cutoff threshold of the high-pass filter, can pass through the filter. Comparatively, only light that locates in the bandwidth of band-pass filter can reach the detector, which is much fewer than the quasi-spectral case. Moreover, the band-pass filter is more expensive. Therefore, quasi-spectral information acquired using several cutoff filters in physical phantom experiment is employed for bioluminescent source reconstruction in this paper for the first time to our best knowledge, which not only has higher transmissivity and signal-to-noise ratio (SNR) than band-pass filter, but also reduces the imaging system cost.
Generally, the finite element method (FEM) is the most popular numerical technique in the inverse source reconstruction. However, the geometric model should be discretized into finite element meshes in FEM analysis, and it is a difficult and time-consuming process [3]. Comparatively, meshless methods use only a set of points and do not require any node connectivity or element information, which can avoid the burdensome mesh generation [14, 15, 16]. Moreover, the reconstruction can be improved greatly using an adaptive technology based on a posteriori error estimation. Thus, this contribution presents an adaptive improved element free Galerkin method (IEFGM) for bioluminescence tomography using quasi- or multi-spectral measurement, and the algorithm performance is validated with numerical and physical phantom experiments.
2. Methodology
2.1. Diffusion equation and Robin boundary condition
In BLT, the propagation of bioluminescent photons in tissue can be modeled by the steady-state diffusion equation and Robin boundary condition under the assumption that light scattering dominates over absorption and the internal source operates in continuous wave (CW) mode during image acquisition [1, 17]. Considering the influence of spectrum λ, the above equations can be expressed as:
where Ω and ∂Ω are the region of interest and the corresponding boundary; Φ(x,λ) denotes the photon flux density [W/mm2]; S(x,λ) represents the internal source density [W/mm3]; µa(x,λ) is the absorption coefficient [mm-1]; D(x,λ)=(3(µa(x,λ)+µ′s(x,λ)))-1 is the optical diffusion coefficient, µ′s(x,λ) the reduced scattering coefficient [mm-1]; v(x) refers to the unit normal vector outward to the surface ∂Ω; A(x;n,n′) is a function to incorporate the mismatch between the refractive indices n within tissue and n′ in the surrounding medium.
2.2. Source reconstruction algorithm
The conventional element free Galerkin method (EFGM) is based on the moving least squares (MLS) approximation, which is constructed by three components: a weight function related to each discretized point, a basis function and a series of non-constant coefficients [14, 15, 16].
where Ni(x) denotes the MLS shape function on the ith point; p(x) is the monomial basis function, p T(x)=[1,x,y,z,x 2,xy,y 2,yz,z 2,zx]; , Nn the number of nodes and w(x-x i) the weight function; The quartic spline function is specified as the weight function in this contribution [14, 16]; H i(x) stands for the ith column of matrix However, the MLS interpolation does not pass through the nodal function values, so the approximation should be performed twice for the actual solution. Furthermore, the Dirichlet and Robin boundary condition can not be imposed directly. In this paper, the MLS shape function is modified to possess the Kronecker delta function property [14, 16]:
where Ml(x) represents the modified MLS shape function on the lth node; Nl(x i)-1 is an element in the inverse of transformation matrix Γ[14]:
In multi-spectral data analysis, the bioluminescence spectrum can be divided into several bands [λk,λ k+1],k=1,2,⋯τ, and the energy contribution of the internal light source in the spectral band [λk,λ k+1] is ωk, satisfying ωk≥0 and ∑τ k=1 ωk≈1, which can be determined according to the literature [11]. Comparatively, the energy percentages ωk in [λk,∞) in the quasi-spectral framework appear in non-increasing order: 1=ω 1>ω 2>⋯>ωτ>0. In the spectral band [λk,λ k+1] or [λk,∞), the photon flux density Φ(x,λk) and the light source density S(x,λk) can be approximated with their corresponding nodal values:
Based on Galerkin method and Gauss theory, Eqs. (1) and (2) can be transformed to the following weak form:
Substituting Eq. (6) into Eq. (7), and then the matrix equation in the spectral band [λk,λ k+1] or [λk,∞) can be obtained as follows:
where M k is a symmetric positive definite matrix. Considering the photon density Φmeask measured on the side surface, the following linear system can be established:
where A k can be obtained through retaining the rows of M -1 k F k corresponding to the boundary detected points. Taking all the spectral bands into account, we have
where
In order to recover bioluminescent source distribution and keep the uniqueness and stability of the BLT solution, the optimization objective function is defined using Tikhonov regularization method:
where S sup is the upper bound of the source density; L stands for the weight matrix, ‖V‖Λ=V TΛV; α is the regularization parameter. In this contribution, L equals to the unit matrix and the regularization parameter can be determined by L-curve method [10, 18]. A modified Newton method and a specific Hessian matrix are adopted to solve the above minimization problem with simple bounds for the optimal regularization parameter, and then the localization and quantification of the reconstructed bioluminescent source can be determined [19, 20].
In the framework of IEFGM algorithm, if node distribution in the actual source region is sparse, distortion in the reconstructed source quantification and localization will become more significant. Therefore, an adaptive technology based on a posteriori error estimation is employed for the above issue. After each level reconstruction, node refinement will be carried out using the direct maximum selection method until the stopping criterion ‖Φmeas-Φcomp‖<ε is met, where Φcomp denotes the surface photon density computed with the reconstructed source and ε is the stopping threshold. If the above stopping criterion is not met, the procedure will select a series of points with larger reconstructed values for node refinement. It is easily conceivable that the nodes with higher values of the optimization results most likely represent actual source locations. After node refinement, the reconstruction can be further improved. To sum up, the corresponding flowchart of proposed algorithm is showed in Fig. 1.
3. Experiments and results
In imaging experiments, a cubic tissue-like heterogeneous phantom of 11 mm side length was designed to evaluate our algorithm, and two same light sources with 0.125 mm3 volume and 1.0 nW/mm3 power density were embedded into the phantom with their centers at (3.25,3.25,5.25) and (7.25,7.25,5.25) respectively, as shown in Fig. 2(a). According to the experimental results of the bioluminescent reporters in the reference [11], we partitioned the spectrum into three wavelength ranges: [400,530]nm, [530,630]nm and [630,750]nm, and the corresponding energy percentages of each spectral band were specified as ω 1=0.29,ω 2=0.48 and ω 3=0.23. In addition, the optical properties at the emission wavelength range between 400 nm and 750 nm are listed in Table 1, which could be obtained from the document [21]. In the adaptive IEFGM analysis, 12×12×12 uniform initial node distribution was arranged. Figure 2(b) presents the reconstruction result without adaptive strategy, and the flux densities of the recovered light sources were 0.125 nW/mm3 and 0.124 nW/mm3 respectively. However, the reconstructed power densities could reach 0.889 nW/mm3 and 0.924 nW/mm3 separately after node refinement, as shown in Fig. 2(c). The corresponding quantitative comparison is summarized in Table 2, and RE in the table represents the relative error between the actual source and the reconstructed result, which can be computed with the formula RE=|Sreal-Srecons|/Sreal,Sreal the real internal bioluminescent source, and Srecons the reconstructed source.
Table 1. Optical properties of the numerical phantom in different multi-spectral ranges.
Then, the above bioluminescent sources were also recovered using quasi-spectral information. According to the separation of aforementioned multi-spectral bands, the quasi-spectrum was divided into three bands: [400,∞)nm, [530,∞)nm and [630,∞)nm, and the corresponding energy contributions in each band were set as: ω 1=1.00, ω 2=0.71 and ω 3=0.23. The optical parameters in each spectrum are summarized in Table 3. The reconstructed power densities were 1.078 nW/mm3 and 0.862 nW/mm3 respectively with adaptive technology, and the centers of the recovered sources located at (3.0,3.0,5.0) and (7.5,7.5,5.5) separately. The reconstructed results are depicted in Fig. 3, and the corresponding quantitative comparison is listed in Table 4. From Table 2 and 4, we can see that the quasi-spectral information is equivalent to multispectral measured data for source reconstruction.
In addition, a polyoxymethylene cube phantom with 20 mm side length was fabricated for further verifying the feasibility of the proposed method, as presented in Fig. 4(a). A hole of 1.25 mm radius was drilled in the phantom to inject luminescent solution used as the light source, which contained peroxide, ester compound and fluorescent dye. The source power density was
Table 2. Comparison between the actual and reconstructed sources in multi-spectral case.
Table 3. Optical properties of the numerical phantom in different quasi-spectrums.
273.89 nW/mm3 measured by a integrating sphere and the main emission spectrum located about between 600 nm and 800 nm. In this experiment, the used imaging system mainly contains a liquid-nitrogen-cooled back-illuminated CCD camera (Princeton Instruments VersArray 1300B), a Nikkor 17–55 mm f/2.8 lens, a rotation stage, a light-tight enclosure (not in this figure) and a group of translation stages, as shown in Fig. 4(b), which can accomplish non-contact light detection. Figure 4(c) depicts the middle cross-section of the phantom, and the four degrees show the direction of the CCD camera during data acquisition. The source spectrum was split into three bands using a group of cutoff filters: [600,∞)nm, [650,∞)nm and [700,∞)nm, and the energy contributions in each spectral band were measured by a spectrometer as follows: ω 1=1.00, ω 2=0.64 and ω 3=0.16. Over the above bands, the absorption coefficient and the reduced scattering coefficient were determined by a time-correlated single photon counting (TCSPC) system [22], as listed in Table 5. Using single wavelength data, the density of the recovered source was 305.78 nW/mm3, and the relative error was 11.64%, as depicted in Fig. 5(a). For two spectrum bands case in Fig. 5(b), the reconstructed flux density was 296.44 nW/mm3 with the relative error about 8.23%. Finally, the source density and the relative error
Table 4. Comparison between the actual and reconstructed sources in quasi-spectral case.
of the reconstructed source were 289.63 nW/mm3 and 5.75% respectively using three bands information, as shown in Fig. 5(c).
Fig. 2. Numerical phantom experiment using multi-spectral data. (a) A heterogeneous phantom with two light sources; (b) Transverse view of the reconstructed result without adaptive technique; (c) The recovered sources after node refinement. The black square indicates the actual source position.
Fig. 3. Numerical phantom experiment using quasi-spectral information. (a) The reconstructed result without adaptive technique; (b) The recovered sources after node refinement.
Fig. 4. Imaging system and resinous phantom. (a) The polyoxymethylene phantom; (b) The non-contact imaging system; (c) The middle cross-section of the phantom.
4. Discussion and Conclusion
An improved element free Galerkin method (IEFGM) based on a posteriori error estimation adaptive strategy with quasi- or multi-spectral information has been proposed for BLT. Compared with conventional reconstruction algorithms, the above method has the following merits.
First, meshless method is applied to inverse bioluminescent source reconstruction for the first time to the best of our knowledge, which uses only a series of nodes without consideration
Table 5. Optical parameters for the physical phantom.
of the interrelationship among the points, so the troublesome mesh generation and data preprocessing can be avoided in comparison with FEM. Furthermore, the traditional MLS shape functions are modified to satisfy the Kronecker delta function property, with which the Dirichlet and Robin boundary condition can be imposed directly.
Fig. 5. Physical phantom experiment. (a), (b) and (c) The reconstructed source using single, two and three bands measurement respectively.
Secondly, in inverse source reconstruction, the internal light source is unknown for us. If we use the meshless method without adaptive technology, the bioluminescent source will not be recovered accurately, which has been demonstrated by the first numerical experiment results. Therefore, the proposed algorithm uses an adaptive technique based on a posteriori error estimation to eliminate the influence of the discretized nodes distribution, and the reconstruction is further improved.
Lastly, multi-spectral information has been widely used in BLT. However, the spectrum measurements are commonly obtained by the band-pass filter in the previous studies, whose transmissivity is too low for the weak surface bioluminescent signal. Comparatively, the transmittance of the cutoff filter is higher and the price is lower. To our best knowledge, quasi-spectral data detected with the cutoff filter is utilized for source reconstruction for the first time, which significantly reduces the ill-posed property of BLT. In practice, the number of spectrums is determined by the number of surface measurable points and the number of internal reconstructed nodes.
In conclusion, we have presented an adaptive improved element free Galerkin method based on quasi- or multi-spectral measurement for inverse bioluminescent source reconstruction. The numerical and physical phantom experimental results demonstrate the feasibility of the proposed algorithm, and the source location and density can be recovered exactly. Our future work will focus on testing the proposed algorithm using the real small animal in vivo experiment data. The corresponding results will be reported later.
Acknowledgments
This work has been supported by the Project for the National Basic Research Program of China (973) under Grant No. 2006CB705700, Changjiang Scholars and Innovative Research Team in University (PCSIRT) under Grant No. IRT0645, CAS Hundred Talents Program, CAS scientific research equipment develop program under Grant No. YZ200766, the Knowledge Innovation Project of the Chinese Academy of Sciences under Grant No. KGCX2-YW-129, the National Natural Science Foundation of China under Grant No. 30672690, 30600151, 60532050, 60621001, 30873462, 60910006, 30970769, 30970771, Beijing Natural Science Fund under Grant No. 4071003, Technology Key Project of Beijing Municipal Education Commission under Grant No. KZ200910005005.
References and links
1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weisslder, “Looking and listening to light: the evolution of whole body photonic imaging,” Nat. Biotechnol. 23, 313–320 ( 2005). [CrossRef] [PubMed]
2. J. Tian, J. Bai, X. P. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, “Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. 27, 48–57 ( 2008). [PubMed]
3. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 ( 2005). [CrossRef] [PubMed]
4. G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. 31, 2289–2299 ( 2004). [CrossRef] [PubMed]
5. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13, 6756–6771 ( 2005), http://www.opticsinfobase.org/abstract.cfm?URI=oe-13-18-6756. [CrossRef] [PubMed]
6. J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16, 15640–15654 ( 2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-20-15640. [CrossRef] [PubMed]
7. H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. 31, 365–367 ( 2006), http://www.opticsinfobase.org/abstract.cfm?URI=OL-31-3-365. [CrossRef] [PubMed]
8. Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finite element analysis: methodology and simulation,” Phys. Med. Biol. 52, 4497–4512 ( 2007). [CrossRef] [PubMed]
9. A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. 50, 5421–5441 ( 2005). [CrossRef] [PubMed]
10. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14, 8211–8223 ( 2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-18-8211. [CrossRef] [PubMed]
11. H. Zhao, T. C. Doyle, O. Coquoz, F. Kalish, B. W. Rice, and C. H. Contag, “Emission spectra of bioluminescent reporters and interaction with mammalian tissue determine the sensitivity of detection in vivo,” J. Biomed. Opt. 10, 041210 ( 2005). [CrossRef]
12. G. Wang, H. Shen, Y. Liu, A. Cong, W. Cong, Y. Wang, and P. Dubey, “Digital spectral separation methods and systems for bioluminescence imaging,” Opt. Express 16, 1719–1732 ( 2008),http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-3-1719. [CrossRef] [PubMed]
13. W. M. Han and G. Wang, “Theoretical and numerical analysis on multispectral bioluminescence tomography,” IMA J. Appl. Math. 72, 1–19 ( 2006). [CrossRef]
14. C. Qin, J. Tian, X. Yang, K. Liu, G. Yan, J. Feng, Y. Lv, and M. Xu, “Galerkin-based meshless methods for photon transport in the biological tissue,” Opt. Express 16, 20317–20333 ( 2008), http://www.opticsinfobase.org/abstract.cfm?URI=oe-16-25-20317. [CrossRef] [PubMed]
15. J. Dolbow and T. Belytschko, “An introduction to programming the meshless element free Galerkin method,” Arch. Comput. Methods Eng. 5, 207–241 ( 1998). [CrossRef]
16. X. Zhang and Y. Liu, Meshless methods, (Tsinghua University Press, Beijing, 2004).
17. 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 conditions,” Med. Phys. 22, 1779–1792 ( 1995). [CrossRef] [PubMed]
18. P. R. Johnston and R. M. Gulrajani, “Selecting the corner in the L-curve approach to Tikhonov regularization,” IEEE T. Bio-Med. Eng. 47, 1293–1296 ( 2000). [CrossRef]
19. P. E. Gill, W. Murray, and M. Wright, Practical optimization, (Academic Press, New York, 1981).
20. C. Qin, J. Tian, Y. Lv, and W. Yang, “Three-dimensional bioluminescent source reconstruction method based on nodes of adaptive FEM,” Proc. SPIE 6916, 69161K ( 2008). [CrossRef]
21. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50, 4225–4241 ( 2005). [CrossRef] [PubMed]
22. D. Qin, H. Zhao, Y. Tanikawa, and F. Gao, “Experimental determination of optical properties in turbid medium by TCSPC technique,” Proc. SPIE 6434, 64342E ( 2007). [CrossRef]
