1. Introduction

Bioluminescence imaging (BLI) is an emerging technique that can be used to monitor molecular events in intact living animals [1–3]. Important applications of this imaging technique include gene therapy and cell traffic studies. Unlike fluorescence approaches which require an external source of light for excitation of fluorophores, BLI generates a two-dimensional (2D) view of gene expression based on the internal light produced by luciferases, cataysts in a light generating reaction, through the oxidation of an enzyme-specific substrate (luciferin) [2]. BLI, however, has limited resolution due to the nature of surface imaging. In addition, BLI is not able to obtain depth information of luciferase expression.

In this paper, model-based bioluminescence tomography is explored in order to extract the depth information and enhance resolution of bioluminescence imaging. A finite element-based reconstruction algorithm is implemented for recovering the three-dimensional spatial map of gene expression. Both the simulated and experimental results obtained show that an arbitrary internal source distribution can be reconstructed using our three-dimensional finite element algorithm.

2. Methods and materials

2.1 Reconstruction algorithm

Bioluminescent light propagation in continuous-wave mode in tissue can be described by the well-known diffusion equation as follows [4]

where Φ(x, y, z) is the photon density; D and µ_a are the diffusion and absorption coefficients, respectively; S(x, y, z) is the source distribution of gene expression.

The following Type III boundary conditions (BCs) are used for the solution of Eq. (1) [5,6]

where î_n is the unit normal vector for the boundary surface and α is a coefficient related to the internal reflection at the boundary.

A set of matrix equations capable of inverse problem solution can be obtained through the finite element discretization of Eq. (1) and its differentiation and regularized Newton’s method [7,8]

where the elements of matrix [A] and vector {b} are, respectively, a_ij=〈-D∇ϕ_j·∇ϕ_i-µ_aϕ_jϕ_i〉 and b_i=〈${\displaystyle \sum _{\mathrm{k}}^{\mathrm{N}}}$S_kϕ_kϕ_i〉 where 〈 〉 indicates integration over the problem domain and N is the total number of nodes in the finite element mesh used; ϕ_i, ϕ_j and ϕ_k are locally spatially varying Lagrangian basis functions at nodes i, j and k, respectively; ℑ is the Jacobian matrix that should be formed by ∂Φ/∂S at the boundary measurement sites; ΔS=(ΔS₁,ΔS₂,⋯,ΔS_N)^T is the update vector for the light source profile; Φ^(m)=(${\mathrm{\Phi }}_1^{\left(\mathrm{m}\right)}$,${\mathrm{\Phi }}_2^{\left(\mathrm{m}\right)}$,⋯,${\mathrm{\Phi }}_{\mathrm{M}}^{\left(\mathrm{m}\right)}$) and Φ^(c)=(${\mathrm{\Phi }}_1^{\left(\mathrm{c}\right)}$,${\mathrm{\Phi }}_2^{\left(\mathrm{c}\right)}$,⋯${\mathrm{\Phi }}_{\mathrm{M}}^{\left(\mathrm{c}\right)}$) where ${\mathrm{\Phi }}_{\mathrm{i}}^{\left(\mathrm{m}\right)}$ and ${\mathrm{\Phi }}_{\mathrm{i}}^{\left(\mathrm{c}\right)}$, respectively, are measured and calculated data for i=1, 2,…, M boundary locations. Note that to estimate S(x,y,z) spatially we have expanded this quantity as a finite sum of unknown coefficient multiplied by the locally defined Lagrangian basis function.

In bioluminescence tomography, the goal is to iteratively update the S(x,y,z) distribution through the solution of Eqs. (3)–(5) so that a weighted sum of the squared difference between computed and measured optical data can be minimized. While the spatial variation of diffusion and absorption coefficients needs to be considered ultimately in bioluminescence tomography, in this study we have assumed these two parameters are constant during the reconstruction of S(x,y,z).

2.2 Simulations

The reconstruction algorithm was first tested using simulated data. The imaging geometry tested was a 30×30×20mm cube. The diffusion and absorption coefficients of this background medium were set to 0.33 mm and 0.005 mm^-1, respectively. One or two approximate cylindrical targets (7.0mm in diameter and 10mm in height) were embedded 7.2mm off-center in the background to simulate the gene expression. In the target volumes, the source coefficient (S_k) was 1.0, whereas it was zero in the background volume. A 3D finite-element mesh with 1749 nodes and 7440 tetrahedral elements was used for the reconstruction results presented. All the reconstructions were results of 50 iterations, after which no noticeable improvement was observed. The computations were conducted in a 2.6 GHz Pentium IV PC.

