## Abstract

Bioluminescence tomography (BLT) is used to localize and quantify bioluminescent sources in a small living animal. By advancing bioluminescent imaging to a tomographic framework, it helps to diagnose diseases, monitor therapies and facilitate drug development. In this paper, we establish a direct linear relationship between measured surface photon density and an unknown bioluminescence source distribution by using a finite-element method based on the diffusion approximation to the photon propagation in biological tissue. We develop a novel reconstruction algorithm to recover the source distribution. This algorithm incorporates *a priori* knowledge to define the permissible source region in order to enhance numerical stability and efficiency. Simulations with a numerical mouse chest phantom demonstrate the feasibility of the proposed BLT algorithm and reveal its performance in terms of source location, density, and robustness against noise. Lastly, BLT experiments are performed to identify the location and power of two light sources in a physical mouse chest phantom.

©2005 Optical Society of America

## 1. Introduction

Small animal imaging has become an important tool for biomedical research at the anatomical, functional, cellular and molecular levels. To study a small animal using molecular imaging techniques, the animal organ/tissue is typically transfected with a reporter gene in a viral promoter. This mechanism has been used, for instance, in the detection of cancer/metastases [1, 2]. In comparison with conventional imaging techniques like X-ray computed tomography (CT), magnetic resonance imaging (MRI) and positron emission tomography (PET), bioluminescent imaging has the capability to reveal molecular/cellular activities directly, and it is also much more sensitive in imaging gene expression [3].

Bioluminescent imaging employs luciferase enzymes for the real-time in vivo detection of tagged cells in living animals. After luciferin is applied to an animal by injection, those cells in the organism that express the luciferase transgene emit photons of light and the broad emission spectra of some luciferases contain significant red components [2]. The resultant light intensity is directly proportional to the number of luciferase molecules and the concentration of the luciferin. Photon propagation in the biological tissue is subject to both scattering and absorption. Scattering is due to changes in the refractive index at the cell membranes and internal organelles [4]. In our experimental spectral range (~650nm), the absorption varies greatly with the type of cell or tissue and is largely governed by the amount of hemoglobin present [4, 5]. In small animal studies, a significant number of bioluminescent photons can escape the attenuating environment, and they can be detected using a highly sensitive charge-coupled device (CCD) camera. Because the biological tissue does not emit photons and no external light source is required for excitation, the background noise in bioluminescent imaging is very low.

Photon propagation in biological tissue is governed by the radiative transfer equation (RTE) [6]. However, the RTE is computationally expensive in practical bioluminescent imaging. Given the dominance of scattering over absorption in this context, diffusion theory provides a quite accurate description of the imaging model, and given appropriate boundary conditions, it can be applied to find good estimates of the surface flux density measures [7, 8, 9]. Recently, the diffusion approximation with a modified diffusion coefficient has been verified through steady-state measurement in highly absorbing and scattering media [10]. The modeling of photon transport in tissue has demonstrated that a bioluminescent cell count as low as a few hundred can be detected at subcutaneous tissue sites and that approximately 10^{6} cells are required to generate signals that are detectable through a tissue of 2cm in thickness[4].

Our bioluminescence tomography (BLT) technology utilizes the imaging modality fusion approach. In combination with bioluminescent data acquisition, the anatomical structures of a small animal and the associated optical properties are also obtained using a CT/MRI scanner or some alternative method. The resultant anatomy of the small animal is then segmented into its major components, such as heart, lungs, liver, stomach, bones, etc. Practical techniques have been developed to determine the optical parameters by using diffuse reflectance measurements from the biological tissues [11, 12]. Diffuse optical tomography (DOT) can also be used to reconstruct the spatially variable optical parameters with *a priori* information [13]. In this feasibility study, published optical parameters (absorption, reduced scattering) for the major anatomical components were used to build a geometrical model of the mouse. Therefore, the BLT problem is here formulated as an inverse source problem based on the diffusion equation. Recently, Wang et al. described the BLT principles and reported the uniqueness of solution dependent on *a priori* knowledge [14, 15]. Partial numerical studies were also reported in the recent literature [16, 17]. In the present work, we develop a novel reconstruction algorithm to identify the bioluminescent source distribution from the measured external photon density. This algorithm incorporates *a priori* knowledge to define the permissible source region to enhance numerical stability and efficiency. In the second section, the diffusion equation is discretized by finite element analysis to yield the corresponding matrix equation. In the third section, the inversion algorithm for BLT is presented in detail. In the fourth section, numerical and experimental results are described to demonstrate the feasibility of our reconstruction method. In the last section, relevant issues are discussed and conclusions drawn.

