## Abstract

Bioluminescence tomography (BLT) poses a typical ill-posed inverse problem with a large number of unknowns and a relatively limited number of boundary measurements. It is indispensable to incorporate *a priori* information into the inverse problem formulation in order to obtain viable solutions. In the paper, Bayesian approach has been firstly suggested to incorporate multiple types of *a priori* information for BLT reconstruction. Meanwhile, a generalized adaptive Gaussian Markov random field (GAGMRF) prior model for unknown source density estimation is developed to further reduce the ill-posedness of BLT on the basis of finite element analysis. Then the distribution of bioluminescent source can be acquired by maximizing the log posterior probability with respect to a noise parameter and the unknown source density. Furthermore, the use of finite element method makes the algorithm appropriate for complex heterogeneous phantom. The algorithm was validated by numerical simulation of a 3-D micro-CT mouse atlas and physical phantom experiment. The reconstructed results suggest that we are able to achieve high computational efficiency and accurate localization of bioluminescent source.

© 2009 Optical Society of America

## 1. Introduction

Molecular imaging, especially small-animal molecular imaging, is rapidly becoming an active research filed [1]. Among all the imaging modalities that are used for molecular imaging, optical techniques stand out because of their ability to reveal molecular and cellular activities in a small animal *in vivo* [1, 2, 3, 4], including bioluminescence tomography (BLT) [5] and fluorescence tomography (FMT) [6], and so on. Their mechanism is to identity light source from the collected emitted photons from multiple 3-D directions with respect to a living mouse marked by reporter luciferases, which has been applied to cancer research including studies of tumor burden, response to therapy, assessment of gene expression, and development of metastatic lesions [7]. In this paper, we focus exclusively on bioluminescence tomography.

There has been a great effort lately devoted to transforming bioluminescence imaging from a 2-D, planar bioluminescent imaging technique into a truly 3-D tomographic imaging modality applicable to small animals, because planar bioluminescent imaging cannot generate depth information [3, 8]. Therefore, bioluminescence tomography (BLT) has become a research hot spot. Compared with other imaging modalities, the advantages of BLT are sensitivity, speed, non-invasiveness, low cost and background noise [9, 10]. The aim of BLT is to reconstruct an internal bioluminescent source distribution based on both the outgoing bioluminescent photons captured using a highly sensitive charge-coupled device (CCD) camera and a prescanned tomographic volume, such as a CT/Micro-CT or MRI volume, of the same mouse [11].

Bioluminescent photons propagation in biological tissue is governed by the radiative transfer equation (RTE) [12]. However, to solve for the distribution of bioluminescent source in a large 3-D volume, as is required in BLT, the RTE is computationally expensive [13]. Given the dominance of scattering over absorption in bioluminescent imaging, RTE can be replaced by diffusion equation and Robin boundary condition, which provide a quite accurate description of propagation of the photon transport. However, a major difficulty in determining the bioluminescent source distribution is imposed by multiple scattering of photons that propagate through heterogeneous biological tissues. In addition, the absence of external illumination accords a high sensitive signal, it complicates the tomographic problem. So far, although BLT has made a great progress in modeling and reconstruction algorithms [9, 11, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], the current technique has not full explored the potential of this approach.

Wang *et al*. [11] has theoretically proved the solution uniqueness for BLT under practical constrains and presented that the BLT problem can be reduced to an inverse source problem by incorporating sufficient *a priori* information. The commonly used *a priori* information is the permissible source region strategy and its importance in BLT reconstruction has been well recognized [9]. Now, various reconstruction methods and real experiments have been developed and reported based on the strategy [10, 18, 19, 21, 22, 24] and promising results have demonstrated the merits and potential of the strategy. Aside from the permissible source region strategy, the anatomic structure information of the small animal which can be imaged by X-ray CT and MRI techniques can be considered as *a priori* information. The geometrical shape of major organ regions can be acquired by 3-D computer graphic techniques. Then, the organ region can be associated with tissue optical parameters which can be determined by diffusion optical tomography (DOT). Based on the permissible source region strategy, anatomic structure information and tissue optical parameters, the ill-posedness of BLT will be further reduced and the viable solutions can be solved [25].

