## Abstract

Multispectral bioluminescence tomography (BLT) attracts increasing more attention in the area of small animal studies because multispectral data acquisition could help in the 3D location of bioluminescent sources. Generally, BLT problem is ill-posed and *a priori* information is indispensable to reconstruction bioluminescent source uniquely and quantitatively. In this paper, we propose a spectrally solved bioluminescence tomography algorithm with an optimal permissible source region strategy. Being the most different from earlier studies, an optimal permissible source region strategy which is automatically selected without human intervention is developed to reduce the ill-posedness of BLT and therefore improves the reconstruction quality. Furthermore, both numerical stability and computational efficiency benefit from the strategy. In the numerical experiments, a heterogeneous phantom is used to evaluate the proposed algorithm and the synthetic data is produced by Monte Carlo method for avoiding the *inverse crime*. The results demonstrate the feasibility and potential of our methodology for reconstructing the distribution of bioluminescent sources.

© 2008 Optical Society of America

## 1. Introduction

The field of molecular imaging has gained considerable attention over the past few years due to its ability to reveal the molecular and cellular information *in vivo*. Hence, it becomes an increasingly important instrument for biomedical researchers to diagnose diseases, evaluate therapies, and facilitate drug development with small animals such as mouse models [1, 2, 3, 4]. As a mode of molecular imaging, optical molecular imaging, especially bioluminescence tomography (BLT), is outstanding because of its high performance, low cost and non-invasion [3, 5]. BLT is an emerging imaging technique which has been only recently developed to recover bioluminescent source distribution inside a living small animal in 3D [4, 6].

The bioluminescent signal is emitted when luciferin is combined with luciferase in the presence of oxygen and ATP, and luciferase enzymes are generally from firefly (FLuc), click beetle (CBGr68, CBRed), and Renilla reniformins (hRLuc) [7]. So signal has different emission spectra, roughly from 400*nm* to 750*nm*, which can be detected using a sensitive low light imaging systems, typically based on charge-coupled device (CCD) camera. Bioluminescent signal observed on the surface of the small animal forms the basis for tomographic reconstruction of internal bioluminescent source. However, three-dimensional (3D) bioluminescent reconstruction from boundary data is not unique and a highly ill-posed inverse problem in the general case [5].

The uniqueness research of BLT shows that *a priori* information has a quite effect on source reconstruction [5]. The commonly used *priori* information is the permissible source region strategy. So far, *a priori* permissible source region strategy has been developed for BLT reconstruction [8], which is inferred by the surface light power distribution and the heterogeneous structure of the phantom. Although promising results have been reported [4, 6, 8], it is not always reliable to infer such a permissible region especially when a single or multiple bioluminescent sources locate at half-radius or deeper depth from the small animal surface [9]. Recently, *a posteriori* permissible source region method has also been developed [9]. However, it takes several hours to select the permissible source region in the case of complex heterogenous phantom. Therefore, there is a critical need to establish a fast mathematical method for selecting the permissible source region.

Aside from the permissible source region strategy, the spectral characteristic of bioluminescent source can be considered as *a priori* information to improve BLT reconstruction [9, 10, 11, 12, 13, 14, 15], so there is an increasing interest in multispectral bioluminescence tomography. However, bioluminescent signal acquisition generally takes much longer time than that for fluorescent imaging, one exposure normally needs 5 to 10 minutes if the bioluminescent source is deep inside a small animal [16]. Although the bioluminescent signal is generally supposed to be stable in previous studies; in fact, it decays when the exposure time exceeds over one hour [16]. Traditionally, to collect the whole surface images around the imaged object, a rotating mechanism is employed to take four views (front, back, and two sides) [4, 6, 8, 16, 17]. So if we take all images in the multispectral case, one hour is not sufficient and it is possible that the observed boundary data is inaccurate. Therefore, in the multispectral case, the exposure time should be reduced with some methods. In addition, despite multispectral data acquisition could enhance the reconstruction result, dimensional disaster also arises from multispectral measurement and the large-scale data set seriously affects the reconstruction speed.