2.3 Experiments

In our experiments, Quantilum^® recombinant luciferase (Promega, Milwaukee, WI), a standard firefly luciferase with a peak emission at 560nm, was used as the source of enzyme, while a luciferase assay system (Promega, Milwaukee, WI) provided luciferin, magnesium, ATP and other co-factors for generating high total light output with a half life of 10 minutes. Quantilum^® recombinant luciferase was diluted 10⁶-fold in 1 x cell culture lysis reagent (supplied by the luciferase assay system) and 1mg/ml acetylated BSA (Sigma-Aldrich, St Louis, MO). The luciferase assay substrate was dissolved into luciferase assay buffer to obtain luciferase assay reagent. Light was generated when combining the diluted luciferase (350µl) with the luciferase assay reagent (350µl).

The luciferase-luciferin reagent prepared was then placed in 9.8×10mm cylindrical hole embedded at a depth of 7.7 mm from the nearest surface in a 30×30×27.5mm solid cubic phantom (1% Intralipid solution +1% Agar powder). The cubic phantom had D=0.33 mm and µ_a=0.0023 mm^-1. In the experiments, the luciferase-luciferin containing cubic phantom was mounted on a 360° rotation stage. By rotating the phantom three 90°, luminescence light emitted from each side of the phantom was collected by a thermoelectrically cooled (-30°C), front illuminated CCD array (Roper Scientific Instrumentation, Trenton, NJ) coupled with an optical lens subsystem (Zoom 7000, Navitar, Rochester, NY). For each side measurement, an image of 100×90 pixels was recorded with an exposure time of 60s. We used only 480 measurement data for 3D image reconstruction in which the finite element mesh had 1668 nodes and 7788 elements.

3. Results and discussion

Let us now examine the simulated and experimental results. Figures 1(a) and (b) show the exact and reconstructed images of source distribution for a single target, while Figs. 1(d) and (e) display the exact and reconstructed images of source for two targets (a typical cross sectional image for both cases are shown in Figs. 1(c) and 1f, respectively). We can see that the targets are clearly detected for each case. We also note that the background is reconstructed accurately without any boundary artifacts. The recovered image of a single target from the experimental data is presented in Fig. 2 where the exact image and a selected cross sectional slice of the reconstructed image are also shown. We immediately note that the target is successfully reconstructed in terms of its location, shape and source strength. The reconstructed background is quite smooth with slight boundary artifacts.

 

Fig 1. Source distribution (S(x,y,z)) images reconstructed from 748 simulated measurement data: (a) Exact image of the single target. (b) Reconstructed image of the single target. (c) A cross sectional image of the single target at X=10mm. (d) Exact image of the two targets. (e) Reconstructed image of the two targets. (f) A cross sectional image at Z=10mm of the two targets. The axes illustrate the spatial scale, in

Download Full Size | PDF

 

Fig. 2. Source distribution (S(x,y,z)) image of a single target reconstructed from 480 experimental data measured from the luciferase-luciferin containing phantom. (a) Exact image. (b) Reconstructed image. (c) A cross sectional image at X=10mm. The axes illustrate the spatial scale, in millimeters, whereas the color bar records the source coefficient (arbitrary units).

Download Full Size | PDF

Quantitative analysis of the experimental image reveals that the target size is somewhat underestimated particularly along the transverse plane (relative to the four data collection planes). In addition, we currently cannot detect a target at the center of the cubic phantom (i.e., at a depth of 1.5cm from the surface). We believe these limitations are largely due to the use of the front illuminated CCD camera at -30°C which provided a relative low signal-to-noise ratio (SNR). As suggested by Rice et al. [9], a back illuminated CCD should give at least 250% better SNR than a front illuminated CCD. Once such a back illuminated CCD is used in BLT, we expect higher quality of image reconstruction can be obtained. Another limiting factor of the current BLT is the effect of background optical heterogeneity, i.e., the impact of spatial distribution of optical properties. This, however, may be overcome by combining BLT with diffuse optical tomography (DOT) which can provide accurate spatial distribution of optical properties; hence background optical heterogeneity could be known as a priori for BLT reconstruction. Finally, the ill-posed nature involved in BLT remains challenging which would require continued effort in searching more effective methods such as advanced regularization techniques.

In summary, a model-based reconstruction algorithm is implemented for bioluminescence tomography. The reconstruction approach presented promises to provide the depth information of gene expression as well as improve the resolution of bioluminescence imaging. In vivo experimental evaluation of the 3D reconstruction algorithm described is underway in our laboratory.