## 2. Diffusion equation and its finite-element discretization

#### 2.1 Diffusion approximation

In bioluminescence imaging, biological entities (e.g., tumor cells, genes) are tagged with luciferase enzymes. When the luciferase is combined with the substrate luciferin, oxygen and ATP, a biochemical reaction occurs that transforms part of the chemical energy into bioluminescent photons with a wavelength of about 600nm [2]. Photon scattering predominates over absorption in the biological tissue. The photons’ propagation can be described by the following steady-state diffusion equation [7, 8, 9]:

where Ω is the region of interest; Φ(**x**) represents the photon density [Watts/mm^{2}] at location x; *S*(**x**) the bioluminescent source density [Watts/mm^{3}]; *μ*_{a}
(**x**) the absorption coefficient [mm^{-1}]; *μ*_{s}
(**x**) the scattering coefficient [mm^{-1}]; and g the anisotropy parameter. Ideally, the optical parameters *μ*_{a}
(**x**), *μ*_{s}
(**x**) and g are obtained from optical tomography/sensing techniques [5, 11, 12, 13].

#### 2.2 Boundary conditions and measurement

Assuming that the bioluminescent imaging experiment is performed in a totally dark environment, no external photon travels into Ω through its boundary ∂Ω. Taking into account the mismatch between the refractive indices *n* for Ω and *n*′ for the surrounding medium, the boundary condition for the diffusion equation (1) can be expressed as [ 18, 19]:

where **v** is the unit outer normal on ∂Ω. In the experiment, the medium surrounding Ω is air, for which *n*′ is approximately 1. Therefore,

where *R* can be approximated with *R* ≈ -1.4399*n*
^{-2} + 0.7099*n*
^{-1} + 0.6681 + 0.0636*n* [19]. The measured quantity is the outgoing photon density on ∂Ω. [18]:

#### 2.3 Finite-element discretization

The governing equations (1)–(2) for Φ(**x**) can be equivalently represented as the following weak form [9] that Φ(**x**) satisfies

$$\phantom{\rule{2em}{0ex}}+\frac{{\int}_{\partial \Omega}\Phi \left(\mathbf{x}\right)\Psi \left(\mathbf{x}\right)}{(2A\left(\mathbf{x};n,n\prime \right))d\mathbf{x}}={\int}_{\Omega}S\left(\mathbf{x}\right)\Psi \left(\mathbf{x}\right)d\mathbf{x}$$

for an arbitrary test function Ψ(**x**). Following the standard finite-element method [20], Ω is discretized with *T* vertex nodes (*N*
_{1},*N*
_{2},…,*N*_{T}
) and *N*
_{e} elements, denoted as Ω^{(l)}(*l* = 1,2,Ω…,*N*
_{e}), such that Ω = ∪^{Ne}
_{l=1}Ω^{(l)}. Then Φ(**x**) can be approximated with a piecewise polynomial function [9, 19–21],

where *ϕ*_{k}
is the nodal value of Φ(**x**) on the *k*-th node *N*_{k}
; and *φ*_{k}
(**x**) the nodal basis function with support over the elements Ω^{(l)} , which have the node *N*_{k}
as a common vertex, that is supp(*φ*_{k}
(**x**)) = ∪_{Nk
∊Ω(l)}Ω^{(l)} . Similarly, the source function *S*(**x**) is approximated as [20],

