1 March 2006 Boundary integral method for bioluminescence tomography
Author Affiliations +
J. of Biomedical Optics, 11(2), 020503 (2006). doi:10.1117/1.2191790
Bioluminescence tomography (BLT) allows in vivo localization and quantification of bioluminescent sources inside a small animal to reveal various molecular and cellular activities. We develop a reconstruction method to identify such a bioluminescent source distribution using the boundary integral method. Based on the diffusion model of the photon propagation in the biological tissue, this method incorporates a priori knowledge to define the permissible source region, and establish a direct linear relationship between measured body surface data and an unknown bioluminescent source distribution to enhance numerical stability and efficiency. The feasibility of the proposed BLT algorithm is demonstrated in heterogeneous mouse chest phantom studies.
Cong and Wang: Boundary integral method for bioluminescence tomography

Bioluminescent imaging has the capability to reveal molecular and cellular activities directly, and can be applied to all disease processes in most small-animal models. It is very sensitive, and helps diagnose diseases, monitor therapies, and facilitate drug development.1, 2 Bioluminescent imaging has been in a planar mode and largely a qualitative imaging tool.3 With bioluminescence tomography (BLT), 3-D localization and quantification is enabled of a bioluminescence source distribution inside a living small animal such as a mouse.4, 5

In bioluminescence imaging, photon scattering predominates over absorption in the biological tissue in this spectral range of interest. As a result, a significant number of bioluminescent photons can escape the attenuating environment, and they can be detected using a highly sensitive charge-coupled device (CCD) camera.3 In this case, the photons propagation in the biological tissue can be well described by steady-state diffusion equation and Robin-type boundary conditions.6, 7 The BLT principles and solution uniqueness for the BLT problem conditions were studied under some practical conditions. 4, 5, 8, 9 Some numerical algorithms for BLT were already reported, including the finite element based BLT algorithms, 10, 11, 12, 13, 14 BLT in combination with PET (OPET),15 multispectral bioluminescence optical tomography,16, 17 and BLT method based on diffusion theory of half infinite medium.18 In this letter, for the first time we develop a BLT algorithm using the boundary integral method. This methodology only requires finite element meshing of structural boundaries instead of complex volumetric finite elements previously used for BLT. Hence, the complexity and stability of the new BLT algorithm are improved as compared to that of the finite-element-based BLT algorithms.

The overall support region of a small animal Ω can be decomposed into a number of subregions Ωj (j=1,2,,τ) with Ω=j=1τΩj ; for example, lungs, heart, liver, bone, muscle, and so on. Applying the Gauss theorem, the steady-state diffusion equation can be transformed to a boundary integral equation on each subregion19:


where G(r,x)=exp(μeffrx)(4πDrx) is a Green function of the steady-state diffusion equation, D is the diffusion coefficient given by D=1[3(μa+μs)] , and μeff=[3μa(μa+μs)]12 where μa is the absorption coefficient (mm1) and μs is the reduced scattering coefficient (mm1) for the subregion Ωj . Equation 1 is a well-posed second kind integral equation. The bounding surface Ωj can then be split into N surface elements Γij (i=1,2,,N) on which the function Φ(x) and Φ(x)n are approximated by use of a set of p interpolation points and interpolation functions ϕk(x) on Γij 19:


The shape of the surface element can be arbitrarily selected, and usually made as a quadrilateral or a triangle. The number of points per surface element depends on the accuracy needed in the interpolation procedure. Let Φk=Φ(xk) and qk=DΦ(xk)n represent the values of the function Φ(x) and DΦ(x)n at an interpolation point xk , respectively. Since the BLT reconstruction is underdetermined and ill-posed, we incorporate a priori knowledge to define a permissible source region Ωjs(ΩjsΩj) for a subregion Ωj without loss of generality, where the bioluminescence source may be distributed. The region Ωjs can be discretized, and the embedded source function S(x) is approximated on Ωjs as


where sk=S(xk) represent the value of the source function S(x) at an interpolation point xk,ψk (x) is the interpolation basis function, and Ns is the number of discrete values of the source function.

Inserting Eqs. 2, 3 into Eq. 1, we obtain the following matrix equation:


where Φ=(Φ1,Φ2,,ΦM)T , Q=(q1,q2,,qM)T , M is the number of nodal points for Ωj , S=(s1,s2,,sNs)T , and Hj is a strictly diagonally dominant matrix, and allows its inversion. Multiplying Eq. 4 with the inverse Hj1 , we have


Furthermore, in terms of the continuity of Φ(x) and DΦ(x)n at the interface and the Robin-type external boundary condition, the above matrix equation for j=1,2,,τ can be assembled into19