Bayesian approach provides a nature framework to incorporate multiple types of *a priori* information. Recently, Bayesian approachs have been applied to nonlinear inverse problems such as diffusion optical tomography (DOT) [26, 27]. But surprisingly, Bayesian approach has not been applied to the BLT problem. In addition, Bayesian reconstruction methods for DOT can not be directly applied to BLT problem because of the differences between DOT problem and BLT problem. Furthermore, the known Bayesian methods used in DOT problem are not appropriate to complex phantoms such as real mouse [27]. Therefore, there is a critical need to derive and construct a method based on Bayesian approach that can accurately reconstruct bioluminescent sources in heterogeneous and complex tissues.

In the work, various *a priori* information is incorporated into a Bayesian network to reduce the ill-posedness of BLT problem and an adaptive Gaussian Markov random field prior model for unknown source density estimation is developed on the basis of finite element theory. Under the framework of Bayesian approach, we develop a reconstruction algorithm to identify the bioluminescent source distribution in complex phantoms. The paper is organized as follows. The next section describes the details of the proposed reconstruction algorithm. In the third section, numerical experiments are presented to demonstrate the feasibility and merits of the algorithm with a 3-D micro-CT mouse atlas. With the model, the reconstructions were carried out from the data on the surface of the phantom, which were synthesized by molecular optical simulation environment (MOSE) developed by our group [28]. Furthermore, real physical phantom experiment was performed to valuate the algorithm. Finally, Section 4 discusses the relevant issues and makes concluding remarks.

## 2. Methods

#### 2.1. Forward model for BLT

In bioluminescence imaging experiment, it is necessary to depict the propagation of the photon transport accurately. Generally, the bioluminescent photon is subject to both scattering and absorption in biomedical tissues. However, in the wavelength range of bioluminescent light, the biomedical tissues are highly scattering, photon scattering predominates over photon absorption. Assuming the bioluminescent source density is stable when photons are collected, the photon propagation in biological tissues can be modeled by steady-state diffusion equation and Robin boundary condition [9, 29, 30]:

where **r**∊ℝ^{3} represents the location vector, Ω and *∂*Ω are the domain and its boundary respectively; Φ(**r**) denotes the photon flux density [*Watts/mm*
^{2}]; **x**(**r**) is the source energy density [*Watts/mm*
^{3}]; *µ _{a}*(

**r**) is the absorption coefficient [

*mm*

^{-1}];

*D*(

**r**)=1/(3(

*µ*(

_{a}**r**)+(1-

*g*)

*µ*(

_{s}**r**))) is the optical diffusion coefficient [

*mm*],

*µ*(

_{s}**r**) the scattering coefficient [

*mm*

^{-1}] and

*g*the anisotropy parameter;

*ν*(

**r**) is the unit outer normal on

*∂*Ω. Given the mismatch between the refractive indices

*n*for Ω and

*n*′ for the external medium,

*A*(

**r**;

*n,n*′) can be approximately represented:

where *n*′ is close to 1.0 when the mouse is in air; *R*(**r**) can be approximated by: *R*(**r**)≈-1.4399*n*
^{-2}+0.7099*n*
^{-1}+0.6681+0.0636*n* [31]. The measured quantity is the outgoing flux density *Q*(**r**) on *∂*Ω, that is:

#### 2.2. Formulation of the Bayesian Problem

We use the vector x to denote the set of unknown source density. Φ* ^{meas}_{k}* is assumed as the measured photon flux density at the

*kth*detector position (

*k*=1,2, …,

*M*) and then organize the measurements as a single column vector

**y**:

**y**=[

**Φ**

^{meas}

_{1},

**Φ**

^{meas}

_{2}, …,

**Φ**

*]*

^{meas}_{M}^{T}. Bayesian methods provide a natural framework for incorporating prior information about unknown source density

**x**. The MAP estimate of x given the measurement vector y can be represented as follows:

$$=\underset{\mathbf{x}\ge 0}{\mathrm{arg}\phantom{\rule[-0ex]{.2em}{0ex}}max}\left\{\mathrm{log}\rho (\mathbf{y}\mid \mathbf{x})+\mathrm{log}\rho \left(\mathbf{x}\right)\right\}$$