where *S*_{k}
are the values of *S*(**x**), and *γ*_{k}
(**x**) are the interpolation basis functions (*k* = 1,2,…,*N*_{s}
); *N*_{s}
is the number of interpolation basis functions for the source. By inserting Eqs. (5)–(6) into Eq. (4) and using the nodal basis functions *φ*_{k}
(**x**) (*k* = 1,2,…,*T*) as the test functions, we obtain the following matrix equation

where the components of the system matrices are given by

According to the convergence theorem for the finite element method [9, 21], the solution of Eq. (7) converges to the exact continuous solution of the diffusion equation (1) when the maximum size of the elements approaches zero. In the BLT reconstruction, mesh sizes and shape functions should be chosen to balance the accuracy and the computational efficiency.

## 3. BLT reconstruction method

To describe the BLT reconstruction method, we rewrite Eq. (7) as

where Φ^{m} represents the measurable nodal photon density on the boundary ∂Ω, and Φ^{*}the photon density on internal nodes. In our experiments, Φ^{m} is computed from the surface flux image captured with a CCD camera. The source vector *S* is divided into two parts: *S*^{p}
in the permissible source region Ω_{s} and *S*
^{*} in the forbidden region. Clearly, *S*
^{*} is zero by definition. Note that the permissible source region Ω_{s} must be selected to be sufficiently large to contain the actual source region where the reporter genes may present. The sub-metrics, M_{11}, M_{12} and M_{22} of [M] and F_{11}, F_{12}, F_{21} and F_{22} of [F], are arranged in consistence with *S*^{p}
,*S*
^{*}, Φ^{m} and Φ^{*}. Thus, Eq. (9) can be reduced to

Since the matrix [M] is symmetric and positive definite, Φ^{m}can be uniquely determined from

assuming that the source *S*^{p}
is known. In a BLT experiment, the output photon flux *Q*(**x**) from the specimen is captured with a CCD camera. By Eq. (3), photon density on the surface of the specimen is obtained from *Q*(**x**), and is discretized as a vector Φ^{meas}, in consistency with vector Φ^{m} in Eq. (11). The purpose of the BLT is to identify the source *S*^{p}
from measured Φ^{meas} . Generally, since the measured data in bioluminescent imaging are corrupted by noise, it is not practical to directly solve for *S*^{p}
from linear system (11) with Φ^{m} replaced by Φ^{meas}. A maximum-likelihood approach is often employed to obtain the optimal solution to Eq. (11). The data Φ^{meas} can be expressed as Φ^{meas} = Φ^{m} + **N** , where Φ^{m} is the prediction of the emission fluence on the small animal body surface given by Eq. (11), and **N** is the noise distribution reflecting all the randomness and errors associated with the imaging process. The inherent data noise **N** is typically characterized by the Poisson distribution. In practice, it is common to model the noise as a Gaussian distribution since the Gaussian distribution is often a very good approximation of the Poisson distribution when the photon rate is sufficiently high. Hence, the Gaussian noise model can be used with both the mean and the variance equal to the actual measurement. The density function for a specific measure is given by [22]:

where ${\mathrm{\Phi}}_{k}^{\mathit{\text{meas}}}$ is the measured photon density at the k-th detector position(*k* = 1,2,…,*M*). We assume that the noise is independent for different detectors so that the covariance matrix C of data Φ^{meas} is diagonal and expressed by *C*_{kk}
= ${\mathrm{\Phi}}_{k}^{\mathit{\text{meas}}}$(*k* = 1,2,…,*M*). To simplify the notation, we define the diagonal matrix **W** as **W** = C^{-1}/2 . The data likelihood is then formulated as a function of multi-variables:

where *S*^{p}
= {*s*
_{1},*s*
_{2},…,*s _{Ns}*}

^{T}is the light source distribution to be reconstructed. This Bayesian framework should also be subject to the previously mentioned constraints, such as that the source distribution must stay in the permissible region Ω

_{s}, which can be determined in reference to the bioluminescent signals and

*a priori*knowledge available from a specific biomedical application. Then, we have the following objective function for the maximum likelihood of the data:

The BLT reconstruction is reduced to minimize the quadratic multivariate function Θ(*S*^{P}
),

In our numerical computation, a modified Newton method and an active set strategy for bound constrained optimization [23] were adopted to solve the minimization problem (15).

