Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Opt Soc Am A Opt Image Sci Vis. Author manuscript; available in PMC 2010 April 13.
Published in final edited form as:
PMCID: PMC2853929

Bioluminescence tomography based on the phase approximation model


A reconstruction method of bioluminescence sources is proposed based on a phase approximation model. Compared with the diffuse approximation, this phase approximation model more correctly predicts bioluminescence photon propagation in biological tissues, so that bioluminescence tomography can accurately locate and quantify the distribution of bioluminescence sources. The compressive sensing (CS) technique is applied to regularize the inverse source reconstruction to enhance numerical stability and efficiency. The numerical simulation and phantom experiments demonstrate the feasibility of the proposed approach.


Bioluminescence is derived from the expression of luciferase enzymes (such as in firefly or Renilla). When cells are encoded with the luciferase enzymes and introduced in a living mouse, they serve as optical probes and allow bioluminescence imaging for biological and preclinical studies [1,2]. Bioluminescent photons transport through biological tissues and can be externally detected using a highly sensitive charge-coupled device (CCD) camera [3,4]. In vivo bioluminescent imaging (BLI) is based on the detection of light emission from the bioluminescence probes in small animals to reveal the molecular and cellular activities sensitively and specifically [1-4]. Compared with fluorescence imaging, in which external excitation light is required to enable fluorescence at the cost of strong tissue autofluorescence, the background noise of BLI is rather low because there is no external illumination source. Therefore, BLI is a low-cost, noninvasive, and important tool for real-time analysis of disease processes at the molecular and cellular levels. This technology has been applied to numerous biological models to detect genetic interactions, track cells, evaluate therapies, and facilitate drug development [1-3].

The two bioluminescent reporters most commonly used are Renilla luciferase (coelenterazine substrate) and firefly luciferase (D-luciferin substrate). Renilla luciferase generally produces blue light, with the peak emission wavelength about 480 nm. Firefly luciferase emits yellow-to-red light, with the peak emission wavelength about 560 nm [1,2]. Bioluminescent-photon propagation in biological tissues is subject to both scattering and absorption and can be well modeled by the radiative transfer equation (RTE) [5,6]. A number of computational schemes have been developed to solve RTE directly, including Monte Carlo simulation techniques [7] and discrete ordinate methods [8,9]. Although the Monte Carlo simulation may provide a highly accurate solution to the RTE, the stochastic pattern and overwhelming computational cost make it inappropriate for bioluminescence tomography (BLT). The discrete ordinate method discretizes the RTE in solid angle directions, and the resultant system of algebraic equations has to be solved iteratively. Hence, it is computationally inefficient in practice. Because of the computational complexity of the RTE, current BLT algorithms are mainly based on the diffusion approximation (DA) model, which is a low-order approximation of the RTE [10-17]. However, the DA model works well only in a highly scattering and weakly absorbing medium [6,18,19]. For quantitative imaging of bioluminescence in the 450–650 nm range, the absorption coefficients are often relatively large in some organs of small animals. Therefore, the DA model is not sufficiently accurate in these cases. The model mismatch can significantly compromise reconstruction quality in this ill-posed inversion source scenario. The phase approximation (PA) model is more accurate than the DA model to describe photon propagation in biological tissues over a broad range of optical parameters [20,21]. In contrast to the DA derived from the first-order approximation of the photon radiance, the PA is based on the exact solution of the photon radiance assuming generalized Delta–Eddington phase functions. The PA model relaxes the limitation of the scattering domination required by the DA and performs well in simulating bioluminescent-photon propagation.