Considering the real physical meaning, **x**≥0 nonnegative constrain is adopted. Taking into the anatomical tissue information **C** which is derived from *a priori* anatomical image, the MAP estimate can be modified as

$$=\underset{\mathbf{x}\ge 0}{\mathrm{arg}\phantom{\rule[-0ex]{.2em}{0ex}}max}\left\{\mathrm{log}\rho (\mathbf{y}\mid \mathbf{x},\mathbf{C})+\mathrm{log}\rho (\mathbf{x}\mid \mathbf{C})\right\}$$

When the permissible source region information **PS** is used as *a priori* information to reduce the ill-posedness of BLT, the MAP estimate can be further modified as

$$=\underset{\mathbf{x}\ge 0}{\mathrm{arg}\phantom{\rule[-0ex]{.2em}{0ex}}max}\{\mathrm{log}\rho (\mathbf{y}\mid \mathbf{x},\mathbf{C},\mathbf{PS})+\mathrm{log}\rho (\mathbf{x}\mid \mathbf{C},\mathbf{PS})\}$$

Where *p*(**y**|**x, C, PS**) is the data likelihood and *p*(**x**|**C, PS**) is the conditional probability density function of **x** given the tissue information **C** and permissible source region information **PS**. Given the forward model Eqs. (1) and (2), the data likelihood is governed mainly by the noise statistics, therefore Eq. (7) reduces to

#### 2.3. The data likelihood model

In a Bayesian framework, the data likelihood *p*(**y**|**x**) is required. The bioluminescence experiment is generally operated at low temperature, photon detection can be modeled using shot noise statistics, then the data likelihood can be given by [26]:

Where *α* is the parameter related to the noise variance, Λ is the diagonal covariance matrix, ‖ω‖^{2}
_{Λ}=*ω ^{T}*Λ

*ω*, and the vector value function

*f*(

**x**) represents the exact value of the outgoing flux for the assumed value of the source density

**x**. For BLT, we can assume that the measurements are statistically independent with the variance of each measurement equal to its mean [27], so Λ is diagonal. In our simulations, we approximately assume Λ as:

To compute the likelihood of the data, we need to solve the diffusion equation. The finite element method (FEM) has been used for the purpose and the detailed formulation of *f* (**x**) can be determined. The FEM method used in the paper is proposed in [19]. In finite element analysis, the given domain Ω can be discretized into *N _{T}* tetrahedron elements and

*N*vertex nodes. Taking the permissible source region information

_{v}**PS**into account, there are

*N*tetrahedron elements in the permissible source region which represent the possible unknown source distribution. Then the column vector

_{p}**x**is reduced to the set of

*N*elements’ source density distribution.

_{p}#### 2.4. Source reconstruction

Bayes requires that we must assign a prior to unknown variable **x**. When the generalized adaptive Gaussian Markov random field (GAGMRF) prior model is used [32, 33],

Where *σ* is a normalization hyperparameter and 1≤*p*≤2 with *p*=2 corresponding to the Gaussian case. Given the *ith* and *jth* tetrahedron elements, their four local vertexes are *E _{ie}* and

*E*(

_{je}*e*=1,2,3,4), respectively. If one of vertex of

*E*is identified to any vertex of

_{ie}*jth*element, we assume that

*jth*element is the adjacent element.

*b*denotes the weighting assigned to be inversely proportional to the pair {

_{ij}*i, j*}, for each

*i*, the weights

*b*sum to 1.

_{ij}*N*is adaptively determined in the algorithm, which consists of all pairs of adjacent element.

*z*(

*p*) is partition function. If

*α*is unknown, take the Eqs. (9) and (11), the source reconstruction problem can be stated as following optimization problem:

In the optimization process, *α* is adaptively estimated, which can be solved by viewing the Eq. (12) as a cost function of *α*, and setting the derivative with respect to *α* equal to zero. Then we can obtain:

By substituting Eq. (13) into Eq. (12), the optimization problem is converted into