## 4. Preliminary results

#### 4.1 Numerical simulation

A cylindrical heterogeneous phantom of 20mm diameter and 12mm height was designed for numerical simulation; it contained four kinds of materials to represent muscle (M), lungs (L), heart (H), and bone (B), respectively. Cross-sections of the phantom are shown in the following Figs.(1~5). Optical parameters documented in the literature [5] were assigned to each of the four components; see the Table 1. The phantom was discretized into 11340 wedge elements and 6576 nodes. Simulated measurement data of photon density from 1024 datum nodes, which are distributed along the entire side surface of the phantom, are collected for the source reconstruction.

### 4.1.1 Light source reconstruction

We designed a number of source configurations to reveal the BLT image resolution. The first model contains only a single source made of six adjacent wedge elements with a total volume of 1.7mm^{3} and a power density of 200pico-Watts/mm^{3}, as shown in Fig. 1. In the second model, the two sources are in the left and right lungs, respectively; each source is made of six adjacent wedge elements with a total volume of 1.7mm^{3} and a density of 200pico-Watts/mm^{3}, as presented in Fig. 2. The third model supports two light sources with a 1.8mm separation embedded in the left lung region of the phantom, as shown in Fig. 3. In the fourth model, two sources are in the left lung and one source is in the right lung, see Fig. 4. Each light source in Figs. 3 and 4 occupies about 1.7mm^{3} and has a photon density of 200pico-Watts/mm^{3}. The photon density data on the side surface of the phantom was generated from a finite-element forward model, with 10% Gaussian noise added to simulate measured data *Q*(**x**) . The reconstruction algorithm was implemented in MATLAB with the computationally intensive parts coded in C. Then, a multivariate optimization procedure subject to non-negativity was conducted to minimize the difference between the predicted and measured photon densities on the phantom surface. The initial density distribution was set to 7.5e-6 pico-Watts/mm^{3}. Figs. 1–4 display the true and reconstructed source distributions for the four source models, respectively. The reconstruction results indicate that the locations of the light sources are accurately identified but that the recovered average source densities are subject to relative errors of about 5%. This finite-element base reconstruction method is computationally tractable. Using a finite-element mesh with a few thousand nodes and a permissible source region with a few hundred elements, the BLT reconstruction program coded in MATLAB/C took about 5 minutes on our desktop computer (AMD Athlon MP 2800+ AT/AT COMPATIBLE 2G RAM).

### 4.1.2 Spatial resolution study

To demonstrate the tomographic capability of our reconstruction method, we embedded four light sources in the left lung region of the phantom with these sources separated by 1mm between the first and the second, 2mm between the second and the third, and 3mm between the third and the fourth sources, as shown in Fig. 5. Each source had a flux density of 200pico-Watts/mm^{3}. To mimic real measurement data in the bioluminescence experiment, the output photon density data were corrupted on the side surface of the phantom at a 10% Gaussian noise level to synthesize noisy photon density data *Q*(**x**) . Then, the proposed reconstruction algorithm was applied to reconstruct the light source distribution from the photon density on the side surface of the phantom. Figure 5 shows the true and reconstructed light source distributions, respectively. The reconstruction results indicate that the positions of the light sources are accurately identified, with a maximum relative error 12.5% in terms of source strength.

### 4.1.3 Permissible region study

To evaluate the effects of the permissible source region and the data noise level on the BLT reconstruction quality, two bioluminescent sources were embedded into the left lung region of the phantom. One source had a flux density of 100pico-Watts/mm^{3} and was uniformly distributed in the 6507-th volume element, whose center located at (-6.59, 1.76, 6.0). The other source had a flux density of 200pico-Watts /mm3 and was uniformly distributed in the 6510-th volume element, whose center located at (-6.60, -1.75, 6.0). To simulate the real measurement data, 5%, 10% and 15% Gaussian noises were added to the datasets on the side surface covering phantom.

Because BLT is a highly under-determined problem, it is very effective to regularize the solution by restricting the source distribution within a permissible source region. To compare BLT outcomes with different permissible source regions, three permissible source regions were utilized, which are expressed as