where S denotes the source distribution in the permissible source region, and Φ a vector consisting of photon density values at the boundary and interface nodes, which can be divided into Φin at internal interface nodes and Φex at exterior boundary nodes. Because M is still a diagonally dominant matrix and invertible, we have B=M1F . After removal of those rows of B that correspond to Φin we obtain Bex and a linear relationship between measurable photon density at boundary nodes and the source distribution


Generally, the measured data in bioluminescence imaging are corrupted by noise, so it is not practical to solve for S directly from Eq.7. Instead, an optimization procedure is employed to find a solution by minimizing the following objective function:


where u is the upper bound that is chosen to be physically meaningful, η(S) a stabilizing function, α the regularization parameter, W the weight matrix, and VW=VTWV .

To demonstrate the feasibility and efficiency of our new algorithm, we carried out an experiment using a heterogeneous mouse chest phantom. The physical phantom of 30-mm height and 30-mm diameter was fabricated. It consisted of four types of high-density polyethylene materials: (8624K16), nylon 6/6 (8538K23), delrin (8579K21), and polypropylene (8658K11) (McMaster-Carr Supply Company, Chicago, Illinois), to represent muscle (M), lungs (L), heart (H), and bone (B) respectively. Based on the diffuse model of photon propagation, the optical parameters of the four materials at wavelength of about 650nm were independently found from transillumination experimental data,12 which are listed in Table 1 . The luminescent light stick (Glowproducts, Victoria, British Columbia, Canada) was selected as the testing source. The stick consisted of a small glass vial containing one chemical solution and a large plastic vial containing another solution, with the former being embedded in the latter. By bending the plastic vial, the glass vial can be broken to mix the two solutions and emit red light around 650nm , which is spectrally similar to bioluminescent light generated by the luciferase. Two small holes of diameter 0.6mm and height 3mm were longitudinally drilled in the left lung region of the phantom with their centers at (9.0,1.5,15.0) and (9.0,1.5,15.0) , respectively, as shown in Fig. 1 . Two catheter tubes about 1.9mm in height were filled with red luminescent liquid, and were placed inside the two holes as light sources, respectively. Emitting light power of two catheter tubes was measured with the CCD camera. They were 105.1 nW and 97.4nW , respectively.

Fig. 1

Heterogeneous mouse chest phantom. (a) A geometrical model of the phantom embedded with two sources; (b) a middle cross-section through two embedded hollow cylinders for hosting luminescent sources in one lung.


Table 1

Optical parameters of the heterogeneous mouse chest phantom.

MaterialMuscle (M)Lung (L)Heart (H)Bone (B)
μa[mm1] 0.0070.0230.0110.001
μs[mm1] 1.0312.0001.0960.060

The experiment was conducted in a totally dark environment to avoid background noise. The flux density on the cylindrical surface of the phantom was recorded with the sensitive CCD camera, along four radial directions separated by 90deg , as schematically shown in Fig. 1b. During each data acquisition session, one luminescent view was taken by exposing the camera for 60s , as shown in Fig. 2

Fig. 2

Luminescent views of the side surface of the cylindrical phantom taken using a CCD camera in four directions of 90deg apart, (a)–(d) Front, back, left, and right views, respectively.


A permissible source region was assigned as Ωs={(x,y,z)x<0,13.0<z<17.0,(x,y,z)L} to regularize the BLT solution. Along the longitudinal direction high signal-to-noise ratios (SNR) were clustered between z=1.9 to 28.1mm relative to the phantom bottom. Beyond the above region, the SNR were insignificant, and ignored to reduce the computational cost. To simulate the photon propagation in the phantom, a geometrical model of diameter 30mm and height 26.2mm was constructed corresponding to a middle section of the physical phantom. Based on this model, a discrete mesh was generated for the external boundary and internal interface. The external boundary mesh consisted of 960 quadrilateral elements and 1024 measurement datum points. The internal interface mesh consisted of 1680 quadrilateral elements and 1792 nodes. The permissible source region Ωs was also discretized into 308 wedge elements, as shown in Fig. 3b. The measured photon density at each detector location was computed from the CCD luminescent image using our calibration formula. Then, the proposed algorithm was applied to reconstruct the light source distribution in the heterogeneous phantom. The reconstructed results correctly revealed that there were two strong light sources in the phantom located at (7.5,2.7,15.0) with flux density 73.6nWmm3 and at (7.5,2.7,15.0) with 56.1nWmm3 , respectively. The former was estimated to yield a total power of 108.2nW (the total power=source volume×source flux density=1.47mm3×73.6nWmm3=108.2nW) , while the latter was computed to have a total power of 82.5nW (1.47×56.1=82.5nW) . Figures 3a and 3b shows the reconstructed source distribution. The differences between the reconstructed and real source positions were 1.9 and 1.9mm for the two sources, respectively. The relative errors in the source strength were 3.0% and 15.3%, respectively. The computed surface photon density profiles based on the reconstructed light sources were in good agreement with the experimental counterparts with an average relative error about 20%.