BLT localizes and quantifies bioluminescence probes inside a small animal from externally measured bioluminescence signals. Since the bioluminescence signal is measured only on a small-animal body surface, the reconstruction of 3D bioluminescence sources in a mouse is a typical underdetermined problem. Compressive sensing (CS) theory asserts that one can reconstruct images with a sparse representation from far fewer samples or measurements than what the Nyquist sampling theorem demands [22-24]. The success of CS relies on both the sparsity of an underlying image and the incoherence of the sensing mode. This premise radically changes the solution scheme for an underdetermined system, and such a solution can be achieved with good accuracy by a convex optimization procedure. The sparseness characteristic of a bioluminescent source distribution suggests that the CS-type approach would be a promising method for BLT. Lu et al. applied the CS technique for the reconstruction of bioluminescence sources, which used an l1 -norm regularization-optimization method based on the DA model [25]. Its feasibility was validated in both numerical simulation and phantom experiments.

In this paper, a reconstruction method of bioluminescence sources is proposed based on a PA model. In Section 2, the PA model is briefly described. In Section 3, the CS technique is applied to regularize the bioluminescence source reconstruction to enhance numerical stability. In Section 4, numerical simulation and phantom experiments are described to demonstrate the feasibility of the proposed method. In Section 5, relevant issues are discussed and conclusions drawn.


Bioluminescence photon propagation in biological tissues is modeled by a steady-state RTE [5,6],

v[center dot][nabla]L(r,v)+(μa+μs)L(r,v)=μs4πL(r,v)p(v,v)dv+14πS(r),r[set membership]Ω,

where Ω is an object support, L(r,v) a photon radiance at a location r in a unit direction v (Watts·mm−2·sr−1), S(r) an isotropic bioluminescence source (Watts·mm−3), μs a scattering coefficient (mm−1), and μa an absorption coefficient (mm−1) at the wavelength λ. A scattering phase function p(v,v′) gives the probability for a photon coming in a direction v′ and scattered into a direction v. Biological tissue scatters photons strongly in the forward direction. Hence, the phase function can be well modeled by a generalized Delta–Eddington phase function [6,20,21,26]:

p(v[center dot]v)=14π(1f)+12πfδ(1v[center dot]v),

where f[set membership][−1,+1] is the weight factor measuring the anisotropy of the photon scattering, which is referred to as the anisotropy weight. The phase function is a linear combination of the isotropic scattering and the peaked forward scattering with the anisotropic weight f. The original Delta–Eddington phase function defines the parameter f as a fixed value g, which is the mean value of the cosine of the scattering angles. The generalized Delta–Eddington function defines the anisotropy weight as a variable, which is related to the photon absorbing and scattering properties of the media. The anisotropy weight can be determined by Monte Carlo simulation with known absorption and scattering coefficients and a conventional anisotropic factor g. Alternatively, the anisotropy weight can be determined along with the absorption and scattering coefficients by diffuse optical tomography techniques [27,28]

Using the generalized Delta–Eddington phase function, RTE can be simplified to the following integral equation with respect to photon fluence rate Φ(r) [20,21],

Φ(r)=ΩG(r,r)[μs(r)Φ(r)+S(r)]drrd1+rd[contour integral operator][partial differential]ΩG(r,r)Φ(r)β[center dot]ndr,

where the integral kernel G(r,r)=exp(|rr|01μtr[(1t)r+tr]dt)/4π|rr|2, unit vector β =(rr′)/|rr′|, and n is the outward unit normal at r on [partial differential]Ω. The internal reflection coefficient rd can be calculated by rd =−1.4399η−2+0.7099η−1+0.6681+0.0636η, η =ntissue/nair, ntissue is the refractive index in biological tissues, and nair the refractive index of the surrounding air. Equation (3) with respect to the photon fluence rate Φ(r) is a well-posed integral equation of the second kind [29]. For simplicity, we call it the phase approximation (PA) model because it is derived from an approximate phase function.


Based on the structural/anatomical information on a living animal and known optical parameters μa, μs, f,or g in every sub-region of Ω, the BLT problem is reduced to determine a bioluminescence source distribution according to Eq. (3) [12-14]. Let Ω be partitioned into finite elements of N vertex nodes and the photon fluence rate be approximated by nodal-based basis functions,


Substituting Eq. (4) into Eq. (3), we obtain the following matrix equation,


where M, B, and F represent the corresponding discrete integral kernels in Eq. (3), with the components of the matrices defined by