$${\Omega}_{s}^{2}=\left\{\left(x,y,z\right)|x<0,5.6<z<8.0,\left(x,y,z\right)\phantom{\rule{.2em}{0ex}}\u220a\phantom{\rule{.2em}{0ex}}L\right\}$$

and ${\mathrm{\Omega}}_{s}^{3}$ ={(*x*,*y*,*z*) |*x* < 0, 4.8 <*z*<8.0,(*x*,*y*,*z*)∊L}. There are 154, 308 and 462 volume elements in these three permissible regions, respectively. The bioluminescent source distribution was reconstructed in these three permissible regions ${\mathrm{\Omega}}_{s}^{1}$, ${\mathrm{\Omega}}_{s}^{2}$ and ${\mathrm{\Omega}}_{s}^{3}$, for different noise levels, respectively. The BLT results associated with the different permissible region sizes and data noise levels are summarized in Table 2.

Table 2 shows the variations in the reconstructed source’s power errors. The smaller permissible region ${\mathrm{\Omega}}_{s}^{1}$ possesses a maximum source power error of 13.0%, while the larger one ${\mathrm{\Omega}}_{s}^{3}$ possesses a maximum source power error of 33%. In the permissible source regions ${\mathrm{\Omega}}_{s}^{1}$ and ${\mathrm{\Omega}}_{s}^{2}$ , the reconstructed source location can be exactly determined, as shown in Fig.6–9. For the largest permissible source region ${\mathrm{\Omega}}_{s}^{3}$ , there exist appreciable discrepancies in source distribution and shape: for instance, in the high noise situation (i.e., 15% noise), in Fig. 9(c), the reconstructed sources are scattered in a small region positioned closely around the actual source position. The small region is approximately bounded by a sphere with a radius of 1.5mm, which clearly suggests the actual location of the source. These results indicate that the smaller the permissible source region, the more accurately we can recover the source position and strength. Based on our numerical results, it is clear that a promising BLT strategy is to use a multi-scale reconstruction procedure. That is, subsequent BLT reconstructions should be carried out within progressively reduced permissible source regions.

#### 4.2 Heterogeneous phantom experiment

### 4.2.1 CCD camera calibration

In bioluminescent imaging, a CCD camera is used for data acquisition on the heterogeneous phantom surface. The collected bioluminescent views need to be transformed from grey-scale pixel values into corresponding numbers in physical units. Hence, camera calibration is a prerequisite for BLT [24]. To do this, we used an absolutely calibrated integrating sphere of 8-inches in diameter, which contains a night vision monitor resolving 10e-7 F-L (~5fempto-Watts/mm^{2}) or the equivalent (LR-8-LC, 8” low level output sphere system, SphereOptics, Contoocook, New Hampshire). The sphere is illuminated with a tungsten lamp. A filter and variable attenuator help to select a particular wavelength with FWHM 20nm and to control the light level entering the sphere. For a selected wavelength, gray levels are associated with varying intensity values. For the wavelength range of interest, 600–650nm, a calibration formula for the CCD camera is given by *φ* = *pix*×0.377 pico-Watts/mm^{2} , where *φ* represents photon density and *pix* the pixel value.

### 4.2.2 Mouse chest phantom

A cylindrical heterogeneous mouse chest phantom of 30mm height and 30mm diameter was designed and fabricated. It consisted of four different materials high-density polyethylene (8624K16), nylon 6/6 (8538K23), delrin (8579K21) and polypropylene (8658K11) (McMaster-Carr supply company, Chicago, IL, US) to represent muscle (M), lungs (L), heart (H) and bone (B), respectively, as shown in Fig. 10. A luminescent light stick (Glowproducts, Canada) was selected as the testing source. The stick consists of a glass vial containing one chemical solution and a larger plastic vial containing another solution with the former being embedded in the latter. By bending the plastic vial, the glass vial can be broken to mix the two solutions after which luminescent light is emitted due to the reaction of the solutions. The particular dye in the chemical solution is red light, and it can last for approximately 4 hours at an emission wavelength range between ~650nm and ~700nm, which is close to that of the red spectral region of the luciferase. Two small holes of diameter 0.6mm and height 3mm were drilled in the phantom with their centers at (-9.0, 1.5, 15.0) and (-9.0, -1.5, 15.0) in the left lung region of the phantom, respectively. Two red luminescent liquid filled catheter tubes of 1.9mm height and 0.56mm diameter were placed inside the two small holes, respectively. We measured the total power of the red luminescent liquid filled polythene tubes with the CCD camera. They were 105.1 nano-Watts and 97.4 nano-Watts, respectively.

