Bioluminescent imaging has the capability to reveal molecular and cellular activities directly, and can be applied to all disease processes in most small-animal models. It is very sensitive, and helps diagnose diseases, monitor therapies, and facilitate drug development.1, 2 Bioluminescent imaging has been in a planar mode and largely a qualitative imaging tool.3 With bioluminescence tomography (BLT), 3-D localization and quantification is enabled of a bioluminescence source distribution inside a living small animal such as a mouse.4, 5
In bioluminescence imaging, photon scattering predominates over absorption in the biological tissue in this spectral range of interest. As a result, 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.3 In this case, the photons propagation in the biological tissue can be well described by steady-state diffusion equation and Robin-type boundary conditions.6, 7 The BLT principles and solution uniqueness for the BLT problem conditions were studied under some practical conditions. 4, 5, 8, 9 Some numerical algorithms for BLT were already reported, including the finite element based BLT algorithms, 10, 11, 12, 13, 14 BLT in combination with PET (OPET),15 multispectral bioluminescence optical tomography,16, 17 and BLT method based on diffusion theory of half infinite medium.18 In this letter, for the first time we develop a BLT algorithm using the boundary integral method. This methodology only requires finite element meshing of structural boundaries instead of complex volumetric finite elements previously used for BLT. Hence, the complexity and stability of the new BLT algorithm are improved as compared to that of the finite-element-based BLT algorithms.
The overall support region of a small animal can be decomposed into a number of subregions with ; for example, lungs, heart, liver, bone, muscle, and so on. Applying the Gauss theorem, the steady-state diffusion equation can be transformed to a boundary integral equation on each subregion19:is a Green function of the steady-state diffusion equation, is the diffusion coefficient given by , and where is the absorption coefficient and is the reduced scattering coefficient for the subregion . Equation 1 is a well-posed second kind integral equation. The bounding surface can then be split into surface elements on which the function and are approximated by use of a set of interpolation points and interpolation functions on 19: and represent the values of the function and at an interpolation point , respectively. Since the BLT reconstruction is underdetermined and ill-posed, we incorporate a priori knowledge to define a permissible source region for a subregion without loss of generality, where the bioluminescence source may be distributed. The region can be discretized, and the embedded source function is approximated on as represent the value of the source function at an interpolation point is the interpolation basis function, and is the number of discrete values of the source function. , , is the number of nodal points for , , and is a strictly diagonally dominant matrix, and allows its inversion. Multiplying Eq. 4 with the inverse , we have and at the interface and the Robin-type external boundary condition, the above matrix equation for can be assembled into19 denotes the source distribution in the permissible source region, and a vector consisting of photon density values at the boundary and interface nodes, which can be divided into at internal interface nodes and at exterior boundary nodes. Because is still a diagonally dominant matrix and invertible, we have . After removal of those rows of that correspond to we obtain and a linear relationship between measurable photon density at boundary nodes and the source distribution directly from Eq.7. Instead, an optimization procedure is employed to find a solution by minimizing the following objective function: is the upper bound that is chosen to be physically meaningful, a stabilizing function, the regularization parameter, the weight matrix, and .
To demonstrate the feasibility and efficiency of our new algorithm, we carried out an experiment using a heterogeneous mouse chest phantom. The physical phantom of 30-mm height and 30-mm diameter was fabricated. It consisted of four types of high-density polyethylene materials: (8624K16), nylon 6/6 (8538K23), delrin (8579K21), and polypropylene (8658K11) (McMaster-Carr Supply Company, Chicago, Illinois), to represent muscle (M), lungs (L), heart (H), and bone (B) respectively. Based on the diffuse model of photon propagation, the optical parameters of the four materials at wavelength of about were independently found from transillumination experimental data,12 which are listed in Table 1 . The luminescent light stick (Glowproducts, Victoria, British Columbia, Canada) was selected as the testing source. The stick consisted of a small glass vial containing one chemical solution and a large 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 and emit red light around , which is spectrally similar to bioluminescent light generated by the luciferase. Two small holes of diameter and height were longitudinally drilled in the left lung region of the phantom with their centers at and , respectively, as shown in Fig. 1 . Two catheter tubes about in height were filled with red luminescent liquid, and were placed inside the two holes as light sources, respectively. Emitting light power of two catheter tubes was measured with the CCD camera. They were 105.1 nW and , respectively.
Optical parameters of the heterogeneous mouse chest phantom.
|Material||Muscle (M)||Lung (L)||Heart (H)||Bone (B)|
The experiment was conducted in a totally dark environment to avoid background noise. The flux density on the cylindrical surface of the phantom was recorded with the sensitive CCD camera, along four radial directions separated by , as schematically shown in Fig. 1b. During each data acquisition session, one luminescent view was taken by exposing the camera for , as shown in Fig. 2
A permissible source region was assigned as to regularize the BLT solution. Along the longitudinal direction high signal-to-noise ratios (SNR) were clustered between to relative to the phantom bottom. Beyond the above region, the SNR were insignificant, and ignored to reduce the computational cost. To simulate the photon propagation in the phantom, a geometrical model of diameter and height was constructed corresponding to a middle section of the physical phantom. Based on this model, a discrete mesh was generated for the external boundary and internal interface. The external boundary mesh consisted of 960 quadrilateral elements and 1024 measurement datum points. The internal interface mesh consisted of 1680 quadrilateral elements and 1792 nodes. The permissible source region was also discretized into 308 wedge elements, as shown in Fig. 3b. The measured photon density at each detector location was computed from the CCD luminescent image using our calibration formula. Then, the proposed algorithm 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 with flux density and at with , respectively. The former was estimated to yield a total power of (the total flux , while the latter was computed to have a total power of . Figures 3a and 3b shows the reconstructed source distribution. The differences between the reconstructed and real source positions were 1.9 and for the two sources, respectively. The relative errors in the source strength were 3.0% and 15.3%, respectively. The computed surface photon density profiles based on the reconstructed light sources were in good agreement with the experimental counterparts with an average relative error about 20%.
In conclusion, we have developed a boundary integral method to reconstruct a 3-D bioluminescent source distribution embedded in a heterogeneous object of interest from measured flux data on the external surface of the object. The feasibility of the new reconstruction method has been demonstrated in a physical heterogeneous phantom experiment. The results have indicated that the method can reliably localize the source locations, and recover the source strength fairly well. The reconstruction results may be further improved by decreasing the boundary element size, lowering the measurement noise, and increasing the accuracy of optical parameters. The method only requires the discretization on the object boundary and structural interface, which is much easier than meshing for volumetric finite elements. Hence, the new method can handle a complex geometrical model much efficiently than finite element based BLT algorithm.
This work is supported by an NIH/NIBIB Grant EB001685 and the U.S. Army's Breast Cancer Research program Concept Award W81XWH-05–1–0461.