Fig. 3

BLT reconstruction of the mouse chest phantom. (a) The 3-D and (b) the 2-D reconstructed source distribution in the permissible region.


In conclusion, we have developed a boundary integral method to reconstruct a 3-D bioluminescent source distribution embedded in a heterogeneous object of interest from measured flux data on the external surface of the object. The feasibility of the new reconstruction method has been demonstrated in a physical heterogeneous phantom experiment. The results have indicated that the method can reliably localize the source locations, and recover the source strength fairly well. The reconstruction results may be further improved by decreasing the boundary element size, lowering the measurement noise, and increasing the accuracy of optical parameters. The method only requires the discretization on the object boundary and structural interface, which is much easier than meshing for volumetric finite elements. Hence, the new method can handle a complex geometrical model much efficiently than finite element based BLT algorithm.


This work is supported by an NIH/NIBIB Grant EB001685 and the U.S. Army's Breast Cancer Research program Concept Award W81XWH-05–1–0461.


1.  V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol.1087-0156 10.1038/nbt1074 23, 313–320 (2005). Google Scholar

2.  C. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng.1523-9829 10.1146/annurev.bioeng.4.111901.093336 4, 235–260 (2002). Google Scholar

3.  W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt.1083-3668 10.1117/1.1413210 6, 432–440 (2001). Google Scholar

4.  G. Wang, E. A. Hoffman, and G. McLennan, “Systems and methods for bioluminescent computed tomographic reconstruction,” Patent disclosure filled in July 2002; US provisional patent application filled in March 2003; US patent application filed in March 2004. Google Scholar

5.  G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescent CT scanner,” Radiology0033-8419 229(P), 566 (2003). Google Scholar

6.  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.0094-2405 10.1118/1.597634 22, 1779–1792 (1995). Google Scholar

7.  J. Ripoll, D. Yessayan, G. Zacharakis, and V. Ntziachristos, “Experimental determination of photon propagation in highly absorbing and scattering media,” J. Opt. Soc. Am. A0740-3232 22, 546–551 (2005). Google Scholar

8.  H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A mouse optical simulation environment (MOSE) to investigate bioluminescent phenomena with the Monte Carlo method,” Acad. Radiol.1076-6332 11, 1029–1038 (2004). Google Scholar

9.  G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys.0094-2405 10.1118/1.1766420 31, 2289–2299 (2004). Google Scholar

10.  W. Cong, D. Kumar, Y. Liu, A. Cong, and G. Wang, “A practical method to determine the light source distribution in bioluminescent imaging,” Proc. SPIE0277-786X 5535, 679–686 (2004). Google Scholar

11.  M. Jiang and G. Wang, “Image reconstruction for bioluminescence tomography,” Proc. SPIE0277-786X 5535, 335–351 (2004). Google Scholar

12.  W. Cong, G. Wang, D. Kumar, Y. Liu, “Practical reconstruction method for bioluminescence tomography,” Opt. Express1094-4087 10.1364/OPEX.13.006756 13, 6756–6771 (2005). Google Scholar

13.  M. Jiang, Y. Li, and G. Wang, “Inverse problems in bioluminescence tomography,” in Frontier and Prospect of Contemporary Applied Mathematics, Ciarlet and Tatsien, Eds., Higher Education Press (Beijing) & World Scientific (2005). Google Scholar

14.  X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express1094-4087 10.1364/OPEX.12.003996 12, 3996–4000 (2004). Google Scholar

15.  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.0031-9155 10.1088/0031-9155/50/17/021 50, 4225–4241 (2005). Google Scholar

16.  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.0031-9155 10.1088/0031-9155/50/23/001 50, 5421–5441 (2005). Google Scholar

17.  A. Cong and G. Wang, “Multi-spectral bioluminescence tomography: methodology and simulation,” Int. J. Biomed. Imag. 1, 73–79 (2006). Google Scholar

18.  T. Troy, O. Coquoz, C. Kuo, D. Zwarg, and B. Rice, “Single-view bioluminescent tomography in small animal models,” Molecular Imaging1535-3508 4, 369–370 (2005). Google Scholar

19.  F. Paris and J. Canas, Boundary Element Method: Fundamentals and Applications, Oxford University Press, New York (l997). Google Scholar

Wenxiang Cong, Ge Wang, "Boundary integral method for bioluminescence tomography," Journal of Biomedical Optics 11(2), 020503 (1 March 2006). http://dx.doi.org/10.1117/1.2191790


Chemical elements




Data modeling



Back to Top