### 4.2.3 Optical parameters

Since optical parameters are needed for BLT, we have to determine them for the four components (M, H, L and B) of the physical phantom. Cylinders of the above-mentioned materials were made with diameter 20mm and height 20mm for determination of the optical parameters. The side surface of the cylindrical homogeneous specimen was blackened. The two opposite bottom surfaces of the specimen were left uncovered. The light output from the exit port of the integrating sphere was guided through the optic fiber and used for illumination. The other end of the optic fiber was inserted into a small hole of 10mm depth at the center of one specimen bottom surface. Then, the specimen was placed on a sample holder in front of the nitrogen-cooled CCD camera (Princeton Instruments VersArray: 1300B, Roper Scientific Inc, Trenton, NJ) in a dark environment. It captured the output photon density on the other bottom surface of the specimen with an exposure time of 30 seconds. After the data acquisition, the surface output photon density was calculated by transforming the pixel gray levels in the CCD image into the light unit according to our experimentally established calibration formula. The optical parameters of each material were computed by an optical tomography procedure. Specifically, the specimen was modeled as a semi-infinite homogeneous medium. Steady-state diffusion theory was applied with the extrapolated boundary condition that the photon density was zero at an artificial boundary parallel to the boundary of the medium. Then, an analytic formula was used to predict the photon density on the bottom surface. Finally, a nonlinear least-square fitting was done to determine the absorption coefficient *μ*_{a}
and the reduced scattering coefficient *μ*_{a}
′ [11, 12, 13]. The calculated optical parameters of the four regions are given in Table 3. Fig. 11 shows the matching results between the experimental and computational output flux profiles. These optical parameters were then used as input to our finite-element based reconstruction method for BLT.

### 4.2.4 Experimental data acquisition

The heterogeneous mouse chest phantom containing the two light sources was placed on a sample holder in front of the CCD camera. The experimental setup is placed in a totally dark environment. The flux density was recorded with the CCD camera on the cylindrical surface of the phantom, along four radial directions separated by 90 degrees, as schematically shown in Fig. 10(b). During each data acquisition, one luminescent view was taken by exposing the camera for 60 seconds, as shown in Fig. 12. Furthermore, the recorded pixel gray levels of the luminescent view were transformed into corresponding light units according to the aforementioned calibration relationship.

### 4.2.5 Permissible source region

To regularize the BLT solution, a permissible source region was assigned by analyzing the four luminescent views taken by the CCD camera. These four planar images show high value clusters near the center of the front view image and a low value distribution for the back view image. In the right view and left view images, it can be seen that one half side exhibits high values while the other half shows low values. From these observations, we infer that the light source region should be in the anterior part of the phantom. Hence, the permissible source region should be in the left part of the phantom. Along the longitudinal direction, high signal-to-noise ratios are clustered between z=1.9mm to 28.1mm relative to the phantom bottom. Beyond the above region, the signal-to-noise ratios were insignificant and were ignored to reduce the computational burden. Finally, the permissible source region was defined as Ω_{s} ={(*x*,*y*,*z*)|*x*<0, 13.5<*z*<16.5, (*x*,*y*,*z*)∊ L} , which is the left blue-colored region in Fig. 13 (b).

### 4.2.6 Light source reconstruction