mi,j=ΩG(ri,r)μs(r)ϕj(r)dr,bi,j=rd(1+rd)[contour integral operator][partial differential]ΩG(ri,r)ϕj(r)(β[center dot]n)dr,fi=ΩG(ri,r)ϕj(r)dr,

{S} denotes the source distribution in a permissible source region, and {Φ} consists of photon fluence rate values at the nodes in Ω. Furthermore, {Φ} can be divided into {Φin} at the internal nodes and {Φex} at the detectable boundary nodes on [partial differential]Ω. By Eq. (3), the matrix IMB is diagonally dominant and invertible. Hence, the photon fluence rate {Φ} can be solved from Eq. (5),


where K=(IMB)−1F. Removing those rows of K that correspond to {Φin}, K becomes . Accordingly, the following linear equation system is obtained to link unknown source density {S} and measurable photon fluence rate {Φex},


Theoretically, measured photon fluence rate is two-dimensional, a bioluminescent source distribution is three-dimensional, and the BLT is underdetermined. An effective way to reduce the number of unknown variables is to impose permissible source regions [13]. Biolumines-cent views taken by a CCD camera indicate high-value clusters, and prominent high-density clusters should be closer to sources than other surface places. In some cases, permissible source regions can be also outlined based on prior knowledge. With these hints it is easy to generate permissible source regions where a bioluminescence source may exist. Additionally, an adaptive procedure can also efficiently determine permissible source regions: an initial permissible source region can be given to perform the reconstruction of the BLT on a low-resolution scale, and high values of reconstructed sources will be very likely located in a neighborhood of or within real light sources. Then, the permissible source regions can be updated in reference to these highly valued nodes. Applying the permissible source region can significantly reduce the number of unknown variables and enhance the stability of the BLT reconstruction.

Many in vivo studies involve isolated bioluminescence sources, such as in vivo tumor or metastasized tumor models, and the resultant bioluminescence source distribution would have a sparse structure. Because of the highly scattering nature of bioluminescence light in biological tissues and the measurement data corrupted with noise, multi-pseudo sources of smaller powers surrounding a true source may be equivalent to the true source from the external surface observation. To exclude these pseudo sources, the BLT model can be regularized to minimize the number of light sources subject to Eq. (8), that is, performing the following optimization:


where {·{0 is the l0 norm, i.e., the number of nonzero entries of a given vector. This is a difficult combinatorial problem and has exponential complexity. Based on the compressive sensing theory, Eq. (9) is equivalent to the following l1-norm optimization [30]:


With corrupted measurement data, Eq. (10) can be translated into the following optimization by relaxing the data fidelity term,


where ε is the acceptable reconstruction error, determined by the noise level. Equation (11) is a convex programming problem and can be solved in polynomial time. Hence, the BLT problem is now put to a CS-type framework [22-24], allowing more reliable and robust reconstruction than the least-squares reconstruction method.


A. Numerical Simulation

With the Monte Carlo method, we computed the anisotropy weight from the optical parameters: absorption coefficient μa=0.20 mm−1, scattering coefficient μs =14.5 mm−1, anisotropic coefficient g=0.9, and relative refractive index of 1.37. As a result, the anisotropy weight f was estimated as 0.925 via fitting to the photon fluence rates based on a phantom [21]. A spherical phantom of radius 10 mm was employed in the numerical simulation. The phantom was discretized into 65,775 tetrahedral elements and 12,044 nodes. A total of 2108 virtual detectors were distributed over the phantom surface to record the photon fluence rates. Then two spherical light sources of radius 0.6 mm were embedded into the phantom. The centers of the two light sources were at (3.5, 0.0, 0.0) and (2.5, 2.5, 0.0), respectively. Each source was assigned a power of 10 nW. The photon fluence rate data at the detectors were computed in the Monte Carlo simulation with the Henyey-Greenstein phase function to mimic real measurement [7]. A permissible source region Ωs ={ (x,y,z)|0≤x≤5,−1≤y≤4,−1.5≤z≤1.5} was assigned for the BLT reconstruction. The PA-based BLT reconstruction was performed using the l1-norm optimization to identify the light source distribution from the generated data. It was found that the localization errors of the light sources reconstructed with the PA model were less than 0.15 mm, and the relative errors in terms of source power were about 5%, as shown in Fig. 1(a). The reconstructed source location was computed as the average of nodal source positions weighted by source densities, and reconstructed intensity was defined as the sum of the nodal source intensities after thresholding. For comparison, the DA-based BLT reconstruction was also performed using the l1-norm optimization and the least-squares optimization, respectively, in the identical setting. It was found that the localization errors of the light sources reconstructed with the DA model exceeded 2.5 mm for both optimization criteria, and the relative errors in terms of source power were about 25% with the least-squares optimization technique and 20% with the l1-norm optimization technique. These results show that the proposed PA-based reconstruction is superior to the DA-based reconstruction, as shown in Fig. 1(b) and 1(c).

Fig. 1
(Color online) Comparison of PA-based and DA-based source reconstructions in the numerical simulation. (a) PA-based reconstruction, (b) DA-based reconstruction with the least-squares optimization method, and (c) DA-based reconstruction with the CS-based ...

B. Phantom experiments

A physical cylindrical phantom of 15 mm radius and 30 mm height was fabricated from high-density nylon to mimic a section of a mouse body, as shown in Fig. 2(a). The optical parameters of the phantom were determined using diffuse optical tomography as μa=0.0172 mm−1, μs=5.393 mm−1, and f=0.933 at a wavelength of 550 nm. The luminescent liquid (GlowProducts, Victoria, British Columbia) had a spectral range [300,800] nm. This liquid was used to fill in a small hole of diameter 0.58 mm and height 2 mm centered at (−5.0, 0.0, 15.0) in the phantom as a testing light source, as shown in Fig. 2(b). The luminescent source was directly measured using a bandpass filter and a CCD camera (Princeton Instruments VersArray: 1300B, Roper Scientific Inc, Trenton, New Jersey), and it was found that its total power was 0.4335 nW at a wavelength of 550 nm. Then the optical imaging experiment was conducted with the phantom in a totally dark environment. The exiting photon on the cylindrical surface of the phantom was filtered through a bandpass filter centered at 550 nm and recorded using the CCD camera along four radial directions 90° apart, as shown in Figs. 3(a)–3(d). The cylindrical phantom was discretized into 58,303 linear tetrahedral elements with 10,801 nodes. A total of 2170 detector positions were defined on the side surface of the phantom. The detector readings were computed from these four CCD views. The permissible source region was set as Ωs={ (x,y,z)|x≤−2,−2≤y≤2,13≤z≤17,x2+y2+z2≤152}. Then the proposed PA-based reconstruction method was performed to identify light source distribution from the measurement data. As a result, the reconstructed source was much closer to the true source: the difference between the reconstructed and the real source positions was only about 1.0 mm. The relative error of the reconstructed source power was about 17.3%, as shown in Figure 4(a). Furthermore, the DA-based reconstruction was also performed using the least-squares optimization and CS-based optimization techniques. The localization error of the reconstructed light source exceeded 2.5 mm for both methods, and the error in term of source power was about 30% and 25% for the least-squares optimization and the CS-based optimization, respectively. The reconstructed sources with the least-squares optimization were dispersed around the true source, as shown in Fig. 4(b), having introduced the pseudo sources. Compared with the two optimization methods based on the DA model, the CS-inspired reconstruction technique was effective in avoiding pseudo sources, as shown in Fig. 4(b) and 4(c).

Fig. 2
Physical phantom. (a) A photograph of the physical phantom, and (b) the cross-section through the luminescence source center, along with the four viewing directions of the CCD camera during the data acquisition process.
Fig. 3
Four luminescent views of the physical phantom taken by a CCD camera in directions 90° apart. (a)–(d) Front, back, left, and right views, respectively, recorded at wavelength 550 nm.
Fig. 4
(Color online) Comparison of PA-based and DA-based source reconstructions in the phantom experiment. (a) The transverse section through the source reconstructed with the PA model, (b) the transverse section through the source reconstructed with the DA ...


BLT is an underdetermined linear inverse source problem, leading to multiple solutions and aberrant source reconstructions. Compressive sensing is a powerful tool for solving the underdetermined problem, which allows that images of sparse structures be exactly recovered with overwhelming probability from far fewer data than what the Nyquist theorem requires. Thus, this technique makes the bioluminescent source reconstruction more accurate and stable than the least-squares optimization method when bioluminescence source distribution consists of isolated clusters. Theoretically, compressive sensing requires an isometry property of the matrix in Eq. (8), which is typically called the restricted isometry property (RIP). This property requires a bounded condition number for all sub-matrices. Although verifying the RIP condition is difficult, the goal of bioluminescence source reconstruction is essentially the same as that of compressive sensing tasks, that is, to find sparse solutions from underdetermined systems of linear equations. Therefore, compressive sensing techniques do help enhance the numerical stability of the BLT reconstruction.

Performing in vivo experiments with BLT involves modeling of the individualized anatomy and determination of the optical parameters of biological tissues. Since the mouse anatomy is rather complex, it is difficult to include all the features in a geometrical model. In practice, a mouse volume should be segmented into major organ regions to simplify the modeling. The optical parameters are considered piecewise constant and can be uniquely determined by diffuse optical tomography techniques [28]. We are currently studying a new modality-fusion methodology. Specifically, we would like to solve these problems utilizing photoacoustic tomography techniques, which constitute a noninvasive imaging method for high-resolution mapping of optical absorption in tissues [31]. Detailed in vivo studies with PA-based BLT will be reported in another paper. Inhomogeneous optical properties of tissues will be reflected in the variation of the weighting matrix in Eq. (8) without any technical difficulty, and the PA-based BLT approach can be applied to complex biological tissues for BLT.

In summary, we have developed a compressive sensing-inspired BLT approach based on the PA model. It has been demonstrated in the numerical simulation and the phantom experiments that our approach can accurately localize and quantify a bioluminescence source distribution from the surface measurement of bioluminescence signals.


This work is supported by a NIH/NCI grant CA135151 and NIH/NIBIB grants EB008476, EB001685, and EB006036.


1. Contag CH, Bachmann MH. Advances in bioluminescence imaging of gene expression. Annu. Rev. Biomed. Eng. 2002;4:235–260. [PubMed]
2. Ray P, Wu AM, Gambhir SS. Optical bioluminescence and positron emission tomography imaging of a novel fusion reporter gene in tumor xenografts of living mice. Cancer Res. 2003;63:1160–1165. [PubMed]
3. Ntziachristos V, Ripoll J, Wang LV, Weissleder R. Looking and listening to light: the evolution of whole-body photonic imaging. Nat. Biotechnol. 2005;23:313–320. [PubMed]
4. Rice W, Cable MD, Nelson MB. In vivo imaging of light-emitting probes. J. Biomed. Opt. 2001;6:432–440. [PubMed]
5. Ishimaru A. Wave Propagation and Scattering in Random Media. Vol. 1. Academic; 1978.
6. Welch AJ, Van Gemert MJC. Optical-Thermal Response of Laser-Irradiated Tissue. Plenum; 1995.
7. Wang LH, Jacques SL, Zheng LQ. MCML—Monte Carlo modeling of photon transport in multi-layered tissues. Comput. Methods Programs Biomed. 1995;47:131–146. [PubMed]
8. Abdoulaev Gassan S., Hielscher Andreas H. Three-dimensional optical tomography with the equation of radiative transfer. J. Electron. Imaging. 2003;12:594–601.
9. Klose AD, Ntziachristos V, Hielscher AH. The inverse source problem based on the radiative transfer equation in optical molecular imaging. J. Comput. Phys. 2005;202:323–345.
10. Wang G, Hoffman EA, McLennan G. Systems and methods for bioluminescent computed tomographic reconstruction. Patent disclosure filed in July 2002; U.S. provisional patent application filed in March 2003; U.S. patent application filed in March 2004.
11. Wang G, Hoffman EA, McLennan G, Wang LV, Suter M, Meinel J. Development of the first bioluminescent CT scanner. Radiology. 2003;229:566.
12. Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. Med. Phys. 2004;31:2289–2299. [PubMed]
13. Cong W, Wang G, Kumar D, Liu Y, Jiang M, Wang LV, Hoffman EA, McLennan G, McCray PB, Zabner J, Cong A. Practical reconstruction method for bioluminescence tomography. Opt. Express. 2005;13:6756–6771. [PubMed]
14. Alexandrakis G, Rannou FR, Chatziioannou AF. Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study. Phys. Med. Biol. 2005;50:4225–4241. [PMC free article] [PubMed]
15. Dehghani H, Davis SC, Jiang S, Pogue BW, Paulsen KD, Patterson MS. Spectrally resolved bioluminescence optical tomography. Opt. Lett. 2006;31:365–367. [PubMed]
16. Chaudhari AJ, Darvas F, Bading JR, Moats RA, Conti PS, Smith DJ, Cherry SR, Leahy RM. Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging. Phys. Med. Biol. 2005;50:5421–5441. [PubMed]
17. Kuo C, Coquoz O, Troy TL, Xu H, Rice BW. Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging. J. Biomed. Opt. 2007;12:24007. [PubMed]
18. Flock ST, Patterson MS, Wilson BC, Wyman DR. Monte Carlo modeling of light propagation in highly scattering tissues-I: model predictions and comparison with diffusion theory. IEEE Trans. Biomed. Eng. 1989;36:1162–1168. [PubMed]
19. Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: Boundary and source conditions. Med. Phys. 1995;22:1779–1792. [PubMed]
20. Cong W, Shen H, Cong A, Wang G. Integral equation of the photon fluence rate and flux based on a generalized Delta-Eddington phase function. J. Biomed. Opt. 2008;13:024016. [PMC free article] [PubMed]
21. Cong W, Shen H, Cong A, Wang Y, Wang G. Modeling photon propagation in biological tissues using a generalized Delta-Eddington phase function. Phys. Rev. E. 2007;76:051913. [PubMed]
22. Candès EJ, Tao T. Decoding by linear programming. IEEE Trans. Inf. Theory. 2005;51:4203–4215.
23. Candes EJ, Romberg J, Tao T. Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information. IEEE Trans. Inf. Theory. 2006;52:489–509.
24. Candes EJ, Romberg J, Tao T. Stable signal recovery from incomplete and inaccurate measurements. Commun. Pure Appl. Math. 2006;59:1207–1223.
25. Lu Y, Zhang X, Douraghy A, Stout D, Tian J, Chan TF, Chatziioannou AF. Source reconstruction for spectrally-resolved bioluminescence omography with sparse a priori information. Opt. Express. 2009;17:8062–8080. [PMC free article] [PubMed]
26. Joseph JH, Wiscombe WJ, Weinman JA. The Delta-Eddington approximation for radiative flux transfer. J. Atmos. Sci. 1976;33:2452–2459.
27. Arridge SR. Optical tomography in medical, imaging. Inverse Probl. 1999;15:R41–R93.
28. Harrach B. On uniqueness in diffuse optical tomography. Inverse Probl. 2009;25:055010.
29. Atkinson KE. Numerical Solution of Integral Equations of the Second Kind. Cambridge Univ. Press; 1997.
30. Donoho DL, Tanner J. Sparse nonnegative solution of underdetermined linear equations by linear programming. Proc. Natl. Acad. Sci. U.S.A. 2005;102:9446–9451. [PubMed]
31. Li C, Wang LV. Photoacoustic tomography and sensing in biomedicine. Phys. Med. Biol. 2009;54:R59–R97. [PMC free article] [PubMed]