The development of a fast and robust BLT reconstruction algorithm based on practical application to biological research is necessary. In practice, in view to physical limitations as X-ray computed tomography (CT), the measurement for BLT is angle-limited and the acquired data are often incomplete [18, 19]. In this situation, partial data acquisition methods such as the single-view image technique could work. Recently, the feasibility of BLT with partial measured data has been proved theoretically and two reconstruction methods are extended to this case, but they are time-consuming and need about 5–6 hours in a single wavelength case [19]. The reconstruction time should be longer in the multispectral imaging because the computation data had increased rapidly.

Based on the consideration of above problems, we propose a reconstruction algorithm for BLT from multispectral partial boundary measurement. The most distinct difference between current and earlier tomography algorithms is the development of an optimal permissible source region strategy for reducing the ill-posedness of BLT, which is selected using iterative approach in combination of multispectral measurement without human intervention. In the algorithm, partial boundary measured data is used to reduce the exposure time which arise from multispectral measurement. Adaptive finite element algorithm is employed to reconstruct the underlying source distribution based on diffusion equation approximation. Then the linear relationship between the measured data and unknown source is established through the optimal permissible source region concept. In the optimization procedure, a spectral projected gradient-based optimization is used to reduce the reconstruction time and improve the robustness. In the numerical experiments, a heterogeneous phantom is designed to evaluate the performance of the algorithm. In the end, we discuss the relevant issues and conclude the paper.

## 2. Methods

#### 2.1. Formulation of BLT

In bioluminescence imaging, it is important to depict the propagation of the photon transport accurately in biological tissue. The bioluminescent photon is subject to both scattering and absorption. However, in the 400*nm*–800*nm*wavelength range, photon scattering predominates over photon absorption. When the bioluminescence imaging experiment is performed in a dark environment, the propagation of photons can be described by the steady-state diffusion equation and Robin boundary condition [8, 20, 21]. Taking the influence of light wavelength *λ* on tissue optical property into account, the following model is given [9, 12]:

where Ω is a bounded smooth domain in the three-dimensional Euclidean space *R*
^{3} that contains an object to be imaged; *∂*Ω is the corresponding boundary; Φ(*x*, *λ*) denotes the photon flux density [*Watts*/*mm*
^{2}]; *S*(*x*, *λ*) is the bioluminescent source density [*Watts*/*mm*
^{3}]; *µ*
_{a}(*x*, *λ*) is the absorption coefficient [*mm*
^{-1}]; *D*(*x*, *λ*)=1/(3(*µ*
_{a}(*x*, *λ*)+(1-*g*)*µ*
_{s}(*x*, *λ*))) is the optical diffusion coefficient, *µ*
_{s}(*x*, *λ*) the scattering coefficient [*mm*
^{-1}], and g the anisotropy parameter; *ν*(*x*) the unit outer normal on *∂*Ω. Given the mismatch between the refractive indices *n* for Ω and *n*
^{′} for the external medium, *A*(*x*;*n*,*n*
^{′}) can be approximately represented:

where *n*
^{′} is close to 1.0 when the mouse is in air; *R*(*x*) can be approximated by *R*(*x*) ≈ -1.4399*n*
^{-2} +0.7099*n*
^{-1}+0.6681+0.0636*n* [22]. Supposed that the measurement is conducted on some disjoint surface patches γ_{j} ⊂ *∂*Ω, *j*=1,2,…,*M*, each is smooth and connected [19]. Let γ=⋃^{M}
_{j=1} γ_{j}, the measured quantity is the outgoing flux density *Q*(*x*, *λ*) on partial boundary γ and it can be expressed:

#### 2.2. The optimal permissible source region strategy