where **x**̂ is an estimate of the unknown source density **x**. After neglecting constant terms, the Eq. (14) can simplified as :

In practical calculation procedure, for maximizing *l*(**x**) by maximizing with respect to *α* and **x** using the following equations:

Eq. (16) is a straight-forward computation, and Eq. (17) is a computationally expensive optimization problem. In order to calculate Eq. (17), the spectral projected gradient-based large-scale optimization algorithm is employed [26, 27, 34, 35, 36, 37]. In the optimization procedure, each tetrahedron element of the phantom is updated sequentially. After every tetrahedron element has been updated, the procedure is repeated, starting from the first tetrahedron again. We define a single update of every tetrahedron as a scan. There are a number of scans in the optimization procedure until the convergence criterion is satisfied [26].

As far as convergence criterion is concerned, we use the discrepancy between the measured and computational boundary nodal flux data, maximum number of scan, or the discrepancy between the current scan and last scan log posterior probability to evaluate if the procedure should be terminated, that is, ‖**Φ**
^{m}-**Φ**
^{c}‖<*ε*
_{Φ}, *k*≥*K _{max}* or ‖

*logk*-

*log*

_{(k-1)}‖<

*ε*.

_{log}## 3. Experiments and results

A series of experiments were designed to evaluate our proposed algorithm, including numerical and physical phantom experiments. Reconstructions were carried out on a personal computer with 2.8 GHz Pentium 4 processor and 1.75 GB RAM.

#### 3.1. Numerical experiments

### 3.1.1. Phantom and Synthetic data

In order to optimally establish and evaluate the proposed algorithm, a 3-D mouse atlas of micro-CT was employed to provide 3-D anatomical information. For this, we first prepared a 20 g male BALB/c mouse because of its broad applications in biological and medical research and convenient tail vein injection in albino mouse. All animal studies were performed according to procedures approved by the Animal Care and Use Committee. After fasting for 24 hours, the mouse tail was immersed in warm water for 30–60 seconds to increase blood flow to the tail and dilate the vessels prior to injection of the contrast agent. Then the Fenestra LC at a dose of 15 mL/kg was slowly injected over a period of 30–60 seconds via the lateral tail vein. Finally the mouse was anesthetized with intraperitoneal injection of a 13% aqueous urethane (0.15 mL/20 g body weight).

The anesthetized mouse was scanned 40 minutes after injection of the contrast agent using a micro-CT system. During scanning, the X-ray source was operated in a continuous mode with a 1 mm aluminum filter, and the typical tube voltage was 50 kVp with tube current 1.4 mA [43]. The source to detector distance (SDD) and source to object distance (SOD) were set to 541.36 mm and 448.50 mm, respectively. 500 projection views with image size of 2240×2344 pixels were obtained by 360° full scan. Finally, the FDK filtered backprojection method using GPU hardware acceleration was performed for 3-D reconstruction [44].

After micro-CT scanning and 3-D reconstruction, the real mouse was manually segmented into different tissue organs, shown in Fig. 1. In the following bioluminescence experiments, only the thorax micro-CT images were used, including lung, bone, heart, liver and muscle as shown in Fig. 2(a). Optical parameters from the literature [15, 19] and the references therein were assigned to each of the five components, as listed in Table 1. For finite-element-based BLT reconstruction, the mouse phantom was discretized into a tetrahedral-element mesh, illustrated in Fig. 2(b).

For inverse problems, numerical simulations of reconstruction methods usually make use of the synthetic data from the forward problem. Therefore, the *inverse crime* is likely to occur when closely correlated computational ingredients are used in the forward solver and the inverse scheme [19]. In order to avoid the *inverse crime* problem, the synthetic data was produced by Molecular Optical Simulation Environment (MOSE) which is developed based on Monte Carlo method [28]. In the single source simulation, the solid spherical source with 1*mm* radius and source density of 0.238*nano*-*Watts/mm*
^{3} was centered at (22.8,28.6,12.5) inside the lung. The bioluminescent source was sampled by 1.0×10^{6} photons and was assumed to obey the uniform distribution. In MOSE, the 3-D atlas was discritized by triangular elements with 49992 triangles and 24998 surface measurement points, shown in Fig. 2(c). In the following experiments, the default value σ and p used in prior model were set 0.1 and 1.1, respectively. The stopping thresholds *ε*
_{Φ}, *K _{max}* and