To simulate the photon propagation in the phantom, a geometrical model of diameter 30mm and height 26.2mm was established corresponding to a middle section of the physical phantom. Based on this model, a finite-element mesh was built consisting of 11340 wedge elements and 6576 nodes with 1024 datum nodes on the phantom surface, as shown in Fig. 13(a). The optical properties of every element were assigned in reference to the optical parameters reported in Section 4.2.3. On the surface of the geometric model, 16 circles, separated by about 1.75 mm, were selected, along each of which 64 detection locations were uniformly distributed. The measured photon density at each detector location was obtained from the CCD luminescent image using our calibration formula. The computed photon density at the corresponding detection point was obtained using Eq. (3) in Section 2. Then, the reconstruction method described in Section 3 was applied to reconstruct the light source distribution in the heterogeneous phantom. The reconstructed results correctly revealed that there were two strong light sources in the phantom located at (-8.6, 2.0, 15.0) with flux density 64.08 nano-Watts/mm^{3} and at (-8.6, -2.0, 15.0) with 54.04 nano-Watts/mm^{3}, respectively. The former was estimated to yield a total power of 94.2 nano-Watts (the total power = source volume × source flux density = 1.47 mm^{3}×64.08 nano-Watts/mm^{3} = 94.2 nano-Watts), while the latter was computed to have a total power of 79.4 nano-Watts (1.47×54.04 = 79.4 nano-Watts). Note that the volumes of the reconstructed sources are different from the actual source volumes, depending on the discretized element size. The smaller the element size, the higher the computational cost, and the closer to the actual source volume, as shown in Section 4.1.1. Fig. 13(b) shows the reconstructed source distribution. The differences between the reconstructed and real source positions were 0.72mm and 0.72mm for the two sources, respectively. The relative errors in the source strength were about 10.4% and 18.5%, respectively. The computed surface photon density based on the reconstructed light sources was in good agreement with the experimental counterpart, with the average relative error being about 13% as shown in Fig. 14.

## 5. Discussion and conclusions

We have developed a reconstruction algorithm to identify a 3D bioluminescent source distribution by incorporating *a priori* knowledge. Based on the finite-element discretization of the diffusion equation, a direct linear relationship has been established between the surface measurement and the underlying source distribution, and put within a Bayesian framework for a linear least square optimization with simple bound constraints [22]. Despite the ill-posed nature of this inverse source problem in the general case, it has been theoretically proven that the solution uniqueness in BLT can be established under practical constraints using *a priori* knowledge [15]. This reconstruction method has incorporated *a priori* knowledge, especially about the permissible source region, to enhance numerical stability and efficiency. The simulation and experiments have shown that the method is computationally efficient and fairly robust with respect to noise, initial distribution, and permissible region size. In addition, the finite-element based reconstruction method can handle a complex geometrical model, and it is suitable with small animals that have complicated anatomies.

The use of the permissible source region is a helpful technique for BLT. Our experiments have indicated that the smaller the permissible source region, the more stable the BLT reconstruction. An interesting observation is that when the permissible source region becomes larger and larger, distortion in the reconstructed source shape and power will be more and more significant. An effective remedy is a multi-scale BLT reconstruction. Initially, low resolution BLT will reliably indicate clusters of bioluminescent sources. Consequently, permissible regions can be re-defined to contain only the bioluminescent clusters. Iteratively, optimal results can eventually be obtained. Our physical phantom experiment has clearly demonstrated the success of BLT, with errors in terms of the source position and strength at about 1mm and less than 20%, respectively.

To perform BLT on small animals, especially mice, key issues must be resolved, including modeling of the individualized anatomy, determination of the optical parameters, and so on. An independent tomography tool(s), such as CT, MRI and/or diffuse optical tomography (DOT) scans can be instrumental in compensating for the heterogeneous structures of the mouse. We are currently improving the accuracy of geometric modeling of the mouse based on CT data and recovering optical parameters using optical means.