In the practical experiment, the light which can reach the body surface of the small animals is separated into *m* bands *τ*
_{1},…,*τ*
_{m} using appropriate filters, with *τ*
_{l}=[*λ*
_{l-1}, *λ*
_{l}), *l*=1,2,…, *m*-1, *τ*
_{m}=[*λ*
_{m-1}, *λ*
_{m}]. Here *λ*
_{0} <*λ*
_{1} <…<*λ*
_{m} is a partition of the spectrum range. Based on the finite element theory [23], the weak solution of the flux density Φ(*x*, *τ*
_{l}) ∈ *H*
^{1}(Ω) on each band *τ*
_{l} is given through the Eqs. (1) and (2):

$$+{\int}_{\partial \Omega}\frac{1}{2A(\mathbf{x};n,n\prime )}\Phi (\mathbf{x},{\tau}_{l})\Psi (\mathbf{x},{\tau}_{l})dx={\int}_{\Omega}S(\mathbf{x},{\tau}_{l})\Psi (\mathbf{x},{\tau}_{l})d\mathbf{x}\phantom{\rule{.2em}{0ex}}(\forall \Psi (\mathbf{x},{\tau}_{l})\in {H}^{1}\left(\Omega \right))$$

where *H*
^{1}(Ω) is the Sobelev space and Ψ(*x*, *τ*
_{l}) is an arbitrary piece-wise test function. In the framework of adaptive finite element analysis, let {*T*
_{1},…,*T*
_{k},…} be a sequence of nested triangulation of the given domain Ω based on adaptive mesh refinement, where the sequence gradually changes from coarse to fine along with the increase in *k* [6]. The space of linear finite elements *V*
_{k} are introduced on the discretized level *T*
_{k}, satisfying *V*
_{1} ⊂ … *V*
_{k} ⊂ … … *H*
^{1}(Ω). Now, we only consider the *k*th discretized level which includes ${N}_{{T}_{k}}$ elements and ${N}_{{P}_{k}}$ vertex nodes.
${\psi}_{1}^{k}$,…, ${\psi}_{{N}_{{P}_{K}}}^{k}$ is the nodal basis of the space *V*
_{k}. Then Eq. (5) can be simplified as following matrix form on the single-band *τ*
_{l} [6, 9]:

where the components of the above matrices *K*
_{k}(τ_{l}), *C*
_{k} (*τ*
_{l}), *B*
_{k}(*τ*
_{l}), *F*
_{k}(*τ*
_{l}) are given by

Let *M*
_{k}(*τ*
_{l})=*K*
_{k}(*τ*
_{l})+*C*
_{k}(*τ*
_{l})+*B*
_{k}(*τ*
_{l}), in view to *M*
_{k}(*τ*
_{l}) is a sparse positive definite matrix, we have:

Let *A*
_{k}(*τ*
_{l})=*M*
^{-1}
_{k} (*τ*
_{l})*F*
_{k}(*τ*
_{l}), by substituting *A*
_{k}(*τ*
_{l}) into the Eq. (8), we obtain the Eq. (9):

Generally, BLT problem is an ill-posed problem. One may consider to utilize the permissible source region as *a priori* knowledge to improve the BLT reconstruction [6, 8]. In order to reduce the ill-posedness of BLT, a novel optimal permissible source region strategy is developed by using iterative approach based algorithm before BLT reconstruction. In the permissible source region, the actual source may present and in the forbidden region, the actual source is zero by definition. Commonly, the energy contribution of each spectral band *τ*
_{l} can be determined by performing a beforehand spectral analysis, that is *S*(*τ*
_{l})=*ω*(*τ*
_{l})*S*, where *ω*(*τ*
_{l}) ≥ 0 and ∑^{m}
_{l =1}
*ω*(*τ*
_{l}) ≈1, *S* denotes the total photon density. At the coarsest level (*k*=1), the whole reconstruction object Ω is assumed to be the permissible source solution *S*
_{1}, taking into account the above spectral distribution, based on Eq. (9), we have:

Where

