Practical reconstruction method for bioluminescence tomography

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 OCIS codes: (110.6960) Tomography; (170.3010) Image reconstruction techniques; (170.6280) Spectroscopy, fluorescence and luminescence. ___________________________________________________________________________________________ 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). 2. C. Contag and M. H. Bachmann, “Advances in Bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. 4, 235-260 (2002). 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). 4. W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light-emitting probes,” J. Biomed. Opt. 6, 432440 (2001). 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). 7. V. Ntziachristos, C. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8, 757-760 (2002). 8. R. Schultz, J. Ripoll, and V. Ntziachristos, “Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imag. 23, 492-500 (2004). 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, 299309 (1993). 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). (C) 2005 OSA 5 September 2005 / Vol. 13, No. 18 / OPTICS EXPRESS 6756 #8070 $15.00 USD Received 11 July 2005; revised 16 August 2005; accepted 18 August 2005 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). 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). 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). 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). 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). 17. X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with modelbased reconstruction,” Opt. Express 12, 3996-4000 (2004). 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, 17791792 (1995). 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-HeidelbergNew York, 1994). 22. J. C. Ye, K. J. Webb, C. A. Bouman, and R. P. Millane, “Optical diffusion tomography by iterativecoordinate-descent optimization in a Bayesian framework,” J. Opt. Soc. Am. A 16, 2400-2412 (1999). 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). ___________________________________________________________________________


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.

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; ( ) represents the photon density [Watts/mm 2 ] at location x; ( ) S x the bioluminescent source density [Watts/mm 3 ]; ( ) ( ) μ the scattering coefficient [mm -1 ]; and g the anisotropy parameter.

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 ν 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

Finite-element discretization
The governing equations ( 1)-( 2) for ( ) can be equivalently represented as the following weak form [9] that ( ) for an arbitrary test function ( ) . Following the standard finite-element method [20], Ω is discretized with T vertex nodes ( ) , ,..., T N N N and e N elements, denoted as , such that ( ) can be approximated with a piecewise polynomial function [9,[19][20][21], where k φ is the nodal value of ( ) ( ) ( ) ( ) where k S are the values of ( ) 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.

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

Since the matrix [ ]
M is symmetric and positive definite, m Φ can be uniquely determined from ( )( ) assuming that the source p S 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 p S from measured meas Φ .Generally, since the measured data in bioluminescent imaging are corrupted by noise, it is not practical to directly solve for p S 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 meas k Φ is the measured photon density at the k-th detector position ( ) . We assume that the noise is independent for different detectors so that the covariance matrix C of data meas Φ is diagonal and expressed by To simplify the notation, we define the diagonal matrix W as 1 The data likelihood is then formulated as a function of multi-variables: where 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 ( ) 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).

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.

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 ( ) 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).

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 ( ) 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

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 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  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 8inches 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

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.

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 s μ′ [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.

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.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

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.08nano-Watts/mm 3 = 94.2 nano-Watts), while the latter was computed to have a total power of 79.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.

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.

Fig. 1 .
Fig. 1.Numerical simulation for BLT reconstruction of one source.(a) The true source distribution in the left lung consisting of 6 volume elements and having a homogeneous density of 200.0pico-Watts/mm 3 , and (b) the counterpart reconstructed from the surface data corrupted by 10% Gaussian noise.The average density error is 0.5%.

Fig. 2 .Fig. 3 .Fig. 4 .
Fig. 2. Numerical simulation for BLT reconstruction of two sources in the left and right lungs, respectively.(a) The true source distribution with a density of 200.0pico-Watts/mm 3 for each source, and (b) the counterparts reconstructed from the surface data corrupted by 10% Gaussian noise.The average density error is 3%.

Fig. 5 .
Fig. 5. Numerical simulation for BLT reconstruction of four sources.(a) Four bioluminescent source separations of 1mm between first and second, 2mm between second and third, and 3mm between third and fourth source with the density of 200.0 pico-Watts/mm 3 for each source, and (b) the counterparts reconstructed from the surface data corrupted by 10% Gaussian noise.The source positions are accurately identified with their density being recovered to 196.1 pico-Watts/mm 3 , 184.7 pico-Watts/mm 3 , 175.1 pico-Watts/mm 3 and 181.5 pico-Watts/mm 3 , respectively.

1 s Ω , 2 s Ω and 3 s
are 154, 308 and 462 volume elements in these three permissible regions, respectively.The bioluminescent source distribution was reconstructed in these three permissible regions Ω , for different noise levels, respectively.The BLT results associated with the different permissible region sizes and data noise levels are summarized in Table

Fig. 10 .
Fig. 10.Mouse Chest phantom.(a) A heterogeneous mouse phantom consisting of bone (B), heart (H), lungs (L), and muscle (M); (b) a middle cross-section through two hollow cylinders for hosting luminescent sources in the left lung.The four arrows show the direction of the CCD camera during data acquisition.
4 nano-Watts (1.47×54.04= 79.4nano-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.

Fig. 14 .
Fig. 14.Comparison between measured and computational photon density profiles along the detection circle on the phantom surface at heights (a) 10.6mm,(b) 15.9mm, and (c) 21.1mm, from the top surface of the model.

Table 1 .
Optical parameters for the numerical phantom

Table 2 .
Relative error (%) with BLT results of total source power 14.