We reiterate that the emphasis here is on the feasibility of the presented reconstruction method and that the reported results are preliminary. To make the numerical simulation more realistic in this study, strong Gaussian noise (10%~15%) was added to the data generated by the diffusion approximation based forward model. More importantly, the data in our physical phantom experiment were produced by physical sources, and hence they are totally free of the well-known “*inverse crime*” [25]. Furthermore, while diffuse optical tomography (DOT) recovers tissue optical parameters and uses nonlinear algorithms, BLT can be used to recover source features and can be formulated as a linear inverse problem [15], which may be advantageous relative to DOT in terms of reconstruction quality. More importantly, we utilized a permissible source region based on *a priori* knowledge to enhance the numerical stability and accuracy of the BLT reconstruction. As a side note, in our numerical and physical phantom studies, the reconstruction errors were found to be about 14% and about 19% of the true source power, respectively. That is, the error difference in the numerical simulation and the phantom experiment is only about 5%. We believe that various systematic and measurement errors in the phantom experiment may have been canceled out to some degree, and the *a priori* knowledge (the permissible source region, *etc*.) must have helped generate good results in this feasibility project. Nevertheless, it is clear that a comprehensive error analysis needs to be performed, involving measurement bias, model mismatch, optical characteristic inaccuracy, geometric errors, data processing related complications, and so on. Relevant data will be collected and reported in follow-up publications.

In conclusion, we have developed a finite-element based reconstruction method for BLT and demonstrated its feasibility in numerical and phantom experiments. The initial results are very encouraging and suggest that BLT has great potential for advancing current planar bioluminescent imaging techniques into a 3D quantitative modality for molecular imaging. Further research activities in this area are being actively performed, including a comprehensive error analysis project.

## Acknowledgments

This work is supported by an NIH/NIBIB grant EB001685, a special grant for bioluminescent imaging from Department of Radiology, College of Medicine, University of Iowa. Dr. Ming Jiang was also supported in part by NKBRSF (2003CB716101) and NSFC (60325101, 60272018 and 60372024), China.

## References and links

**1. **S. Bhaumik and S. S. Gambhir, “Optical imaging of Renilla luciferase reporter gene expression in living mice,” Proc. Natl. Acad. Sci. USA **99**, 377–382 (2002). [CrossRef]

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

**3. **P. Ray, A.M. Wu, and S.S. Gambhir, “Optical bioluminescence and positron emission tomography imaging of a novel fusion reporter gene in tumor xenografts of living mice,” Cancer Res. **63**, 1160–1165 (2003). [PubMed]

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

**5. **J. Welch and M. J. C. van Gemert, Optical and Thermal response of laser-irradiated tissue (Plenum Press, New York,1995).

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

**7. **V. Ntziachristos, C. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. **8**, 757–760 (2002). [CrossRef] [PubMed]

**8. **R. Schultz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with non-contact measurements,” IEEE Trans. Med. Imag. **23**, 492–500 (2004). [CrossRef]

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

**10. **J. Ripoll, D. Yessayan, G. Zacharakis, and V. Ntziachristos, “Experimental determination of photon propagation in highly absorbing and scattering media,” J. Opt. Soc. Am. A **22**, 546–551 (2005). [CrossRef]

**11. **T. J. Farrell, M. S. Patterson, and B. Wilson, “A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo,” Med. Phys. **19**, 879–888 (1992). [CrossRef] [PubMed]

**12. **M. Gurfinkel, T. S. Pan, and E. M. Sevick-Muraca, “Determination of optical properties in semi-infinite turbid media using imaging measurements of frequency-domain photon migration obtained with an intensified charge-coupled device,” J. Biomed. Opt. **9**, 1336–1346 (2004). [CrossRef] [PubMed]

**13. **M. Guven, B. Yazici, X. Intes, and B. Chance, “Diffuse optical tomography with *a priori* anatomical information,” Phys. Med. Biol. **50**, 2837–2858 (2005). [CrossRef] [PubMed]

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

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

**16. **W. Cong, D. Kumar, Y. Liu, A. Cong, and G. Wang, “A practical method to determine the light source distribution in bioluminescent imaging,” Proc. SPIE **5535**, 679–686 (2004). [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]

**18. **J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor analysis (
Wiley, New York,1976).

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

**20. **S. S. Rao, The finite element method in engineering (Butterworth-Heinemann, Boston,1999).

**21. **S. C. Brenner and R. L. Scott, The Mathematical Theory of Finite Elements (Springer, Berlin-Heidelberg-New York,1994).

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

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

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

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