$A=\left[\begin{array}{c}\omega \left({\tau}_{1}\right){A}_{1}\left({\tau}_{1}\right)\\ \omega \left({\tau}_{2}\right){A}_{1}\left({\tau}_{2}\right)\\ \vdots \\ \omega \left({\tau}_{m}\right){A}_{1}\left({\tau}_{m}\right)\end{array}\right],\Phi =\left[\begin{array}{c}{\Phi}_{1}\left({\tau}_{1}\right)\\ {\Phi}_{1}\left({\tau}_{2}\right)\\ \vdots \\ {\Phi}_{1}\left({\tau}_{m}\right)\end{array}\right]$

*A*
^{T} is the transpose of *A* and *A*
^{T}
*A* is an $m{N}_{{T}_{1}}\times m{N}_{{T}_{1}}$ symmetric matrix, then Eq. (10) has become a standard linear equation. Starting from an initial guess ${S}_{1}^{0}$
, the solution *S*
_{1} can be updated iteratively by:

Where *n* is iteration number and *β*
_{n} is the gradient of each iteration,

*α*
_{n} is the step size, in order to get *α*
_{n}, we construct a nonlinear least-squares optimization problem:

And

*α*
_{n} can be gotten by directly solving the Eq. (14), that is:

In the iteration process, we directly use the measured data Φ^{meas}
_{1} (*τ*
_{l}) to substitute Φ_{1}(*τ*
_{l}), because the effect of noise is low originally. With the increase of the number of iterations, the permissible source region will become smaller and smaller. When ‖β*‖*
_{n}≤ *δ* or iteration number *n*> *N*
_{max}, the iteration is terminated and the rough region of the optimal solution *S** is obtained, that is Ω*. We name it the optimal permissible source region and the sketch is shown in Fig. 1. Otherwise, *S*
_{1} is updated by the Eq. (11). Here, *N*
_{max} is maximum iteration number and the experiential typical value for δ is between 10^{-2} and 10^{-5}. Generally, for *N*
_{max} and *δ*, we choose 500 and 5×10^{-3} for practical use, respectively.

In BLT, taking into account the real physical meaning, the source density constrain can be taken as *a priori* information. Therefore, the nonnegative penalty is adopted during the iteration process.

#### 2.3. Source reconstruction

After establishing the optimal permissible source region, reformulate the Eq. (8): Φ*k*(*τ*
_{l})=*M*
^{-1}
_{ei
} (*τ*
_{l})*F*
_{k}(*τ*
_{l})*S*
_{k}(*τ*
_{l}), taking the linear relation between the unknown source value *S*
_{k} and boundary measured points
${\Phi}_{k}^{meas}$
into consideration, we have:

Where *G*
_{k}(*τ*
_{l}) can be established by deleting the rows of [*M*
_{k}(*τ*
_{l})^{-1}
*F*
_{k}(*τ*
_{l})] corresponding to unmeasured points. Incorporating the optimal permissible source region, we have:

Where

${\Phi}_{k}^{\mathrm{meas}}=\left[\begin{array}{c}{\Phi}_{k}^{\mathrm{meas}}\left({\tau}_{1}\right)\\ {\Phi}_{k}^{\mathrm{meas}}\left({\tau}_{2}\right)\\ \vdots \\ {\Phi}_{k}^{\mathrm{meas}}\left({\tau}_{m}\right)\end{array}\right],{G}_{k}=\left[\begin{array}{c}\omega \left({\tau}_{1}\right){G}_{k}\left({\tau}_{1}\right)\\ \omega \left({\tau}_{2}\right){G}_{k}\left({\tau}_{2}\right)\\ \vdots \\ \omega \left({\tau}_{m}\right){G}_{k}\left({\tau}_{m}\right)\end{array}\right]$

And *W*
_{k} is a diagonal matrix for selecting the permissible region, that is:

${W}_{k}=\mathrm{Diag}\left({w}_{k\left(11\right)},{w}_{k\left(22\right)},\dots ,{w}_{k\left(\mathrm{ii}\right)},\dots ,{w}_{k\left({N}_{{P}_{l}}{N}_{{P}_{l}}\right)}\right)$ ${w}_{k\left(\mathrm{ii}\right)}=\{\begin{array}{cc}1& \{{s}_{k\left(i\right)}\ge {\gamma}^{k}{s}_{k}^{\mathrm{max}}\}\\ 0& \{{s}_{k\left(i\right)}<{\gamma}^{k}{s}_{k}^{\mathrm{max}}\}\end{array}$

When *k*=1, *s*
_{k}(*i*) is the optimal solution *S*
^{*} and *s*
^{max}
_{k} is its maximum. *k* ≥ 2, *s*
_{k(i)} and *s*
^{max}
_{k} are the reconstructed results prolonged from *k*-1th level and the corresponding maximum, that is

where I^{k}
_{k}-1 is the ratio operator from *k*-1 to *k*, *γ*
^{k} is the ratio factor, *S*
_{k-1} is the reconstruction result on (*k*-1)th level. By retaining the columns of *G*
_{k}
*W*
_{k} corresponding to the permissible region *S*
^{p}
_{k}, the final form of the linear system between the measurable boundary flux Φ_{meas}^{k} and *S ^{p}_{k}* is obtained:

Since the matrix *A*
_{k} is a severely ill-conditioned matrix because of the ill-posedness of BLT, the surface measured data is corrupted by noise, it is not practical to directly solve for *S*
^{p}
_{k} from linear system (19). Then the following *k*th objective function is established:

where *s*
^{k}
_{inf} and *s*
^{k}
_{suf} are the kth level lower and upper bounds of the source density; Λ is the weight matrix, *V*Λ=*V*
^{T}Λ*V*; *ρ*
_{k} the regularization parameter. *S*
^{init}
_{k} is initial value at the *k*th level. For BLT, Θ*k*(*S*
^{p}
_{k}) is a large-scale optimization problem with box-bound constrains, there-fore, a spectral projected gradient-based large-scale optimization algorithm is modified to solve the least square problem [9, 24, 25, 26]. In the algorithm, as for *posteriori* error estimation and local mesh refinement, similar disposal methods are adopted according to [6, 9]. The current gradient norm $\parallel {g}_{\underset{.}{\Theta}k}^{i}\parallel $ and the iteration number *N*
^{i}
_{k} on each level are selected as switch condition indexes for triggering local mesh refinement. The discrepancy between the measured and computational boundary nodal flux data and the number *k* of mesh refinement are used as stop criterion. Finally, the flow chart of the algorithm is shown in Fig. 2.

## 3. Numerical simulation

As far as phantom is concerned, it plays an important role in evaluating reconstruction algorithm and various phantoms have been designed for bioluminescence tomography. In previous studies, the assumption of a homogeneous optical background had demonstrated the disadvantage and limitation for BLT reconstruction [13]. It’s not proper and accurate to evaluate the reconstruction algorithm.

Therefore, a heterogeneous phantom of 30*mm* height and 10*mm* radius was designed to evaluate the proposed algorithm [6]. It was made up of four ellipsoids and one cylinder to represent muscle (M), lungs (L), heart (H), bone (B) and liver (Li), as shown in Fig. 3(a). Based on the emission spectral distribution, the spectrum [500*nm*–750*nm*] can be divided into five discrete bins with steps of 50*nm*. In the experiments, the optical properties of each component are assumed as *priori* information and optical property parameters of each bin are compiled in Table 1 [9].

Because BLT is an inverse source problem, it is necessary to acquire the surface measured data objectively and reliably. In view to *inverse crime*, the synthetic data was produced by molecular optical simulation environment (MOSE) which was developed using Monte Carlo method [27]. In the single source simulation, the solid spherical source with 1*mm* radius and source intensity of 0.238*nano*-*Watts*/*mm*
^{3} was centered at (-3,5,15) inside the right lung as shown in Fig. 3(a). In MOSE, the heterogeneous phantomwas discritized by triangular elements with an average element diameter of about 1*mm*, shown in Fig. 3(b). In the simulation, the bioluminescent source of each band was sampled by 1.0×10 ^{6} photons and the default initial guess *S*
^{0}
_{1} for searching the optimal permissible source region was set 0.01. When *k*=1, the ratio factor *γ*
^{k}=0; *k* > 1, *γ*
^{k} was initially set 10^{-1} for selecting permissible source region and was changed by multiplying a factor of 10.0 with the increase of level *k*. In the optimization process, *s*
^{k}
_{inf}, *s*
^{k}
_{suf} and *S*
^{init}
_{k} are equal to 0.0, 10000 and 1.0×10^{-5}, respectively. The stopping threshold ε_{Φ} and the maximum number of mesh refinement *K*
_{max} were set to 1.0×10^{-8} and 3, respectively. Noted that in the reconstruction procedure, the initial discretization of phantom was coarse volumetric mesh shown in Fig. 3(c). In terms of the present condition, a single-view image technique is feasible for acquiring partial measured data. So in our simulations, the data only from the front view shown in Fig. 3(d) is used in case of partial measurement.

#### 3.1. Experiment Results

In the partial measurement experiment, the default initial guess *S*
^{0}
_{1} was used. As for switch condition indexes, the terminated gradient norm *ε*
_{g} and the maximum iteration number *N*
^{max}
_{k} were set to 1.0×10^{-7} and 5000, respectively. *ρ*
_{k} was 1.16×10^{-9} at each level. First, using partial measured data in [500*nm*,550*nm*], the optimal permissible source region was established and the reconstruction was performed, corresponding result was shown in Fig. 4(a). Despite the use of optimal permissible source region, the monochromatic measurement-based tomographic reconstruction result was still undesirable no matter how we changed the initial guess *S*
^{0}
_{1}, *δ*, *ρ*
_{k} or maximum iteration number *N*
_{max}. However, using multispectral partial measured data, BLT reconstruction showed preferable result with the proposed algorithm after one level refinement and the result was shown in Fig. 4(c). By contrast, when complete measured data from four views was used for selecting the optimal permissible source region, the corresponding result was demonstrated in Fig. 4(e).

In order to qualify the results, both the absolute source position error (ε) and the relative source density error (*ξ*) are defined: $\epsilon =\sqrt{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}+{\left(z-{z}_{0}\right)}^{2}}$, *ξ*=|*S*
_{recon}-*S*
_{real}|/*S*
_{real}, where (*x*, *y*, *z*) is the reconstructed center of source and (*x*
_{0}, *y*
_{0}, *z*
_{0}) is the actual center of source. *S*
_{recon} and *S*
_{real} are reconstructed source density and actual source density, respectively, the unit is *nano*-*Watts*/*mm*
^{3}. Quantitative results about both the location and density of source were shown in Table 2.

In order to compare with other algorithm, BLT reconstruction was also conducted with the

method proposed in literature [9], the result was demonstrated in Fig. 4(g) and quantitative result in Table 2. Quantitative comparison shows that the proposed algorithm has a better performance. Otherwise, the reconstruction speed is also a standard for evaluating algorithm. In order to get the result in Fig. 4(g), the total time cost used the algorithm mentioned in literature [9] was about 11 hours, but our method took much less time, about 3 hours. This was because the system matrix at first level was reduced from 3350×6878 to 1690×5397 with the employment of partial measurement and the optimal permissible source region. These work was all implemented on the same desktop computer with Pentium 4 3.4GHz and 1GB RAM.

From the above comparisons, it is obvious that the proposed algorithm can recover the bioluminescent source distribution relative accurately. This means that the measured data can be acquired only on part body of a small animal, and then the time for acquiring data can be dropped. Meanwhile, the system matrix dimension can also be reduced. Furthermore, our algorithm needs relative shorter reconstruction time and it is more appropriate for practical application. Nevertheless, spectrally resolved BLT is indispensable and provides more sufficient *a priori* information.

#### 3.1.2. BLT reconstruction with noisy data

In order to evaluate the stability and robust of the proposed algorithm, noise effect was considered. In the experiment, white Gaussian noise with different level was added to the synthetic data and corresponding reconstructions were carried out. The corresponding results were shown in Fig. 5 and Table 3. The results reveal the stability of proposed algorithm under noise. When the noise level is lower (5% and 10%), the method has a better performance not only in reconstruction position but also source density. However, the noise level is higher (20% and 50%), although the reconstructed source density is perturbed, the position is not affected. After we have carried out intensive numerical simulation, we find two reasons to explain this. One reason is that the noise effect is low originally for selecting the optimal permissible source region. The other is that source reconstruction benefits from the adaptive from-coarse-to-fine mesh sequence.

#### 3.1.3. Consideration of different initial guesses

It is well known that a big disadvantage of most iterative approach based algorithm is sensitive to the initial guess (*S*
^{0}
_{1}). It is regarded as a criterion to test the stability of the proposed algorithm. Due to the limited space, representative results using different initial guesses of iterative approach are given. When the initial guess was 0.01, the result had been reported, shown in Fig. 4(c). When the initial guesses were set 0 and 0.1, the reconstructed results all indicated the algorithm was stable and robust to initial guess. The results were shown in Fig. 6 and Table 4. When the initial guess was equal to zero, the maximum reconstruction density was 0.142*nano*-*Watts*/*mm*
^{3}. Although the density had a great difference with real density, the recovered center position of bioluminescence source was satisfactory, that was (-1.54,5.35,14.57).

#### 3.1.4. Optical property errors consideration

Currently, with further investigation of BLT, multimodality imaging fusion has been become a hot spot [4, 9, 13, 28]. Besides those, Diffusion Optical Tomography (DOT) is also appropriate for BLT which can be used to acquire optical properties of mouse tissues. In previous algorithms, optical property parameters of biomedical tissues were generally taken as *a priori* information to reduce the ill-posedness of BLT and deal with the nonuniqueness of BLT [5, 9, 12]. With DOT in combination with BLT, the problem can be solved. In view to the error range of DOT reconstruction quantification accuracy [9, 29], ±50% errors for all tissues were used to evaluate the proposed algorithm and reconstruction results were demonstrated in Fig. 7. The reconstructed position of source for -50% and +50% were (-3.34,3.87,14.08) and (-1.64,5.55,15.18), respectively. The maximum absolute source position error is 1.50*mm*, which accounts for the capability for tolerating optical property errors.

Furthermore, the possible effect of two sources closer togetherwas studied. The two spherical sources of 1.0*mm* radius and 0.238*nano*-*Watts*/*mm*
^{3} power density were placed in the right lung with edge-to-edge distance of 2*mm*. Reconstruction results were shown in Talbe 5. In the case of +50% error for all tissues, the center position of two sources both had an offset of about 2*mm*, but the sources could be distinguished. However, when the error for all tissues is -50%, the reconstructed center position of a source had a large offset and the two sources could not be resolved. The effect of two sources with different edge-to-edge distances is necessary to further explore.

## 4. Discussions and conclusion

Given the difficulty that there is no unique solution to the BLT problem in general case and the need for fast BLT reconstruction algorithm, we have developed the novel optimal permissible source region strategy and present a BLT reconstruction algorithm to identify a 3D bioluminescent source distribution, finally evaluate its performance in numerical simulation with quantitative results. Simulation studies have indicated that the algorithm can exhibit very favorable performance. In the reconstruction process, the optimal permissible source region has played an important role to reduce the ill-posedness of BLT as well as enhance the numerical stability and computational efficiency. In addition, in the case of hyper- and multi-spectral measurement, the exposure time increases and dimension disaster arises. With the proposed algorithm, the problems can be dealt successfully.

A challenge problem for BLT is the reconstruction speed. In the algorithm, an optimal permissible source region is firstly established before reconstruction, so reconstruction time and memory cost are reduced. Secondly, adaptive finite element algorithm improves the reconstruction speed. Furthermore, the finite-element based reconstruction method makes it possible to handle a more complex geometrical model such as the real mouse phantom. The heterogeneous phantom experiments have illustrated that the method is fairly robust with respect to both strong Gaussian noise and initial guess. More importantly, the data in the algorithm were produced using Monte Carlo-based method, and hence they are totally free of the well-known *inverse crime* [8, 30]. Therefore, the presented algorithm can be properly evaluated. In the BLT prototype, multimodality fusion such as DOT, Micro-CT and BLT has become a research hot spot. Taking into account the high tolerance of optical property errors, the presented algorithm is suitable for future multimodality BLT system. However, because the algorithm depends on multiple factors, it is possible to get favorite results with try and error method and the optimal setting is still under investigation. As measure of performance, we considered the difference of locations of the center as well as the difference in intensity between the reconstructed and real source. Our algorithm demonstrates a satisfactory performance in this case, however, considering the difference in volume as well would be as important. In other words, there are remaining issues needed to be considered.

In summery, we have presented a novel strategy for selecting the permissible source region, and proposed a BLT reconstruction algorithm for practical application. It can provide a better performance in view of reconstruction quality, robustness and speed. Now, we are making an effort to construct a multimodality BLT system based on the proposed reconstruction algorithm, the new system will provide the potential to enhance our understanding of disease and drug activity during preclinical and clinical drug development, which could aid decisions to select candidates that seem most likely to be successful or to halt the development of drugs that seem likely to ultimately fail [31]. Corresponding 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, 60621001 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 (YZ0642, 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. 4071003.

## 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. **C. H. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng.**4**, 235–260 (2002).

**3. **G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. F. Meinel, “Development of the first bioluminescence ct scanner,” Radiology **229(P)**, 566 (2003).

**4. **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). [CrossRef] [PubMed]

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

**6. **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). [CrossRef] [PubMed]