*ε*were 1.0×10

_{log}^{-8}, 10 and 1.0, respectively.

### 3.1.2. Experiment results

In the experiment of single bioluminescent source, *a priori* permissible source region can be inferred from the surface light power distribution. Figure 3 shows four views of the surface light power distribution with an angular increment of 90 degrees, and red lines denotes the isoline of light power distribution on the phantom surface. From the Fig. 3, the permissible source region can be inferred as:

*PS*={(*x,y, z*)|22<*y*<30,7<*z*<15, (*x,y, z*)∊*lung*}

when default value σ and *p* were used, the BLT reconstruction was performed and the procedure was terminated after 3 scans, the computational time was about 4 minutes. The corresponding results are shown in Fig. 4 and the reconstruction result indicates that the location of light source is accurately identified. The reconstruction central position and maximum source density were (22.88,28.93,12.87) and 0.234*nano*-*Watts/mm*
^{3}, respectively. The relative error (RE) between the real source and reconstruction result was 1.7%, which was calculated according to *RE*=‖*S _{recons}*-

*S*‖/

_{real}*S*.

_{real}*S*and

_{recons}*S*are reconstructed source density and real source density, respectively.

_{real}In Table 2, the actual bioluminescent source and the bioluminescent source reconstructed with and without Bayesian approach are compiled together for easy comparison. From the Table 2, it can be seen that when the Bayesian approach is used, the bioluminescent source is more accurately reconstructed in terms of the position and source density.

In order to evaluate the robustness of the algorithm, different values of *σ* and *p* were set to perform BLT reconstruction when Gaussian noise with different levels was added to the

synthetic data. The representative results were demonstrated in Fig. 5 and corresponding quantitative results are compiled in Table 3. The results show that the method is fairly robust with respect to both *σ* and *p*.

Furthermore, the proposed algorithm was evaluated in the case of multiple bioluminescent sources. In the experiments, two light sources were placed in the lung with different edge-to-edge distance and 20*dB* Gaussian noise was added to the synthetic data. The reconstruction results were demonstrated in Fig. 6 and quantitative results were compiled in Table 4. In both cases, the maximum relative error between the actual and reconstructed source densities were 42%, but the reconstructed center positions were around the centers of the actual light sources.

#### 3.2. Physical phantom experiment

### 3.2.1. Experiments preparation

The physical phantom measurements were acquired based on an imaging configuration of the free-space BLT with 360° geometry projections [38, 39]. The sketch of the system structure is shown in Fig. 7. In the experiment, a cube resinous phantom with 20 mm side length was designed and manufactured. The phantom is made from polyoxymethylene, and one hole of 1.25 mm radius was drilled in the phantom to emplace the light source, as shown in Fig. 7(b). According to luminescence principle of luminescent light stick, peroxide, ester compound solutions and fluorescent dye were injected into the hole in the phantom after thorough mix, and then the red light whose central wavelength located about 660 nm was emitted due to the chemical reaction of the mixed resolutions. In bioluminescent imaging, a liquid-nitrogen-cooled back-illuminated CCD camera (Princeton Instruments VersArray 1300B) and a Micro-Nikkor 55 mm f/2.8 lens are used for data acquisition by multi-view noncontact detection mode in a dark environment. The sensitive CCD can generate 1340×1300 pixels and 16 bits dynamic range images with 20*µm*×20*µm* sized pixel even if the detectable photon flux density of bioluminescent source on the phantom surface was very weak, reaching *pico*-*Watts/mm*
^{2}. The collected bioluminescent views need to be transformed from grey-scale pixel values into corresponding numbers in physical units. Therefore, camera calibration is a prerequisite for BLT [40]. In order to do this, we used an absolutely calibrated integrating sphere of 8-inches in diameter. A filter and variable attenuator help to select a particular wavelength and to control the light level entering the sphere. For a given wavelength, gray levels are associated with varying intensity values. In our experiment, the wavelength range is about 660 nm, and the calibration formula for the CCD camera is given by the formula

