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.³ 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.³ 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),¹⁵ multispectral bioluminescence optical tomography,^{16, 17} and BLT method based on diffusion theory of half infinite medium.¹⁸ 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 $\Omega$ can be decomposed into a number of subregions ${\Omega }_j$ $(j=1,2,\dots ,\tau )$ with $\Omega ={\sum }_{j=1}^{\tau }{\Omega }_j$ ; 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 subregion¹⁹:

where $G(\mathbf{r},\mathbf{x})=\exp (-{\mu }_{eff}\mid \mathbf{r}-\mathbf{x}\mid )∕(4\pi D\mid \mathbf{r}-\mathbf{x}\mid )$ is a Green function of the steady-state diffusion equation, $D$ is the diffusion coefficient given by $D=1∕\left[3\right({\mu }_a+{\mu }_s^{\prime }\left)\right]$ , and ${\mu }_{eff}=\left[3{\mu }_a\right({\mu }_a+{\mu }_s^{\prime }){]}^{1∕2}$ where ${\mu }_a$ is the absorption coefficient $\left({\mathrm{mm}}^{-1}\right)$ and ${\mu }_s^{\prime }$ is the reduced scattering coefficient $\left({\mathrm{mm}}^{-1}\right)$ for the subregion ${\Omega }_j$ . Equation 1 is a well-posed second kind integral equation. The bounding surface $\partial {\Omega }_j$ can then be split into $N$ surface elements ${\Gamma }_i^j$ $(i=1,2,\dots ,N)$ on which the function $\Phi \left(\mathbf{x}\right)$ and $\partial \Phi \left(\mathbf{x}\right)∕\partial n$ are approximated by use of a set of $p$ interpolation points and interpolation functions ${\varphi }_k\left(\mathbf{x}\right)$ on ${\Gamma }_i^j$ ¹⁹:The shape of the surface element can be arbitrarily selected, and usually made as a quadrilateral or a triangle. The number of points per surface element depends on the accuracy needed in the interpolation procedure. Let ${\Phi }_k=\Phi \left({\mathbf{x}}_k\right)$ and $q_k=D\partial \Phi \left({\mathbf{x}}_k\right)∕\partial n$ represent the values of the function $\Phi \left(\mathbf{x}\right)$ and $D\partial \Phi \left(\mathbf{x}\right)∕\partial n$ at an interpolation point ${\mathbf{x}}_k$ , respectively. Since the BLT reconstruction is underdetermined and ill-posed, we incorporate a priori knowledge to define a permissible source region ${\Omega }_j^s({\Omega }_j^s\subseteq {\Omega }_j)$ for a subregion ${\Omega }_j$ without loss of generality, where the bioluminescence source may be distributed. The region ${\Omega }_j^s$ can be discretized, and the embedded source function $S\left(\mathbf{x}\right)$ is approximated on ${\Omega }_j^s$ aswhere $s_k=S\left({\mathbf{x}}_k\right)$ represent the value of the source function $S\left(\mathbf{x}\right)$ at an interpolation point ${\mathbf{x}}_k,{\psi }_k$ $\left(\mathbf{x}\right)$ is the interpolation basis function, and $N_s$ is the number of discrete values of the source function.

Inserting Eqs. 2, 3 into Eq. 1, we obtain the following matrix equation:

where $\Phi =({\Phi }_1,{\Phi }_2,\dots ,{\Phi }_M{)}^T$ , $\mathrm{Q}=(q_1,q_2,\dots ,q_M{)}^T$ , $M$ is the number of nodal points for $\partial {\Omega }_j$ , $S=(s_1,s_2,\dots ,s_{N_s}{)}^T$ , and ${\mathbf{H}}_j$ is a strictly diagonally dominant matrix, and allows its inversion. Multiplying Eq. 4 with the inverse ${\mathbf{H}}_j^{-1}$ , we haveFurthermore, in terms of the continuity of $\Phi \left(\mathbf{x}\right)$ and $D\partial \Phi \left(\mathbf{x}\right)∕\partial n$ at the interface and the Robin-type external boundary condition, the above matrix equation for $j=1,2,\dots ,\tau$ can be assembled into¹⁹ where $S$ denotes the source distribution in the permissible source region, and $\Phi$ a vector consisting of photon density values at the boundary and interface nodes, which can be divided into ${\Phi }^{\mathrm{in}}$ at internal interface nodes and ${\Phi }^{\mathrm{ex}}$ at exterior boundary nodes. Because $\mathbf{M}$ is still a diagonally dominant matrix and invertible, we have $\mathbf{B}={\mathbf{M}}^{-1}\mathbf{F}$ . After removal of those rows of $\mathbf{B}$ that correspond to ${\Phi }^{\mathrm{in}}$ we obtain ${\mathbf{B}}^{ex}$ and a linear relationship between measurable photon density at boundary nodes and the source distributionGenerally, the measured data in bioluminescence imaging are corrupted by noise, so it is not practical to solve for $S$ directly from Eq.7. Instead, an optimization procedure is employed to find a solution by minimizing the following objective function:where $u$ is the upper bound that is chosen to be physically meaningful, $\eta \left(S\right)$ a stabilizing function, $\alpha$ the regularization parameter, $\mathbf{W}$ the weight matrix, and $\parallel V{\parallel }_W=V^{\mathrm{T}}\mathbf{W}V$ .

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 $650\mathrm{nm}$ were independently found from transillumination experimental data,¹² 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 $650\mathrm{nm}$ , which is spectrally similar to bioluminescent light generated by the luciferase. Two small holes of diameter $0.6\mathrm{mm}$ and height $3\mathrm{mm}$ were longitudinally drilled in the left lung region of the phantom with their centers at $(-9.0,1.5,15.0)$ and $(-9.0,-1.5,15.0)$ , respectively, as shown in Fig. 1 . Two catheter tubes about $1.9\mathrm{mm}$ 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 $97.4\mathrm{nW}$ , respectively.

Table 1

Optical parameters of the heterogeneous mouse chest phantom.

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 $90\mathrm{deg}$ , as schematically shown in Fig. 1b. During each data acquisition session, one luminescent view was taken by exposing the camera for $60\mathrm{s}$ , as shown in Fig. 2

A permissible source region was assigned as ${\Omega }_s=\left\{\right(x,y,z)\mid x<0,13.0<z<17.0,(x,y,z)∊L\}$ to regularize the BLT solution. Along the longitudinal direction high signal-to-noise ratios (SNR) were clustered between $z=1.9$ to $28.1\mathrm{mm}$ 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 $30\mathrm{mm}$ and height $26.2\mathrm{mm}$ 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 ${\Omega }_s$ 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 $(-7.5,2.7,15.0)$ with flux density $73.6\mathrm{nW}∕{\mathrm{mm}}^3$ and at $(-7.5,-2.7,15.0)$ with $56.1\mathrm{nW}∕{\mathrm{mm}}^3$ , respectively. The former was estimated to yield a total power of $108.2\mathrm{nW}$ (the total $\mathrm{power}=\mathrm{source}$ $\mathrm{volume}\times \mathrm{source}$ flux $\mathrm{density}=1.47{\mathrm{mm}}^3\times 73.6\mathrm{nW}∕{\mathrm{mm}}^3=108.2\mathrm{nW})$ , while the latter was computed to have a total power of $82.5\mathrm{nW}$ $(1.47\times 56.1=82.5\mathrm{nW})$ . Figures 3a and 3b shows the reconstructed source distribution. The differences between the reconstructed and real source positions were 1.9 and $1.9\mathrm{mm}$ 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.