**7. **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-1–9 (2005). [CrossRef]

**8. **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). [CrossRef] [PubMed]

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

**10. **W. Han, W. Cong, and G. Wang, “Mathematical Study and Numirical Simulation of Multispectral Bioluminescence Tomography,” Int. J. Biomed. Imaging 2006, 54390 (2006).

**11. **O. Coquoz, T. L. Troy, D. Jekic-McMullen, and B. W. Rice, “Determination of depth of in vivo bioluminescent signals using spectral imaging techniques,” Proc. SPIE **4967**, 37–45 (2003). [CrossRef]

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

**13. **G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. **50**, 4225–4241 (2005). [CrossRef] [PubMed]

**14. **A. X. Cong and G. Wang, “Multispectral Bioluminescence Tomography: Methodology and Simulation,” Int. J. Biomed. Imaging **2006**, Article ID 57614, 7 pages, 2006. doi:10.1155/IJBI/2006/57614. [CrossRef]

**15. **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). [CrossRef] [PubMed]

**16. **G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of Multiview and multispectral Data,” Int. J. Biomed. Imaging, **2006**. Article ID 58601, 8 pages, 2006. doi:10.1155/IJBI/2006/58601. [CrossRef]

**17. **X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express **12**, 3996–4000 (2004). [CrossRef] [PubMed]

**18. **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).

**19. **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). [CrossRef] [PubMed]

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

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

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

**23. **S. S. Rao, *The Finite Element Method in Enginering*, (Butterworth-Heinemann, Boston, 1999).

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

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

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

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

**28. **Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental Study on Bioluminescence Tomography with Multimodality Fusion,” Int. J. Biomed. Imaging **2007**, 86741 (2007). [CrossRef]

**29. **V. Ntziachristos, A. H. Hielscher, A. G. Yodh, and B. Chance, “Diffuse Optical Tomography of Highly Heterogeneous Media,” IEEE Transactions on Medical Imaging **20**, 470–478 (2001). [CrossRef] [PubMed]

**30. **S. Holder, *Electrical Impedance Tomography* (Institute of Physics Publishing, Bristol and Philadelphia, 2005).

**31. **J. K. Willmann, N. V. Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery **7**, 591–607 (2008). [CrossRef]