where *φ* represents phantom density (*nano*-*Watts/mm*
^{2}) and pixel the pixel value. In addition, the optical parameters of the phantom were determined by a time-correlated single photon counting (TCSPC) system specifically constructed for the optical properties of the turbid medium [41]. The final measured optical parameters of the phantom at the wavelength around 660 nm were given as follows: the absorption coefficient *µ _{a}*≈0.0002

*mm*

^{-1}and the reduced scattering coefficient

*µ*′

_{s}=(1-

*g*)

*µ*≈1.0693

_{s}*mm*

^{-1}[42].

### 3.2.2. Experimental data acquisition and analysis

The physical phantom containing one source centered at (8.75,8.75,10) was mounted on a sample stage in front of the CCD camera and the experimental setup is placed in a totally dark imaging chamber to avoid external disturbance and light leaking. Under the computer control, a motorized rotation stage was used to rotate the phantom for recording the flux density with the CCD camera on the four surfaces of the phantom, as schematically shown in Fig. 7(c). During each data acquisition, one luminescent view was taken by exposing the camera for 10 seconds. The collected pixel gray levels of each luminescent view were transformed into corresponding light units according to the aforementioned calibration formula Eq. (18), and the light distribution of each view is demonstrated in Fig. 8. For BLT reconstruction, a finite-element mesh was built consisting of 1934 tetrahedron elements and 495 nodes with 225 datum nodes on the phantom surface, as shown in Fig. 9(a). Then, the proposed algorithm was applied to reconstruct the light source distribution in the phantom, the result is presented in Fig. 9(b). The reconstructed source located at (8.03,8.40,9.31), the computational time was about 37 seconds. The reconstruction results indicates that the position of the light source is accurately identified.

## 4. Discussion and conclusions

We have developed a novel BLT reconstruction algorithm to identify the bioluminescent source distribution based on Bayesian approach. The proposed algorithm was not only evaluated with a 3-D micro-CT atlas of inhomogeneous optical properties but also tested through a physical phantom experiment. All results demonstrate preferred localization and quantification both in source center and power. Our work shows the merits and potential of the Bayesian approach for optical molecular tomography. Taking into account the whole computational framework, there are several novel features which have been applied to the proposed algorithm.

First, the Bayesian approach for BLT is firstly proposed. BLT attempts to reconstruct the distribution of the bioluminescent source from boundary measurements. While Bayesian approaches are well suited to the ill-posed problem and it provides a framework to incorporate multiple types of *a priori* information to reduce the ill-posedness of BLT such as anatomical tissue information, tissue optical parameter information and permissible source information. The algorithm maximizes the log posterior probability with respect to a noise parameter and the unknown source. In each scan, the spectral projected gradient-based optimization algorithm is used. Compared with known EM algorithm for BLT, the reconstruction time reported for EM is about several hours [21], but our algorithm needs much less time than that, about tens or hundreds of seconds. Thus, the proposed algorithm is relative computationally inexpensive. Our study shows that the Bayesian approach provides an effective framework to reconstruct the distribution of bioluminescent source. Furthermore, the proposed Bayesian formulation can be extended to incorporate spectral character of bioluminescent source. Recently, the fusion of multi-modality is becoming an active research field. The Bayesian framework will present a new avenue to deal with the multi-modality information.

Second, in the Bayesian framework, a generalized adaptive Gaussian Markov random field (GAGMRF) prior model for unknown source density estimation is developed. The prior model introduce a prior information for unknown source, therefore, the ill-posedness of BLT can be further reduced. In addition, the prior model is adaptively determined on the basis of finite element analysis, which makes the proposed algorithm possible to handle a more complex geometrical model such as the real mouse phantom.

Another advantage of the proposed algorithm is its regularization property. For most algorithms existed for BLT, Tikhonov approaches are adopted and the regularization parameter has an important role in reconstruction result. Although regularization parameter can be determined with the well-known L-curve methods, the computational burden is expensive. Furthermore, the regularization parameter is not often selected accurately because of the ill-posedness of BLT. However, our algorithm is iterative and imposes regularization to the BLT problem by setting iteration numbers.

With the development of *in vivo* molecular imaging for small animals, it is necessary for reconstruction using a real animal-shape model. Furthermore, there is a need to take into consideration the heterogeneous optical properties of biomedical tissues. In our algorithm, a geometrical model is generated from a 3-D micro-CT mouse atlas, which is more accurate to evaluate the proposed algorithm than previous regular models. In addition, in order to avoid *inverse crime*, a Monte Carlo method is adopted to synthesize the measured data, which further ensures the presented algorithm evaluated properly.

In summary, we have presented a novel BLT reconstruction algorithm based on Bayesian approach and demonstrated its feasibility and potential with a micro-CT atlas. It can provide a preferred performance in view of reconstruction quality. Future work will focus on extending our algorithm to real mouse experiments to evaluate its performance as well as developing adaptive finite element meshes to enable increased resolution without overly increasing the computation burden. Relevant results will be reported later.

## Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant No. 60672050, 30672690, 30600151, 30500131, 60532050 and Funding Project for Academic Human Resources Development (PHR) under Grant No. 00627, 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, 863 Program under Grant No. 2006AA04Z216, the Joint Research Fund for Overseas Chinese Young Scholars under Grant No. 30528027, Beijing Natural Science Fund under Grant No. KZ200910005005, 4071003.

## References and links

**1. **V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listensing to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol **23**, 313–320 (2005).
[CrossRef] [PubMed]

**2. **C. H. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. **4**, 235–260 (2002).
[CrossRef] [PubMed]

**3. **B. W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt. **6**, 432–440 (2001).
[CrossRef] [PubMed]

**4. **Z. Paroo, R. A. Bollinger, D. A. Braascb, E. Ricber, and D. R. Corey, “Validating Bioluminescence Imaging as a High-Throughput, Quantitiative Modality for Assesing Tumor Burden,” Mol. Imaging **3**, 117–124 (2004).
[CrossRef] [PubMed]

**5. **G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. Meinel, “Development of the first bioluminescenct CT Scanner,” Radiology **229(P)**, 566 (2003).

**6. **E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular imaging system for small animal imaging,” Med. Phys. **30**901–911 (2003).
[CrossRef] [PubMed]

**7. **T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. **17**, 545–580 (2003).
[CrossRef] [PubMed]

**8. **M. Allard, D. Côté, L. Davidson, J. Dazai, and R. M. Henkelman, “Combined magnetic resonance and bioluminescence imaging of live mice,” J. Biomed. Opt. **12(3)**, 034018-1-11 (2007).

**9. **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**, 6765–6771 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-13-18-6756.
[CrossRef]

**10. **W. Cong and G. Wang, “Iterative method for bioluminescence tomography based on the radiative transport equation,” Proc. SPIE **6318**, 631826-1-7 (2006).

**11. **G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. **31**, 2289–2299 (2004).
[CrossRef] [PubMed]

**12. **A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. **202**, 323–345 (2005).
[CrossRef]

**13. **A. P. Gibson, J. C. Hebden, and S. R. Arridge“1,” Phys. Med. Biol. **50**R1–R43 (2005).
[CrossRef] [PubMed]

**14. **X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express **12**, 3996–4000 (2004), http://www.opticsinfobase.org/oe/abstract.cfm?URI=OPEX-12-17-3996.
[CrossRef] [PubMed]

**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. **50**, 4225–4241 (2005).
[CrossRef] [PubMed]

**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. **50**, 5421–5441 (2005).
[CrossRef] [PubMed]

**17. **N. V. Slavine, M. A. Lewis, E. Richer, and P. P. Antich, “Iterative reconstruction method for light emitting sources based on the diffusion equation,” Med. Phys. **33**, 61–68 (2006).
[CrossRef] [PubMed]

**18. **G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express **14**, 7801–7809 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-17-7801.
[CrossRef] [PubMed]

**19. **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]

**20. **G. Wang, H. Shen, W. Cong, S. Zhao, and G. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express **14**, 7852–7871 (2006), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-17-7852.
[CrossRef] [PubMed]

**21. **M. Jiang, T. Zhou, J. Cheng, W. Cong, and G. Wang, “Image reconstruction for bioluminescence tomography from partial measurement,” Opt. Express **15**, 11095–11116 (2007), http://www.opticsinfobase.org/abstract.cfm?URI=oe-15-18-11095.
[CrossRef] [PubMed]

**22. **Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finite element: methodology and simulation,” Phys. Med. Biol. **52**4497–4512 (2007).
[CrossRef] [PubMed]

**23. **Q. Zhang, L. Yin, Y. Tan, Z. Yuan, and H. Jiang, “Quantitative bioluminescence tomography guided by diffuse optical tomography,” Opt. Express **16**, 1481–1486 (2008), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-16-3-1481.
[CrossRef] [PubMed]

**24. **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]

**25. **J. Tian, J. Bai, X. Yan, S. Bao, Y. Li, W. Liang, and X. Yang, ”Multimodality molecular imaging,” IEEE Eng. Med. Biol. Mag. **27**, 48–57 (2008).
[CrossRef] [PubMed]

**26. **J. C. Ye, K. J. Webb, and C. A. Bouman, “Optical diffusion tomography by iterative coordinate descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A **16**, 2400–2412 (1999).
[CrossRef]

**27. **J. C. Ye, C. A. Bouman, K. J. Webb, and R. P. Millane, “Nonlinear multigrid algorithms for bayesian optical diffusion tomography,” IEEE Trans. Image Process. **10**, 909–922, (2001).
[CrossRef]

**28. **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 in the Living Mouse with Monte Carlo Method,” Acad. Radiol. **11**, 1029–1038 (2004).
[CrossRef] [PubMed]

**29. **S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. **20**, 299–309 (1993).
[CrossRef] [PubMed]

**30. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, 41–93 (1999).
[CrossRef]

**31. **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]

**32. **C. A. Bouman and K. Sauer, “A generalized Gaussian imaging model for edge-preserving MAP estimation,” IEEE Trans. Image Process. **2**, 296–310 (1993).
[CrossRef] [PubMed]

**33. **S. S. Saquib, K. M. Hanson, and G. S. Cunningham, “Model-based image reconstruction from time-resolved diffusion data,” Proc. SPIE **3034**, 369–380 (1997).
[CrossRef]

**34. **A. Mohammad-Djafari, “Joint estimation of parameters and hyperparameters in a Bayesian approach of solving inverse problems,” In Proceedings of IEEE Internation Conference on Image Processing **2**, 473–476 (1996).
[CrossRef]

**35. **P. E. Gill, W. Murray, and M. Wright, *Practical optimization*, (Academic Press, New York, 1981).

**36. **E. G. Birgin and J. M. Martinez, “A box-constrained optimization algorithm with negative curvature directions and spectral projected gradients,” Computing, Sup. **15**, 49–60 (2001).
[CrossRef]

**37. **E. G. Birgin and J. M. Martinez, “Large-scale Active-Set Box-Constrained Optimization Method with Spectral Projected Gradients,” Comput. Optim. Appl. **23**, 101–125, (2002).
[CrossRef]

**38. **N. Dehiolanis, T. Lasser, D. Hyde, A. Soubret, J. Ripoll, and V. Ntziachristos, “Free-space fluorescence molecular tomography utilizing 360° geometry projections,” Opt. Lett. **32**, 382–384 (2007).
[CrossRef]

**39. **H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V. Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic surface extraction,” Appl. Opt. **46**, 3617–3627 (2007).
[CrossRef] [PubMed]

**40. **T. Chen, “Digital Camera System Simulator and applications,” Ph. D. Thesis, Stanford University (2003).

**41. **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]

**42. **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/oe/abstract.cfm?URI=oe-16-25-20317.
[CrossRef] [PubMed]

**43. **S. Zhu, J. Tian, G. Yan, C. Qin, and J. Liu, “An experimental cone-beam micro-CT system for small animal imaging,” Proc. SPIE **7258**, 72582S (2009).
[CrossRef]

**44. **G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” IJ. X-Ray Sci. Technol. **16**, 225–234 (2008).