Open Access
1 March 2007 Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging
Author Affiliations +
Abstract
A new method is described for obtaining a 3-D reconstruction of a bioluminescent light source distribution inside a living animal subject, from multispectral images of the surface light emission acquired on charge-coupled device (CCD) camera. The method uses the 3-D surface topography of the animal, which is obtained from a structured light illumination technique. The forward model of photon transport is based on the diffusion approximation in homogeneous tissue with a local planar boundary approximation for each mesh element, allowing rapid calculation of the forward Green's function kernel. Absorption and scattering properties of tissue are measured a priori as input to the algorithm. By using multispectral images, 3-D reconstructions of luminescent sources can be derived from images acquired from only a single view. As a demonstration, the reconstruction technique is applied to determine the location and brightness of a source embedded in a homogeneous phantom subject in the shape of a mouse. The technique is then evaluated with real mouse models in which calibrated sources are implanted at known locations within living tissue. Finally, reconstructions are demonstrated in a PC3M-luc (prostate tumor line) metastatic tumor model in nude mice.

1.

Introduction

Bioluminescent imaging is a noninvasive technique for performing in vivo diagnostic studies on animal subjects in the areas of medical research, pathology, and drug discovery and development. Bioluminescence is typically produced by cells that have been transfected with a luminescent reporter such as luciferase, and can be used as a marker to differentiate a specific tissue type (e.g., a tumor), monitor physiological function, or follow the progression of a disease. A wide range of applications has been demonstrated, including areas of oncology,1, 2 infectious disease,3 inflammation, and metabolic disease.4, 5

Photons emitted by bioluminescent cells are strongly scattered in the tissue of the subject such that propagation is diffusive in nature.6 As photons diffuse through tissue, many are absorbed, but a fraction reach the surface of the subject and can be detected. In general, absorption in mammalian tissues is high in the blue-green part of the spectrum (<565nm) and lower in the red and near-infrared (NIR) part of the spectrum (600to900nm) .6 Firefly luciferase has a rather broad emission spectrum ranging from 500to700nm , of which part of the emission spectrum is in the low-absorption region.7, 8

Since the mean free path for scattering in tissue is short, on the order of 0.5mm , photons from deep sources are scattered many times before reaching the surface. Bioluminescent imaging systems effectively record the spatial distribution of these scattered photons emitted from the surface of the subject. However, the most important quantitative information is not directly related to the surface emission, but rather pertains to the bioluminescent source inside the subject. Important parameters are the source brightness (related to the number of light-emitting cells), position, and geometry. Most of the in vivo bioluminescent imaging work published to date involves use of single-view 2-D imaging systems,1, 2, 3, 4, 5, 7 for which image analysis involves quantifying a light-emitting region of interest (ROI) on the subject surface. While this analysis methodology is simple and provides a good relative measure of the internal light source brightness, it does not take into account the source depth and resulting attenuation through tissue necessary to obtain an absolute source brightness measurement. Hence, there is interest in developing imaging systems and algorithms that would provide the 3-D tomographic reconstruction of the distribution of photon emission inside the animal from images of radiance measured on the animal surface.

Other groups have reported on tomographic techniques for bioluminescent source location and power, considering tissue optical properties at a single wavelength. Wang, Li, and Jiang discuss issues of uniqueness of the tomographic problem.9 Gu present an image reconstruction technique for bioluminescent sources using a finite-element approach to compute the forward model in a cubic homogeneous phantom.10 Cong demonstrate the performance of a finite-element-based bioluminescence tomography algorithm in a cylindrical heterogeneous tissue phantom, which requires a priori information on a feasible source distribution.11 This group expanded their research to simplify and linearize the forward problem by reducing the finite element mesh resolution to large-scale organ structures,12 and applying the Born theory.13 Tomographic algorithms for small animal fluorescence imaging in a slab geometry have been developed by Ntziachristos 14, 15, 16

The diffusive nature of light in tissue does pose a challenge for determining unique solutions of the tomographic problem. With bioluminescent imaging, uniqueness can be improved by taking images of the subject from multiple views,17, 18 but this can add to the complexity of instrumentation and the length of imaging time. We propose a simpler and faster method that uses spectral information in the image data to constrain the tomographic solution. Bioluminescent images acquired through a number of bandpass filters were analyzed by Coquoz for single source depth and flux, in which the tissue surface is assumed to be planar.19 Examining phantom subjects with cylindrically shaped boundaries, Dehghani present reconstructions of bioluminescent sources from multispectral data.20 Chaudhari discuss their efforts to extract 3-D source distribution from multispectral bioluminescent images of mice by incorporating computed tomography (CT) data to determine the tissue surface topography with finite-element analysis to predict the photon transport in tissue.21

In this study, a fast (<1min) tomographic analysis of single-view multispectral images of bioluminescence is presented. The diffuse luminescence tomographic algorithm (DLIT™) provides estimates on the 3-D location and the photon flux of the sources in the tissue. We show that the kernel matrix can be computed efficiently for complex-shaped objects such as small animals, and is initiated naïve to a source distribution. Unique to our work, the surface topography is measured for each animal subject and is used in the calculation of the kernel matrix.

The intent of this work is to outline the experimental and computational procedures of the diffuse tomography reconstruction technique. Section 2 describes the theory and approximations used in our forward model of photon propagation, including the additional information offered by multispectral images of bioluminescence from animal subjects. The experimental setup for collecting bioluminescent images, as well as the characterization of the surface topography for boundary treatment, is presented in Sec. 3. The tomographic reconstruction technique is validated by analyzing images from a homogeneous phantom mouse in Sec. 4. Results are also presented for luminescent beads implanted in living mice and a metastatic mouse model with luciferase-labeled prostate tumor cells.

2.

Theory

2.1.

Photon Diffusion in Tissue

Light transport in turbid media such as tissue is dominated by scattering and is essentially diffusive in nature.22 The condition for diffusive transport is that the isotropic scattering coefficient μs be much greater than the absorption coefficient μa , so that the change in the photon density is very small between scattering events. In this case, the photon density ρ(x) in a homogeneous medium is governed by the steady-state diffusion equation:

Eq. 1

D2ρ(x)μacρ(x)=U(x),
where U is the photon rate density (photonsscm3) and D is the diffusion coefficient,
D=c3(μa+μs).
The source flux s(x) in a volume element δV is given by
s(x)=δVU(x).
The Green’s function is the solution to Eq. 1, subject to the boundary condition imposed by the surface of the object. In general, numerical techniques such as Monte Carlo or finite-element modeling (FEM) are required to find a rigorous solution to Eq. 1 with a complex 3-D boundary. However, for the case of a homogeneous object, a tangential planar (TP) approximation can be made to reduce computational cost. In this approximation, shown schematically in Fig. 1 , the surface of the object is treated locally as an infinite plane boundary oriented tangentially to the surface at xi . The orientation of the plane changes for each surface element. This approximation is generally valid when the radius of curvature of the surface is greater than the effective absorption length,
1μeff=1[3μa(μa+μs)]12,
a condition that is easily satisfied in most practical cases. A similar technique has been proposed by Ripoll 23

Fig. 1

A schematic diagram illustrating the TP approximation where a plane boundary is drawn tangent to the surface at xi . The photon density is the solution for a point source at xj in a semi-infinite slab defined by the plane boundary, subject to the partial current boundary condition. Note that the orientation of the plane boundary is different for each surface element.

024007_1_031702jbo1.jpg

The Green’s function is the following form for a point source in the semi-infinite slab using the partial-current boundary condition24, 25, 26, 27 [see Eq. 2.4.2 in Ref. 27]:

Eq. 2

Gij=12πD(exp(μeffrij)rij1zb0dl exp(lzb)exp{μeff[(dj+l)2+pij2]12}[(dj+l)2+pij2]12).
Here, rij=[dj2+pij2]12 , where dj is the perpendicular distance from the source sj to the tangent plane, and pij is the distance along the plane from the surface element at xi to the perpendicular bisector of length dj , and
zb=2Dc1+Reff1Reff.
The parameter Reff is the average internal reflection coefficient,27 depends only on the index of refraction of the tissue underneath a surface element, and is typically in the range of 0.3 to 0.5. A simple analytic approximation (without integral) for Eq. 2 has been developed, in which the extrapolated boundary condition27 and Eq. 2 with pij=0 are combined by giving the former more weight for large rijpij , and giving the latter more weight for small rijpij . This simplified expression has excellent agreement ( <1% error in photon density predictions) with Eq. 2 for typical source depths and optical parameters, and can be calculated very quickly compared to evaluating the infinite integral.

2.2.

Converting Light Emission to a Photon Density Map

Photons emitted from the tissue surface are detected by a charge-coupled device (CCD)-based imaging system, as shown in Fig. 2 . The surface radiance L (photonsscm2steradian) is directly related to the photon density ρ (photonscm3) just inside the surface element. The exact form of the relationship depends on the model used to describe the transport of photons across the surface boundary, from tissue to air. Derived from the partial-current boundary condition [see Eq. 2.4.6 in Ref. 27], the relationship between L and ρ is

Eq. 3

L(θe)=c4πn2T(θ)[1+321Reff1+Reffcosθ]ρ,
where c is the speed of light, n is the index of refraction of the object medium, T is the transmission coefficient for light exiting the object through the surface element, and θ is the internal emission angle, which is related to the external emission angle θe through Snell’s law. The imaging system is absolutely calibrated such that electron counts on each CCD pixel can be mapped back onto the surface of the object, producing an absolute value of the surface radiance L from each imaged surface element, as illustrated schematically in Fig. 2. The imaging system collects the light emitted from the surface element at an angle θe (measured with respect to the normal to the surface element) into the solid angle dΩ subtended by the entrance pupil. We can apply Eq. 3 to convert the surface radiance measured at each surface element to values of the photon density directly beneath the surface.

Fig. 2

The object is divided into a solid mesh of voxels, one of which is shown here. Each voxel contains a point light source sj that contributes to the photon density ρi at the surface. The light emission from a surface element passes through the entrance pupil and is recorded in the image. The angle of emission with respect to the surface normal is θe . The entrance pupil subtends a small solid angle dΩ .

024007_1_031702jbo2.jpg

2.3.

Discretizing the Diffusive Transport Problem

Parameterization of the source intensity solution space is in the form of cubic voxels with a point light source at the center. We assign sj the value of the strength (or flux in photons/sec) of the point source inside the j ’th voxel. A 3-D grid of evenly spaced points defines the basis for source intensities estimated to approximate the actual source distribution. The reconstruction method takes advantage of the linear relationship between the source strength sj in each voxel and the photon density ρi measured at surface elements described by the Green’s function Gij [Eq. 2]. For m observations of photon density comprising column vector ρ , and n source points comprising column vector s , G represents the data kernel matrix of dimension m×n

Eq. 4

ρ=Gs.

2.3.1.

Wavelength dependence of Green’s functions

In the homogeneous tissue approximation, the inverse problem for source intensities comprised within the boundaries of a closed surface can be nonunique for surface radiance measured using a single wavelength band, particularly for a single-view perspective. In the case of a single point source embedded deep within a semi-infinite tissue slab, the resulting surface radiance could be fit equally well by a point source at the correct depth, or a large array of closely spaced sources at a shallower depth whose intensities are adjusted to match the surface radiance profile. Taking the surface shape into account, as in a mouse, helps to improve the resolution of the solution, but this is not a very strong constraint.

The resolution of the source parameters can be significantly improved using image data measured at different wavelengths. The Green’s function in Eq. 2 has a strong wavelength dependence due to the dependence on μeff(λ) . In the wavelength range of interest (λ550to700nm) , G decreases monotonically with wavelength as both μa(λ) and μs(λ) [and hence μeff(λ) ] generally decrease with wavelength in this window. If we measure the surface radiance at two or more wavelengths with significantly different G functions, and the wavelength dependence of G and s is known, then resolution of the solution to Eq. 4 is enhanced. In the semi-infinite slab example stated before, the incorrect shallow-source solution would not be allowed because this solution cannot fit the data for multiple wavelengths simultaneously.

Equation 4 can be rewritten to include multiple wavelengths as

5.

[ρ(λ1)ρ(λk)]=[η(λ1)G(λ1)η(λk)G(λk)][s]
ρ̃=G̃s,
where the wavelength dependence in elements of G(λk) comes into the terms D(λk) , μeff(λk) , and zb(λk) of Eq. 2. η(λk) is the relative fraction at which the wavelength λk contributes in the source emission spectrum, and is given by

Eq. 6

η(λk)=λkloλkupξ(λ)dλ0ξ(λ)dλ,
where λklo and λkup denote the lower and upper limits of the bandpass filter centered on wavelength λk , and ξ(λ) is the light source emission spectrum. η(λk) is introduced so that source flux unit (photons/second) is integrated over all wavelengths.

To illustrate improved resolution in a single-view system afforded by adding spectral information, singular value decomposition was performed for Green’s function kernel matrices for a source-detection geometry in a 2-D tissue slab oriented in the x-z plane. The 1×1-cm slab was pixelized into a 4×4 grid, with the pixel centers denoting source locations. Detection points were evenly spaced across the top boundary to emulate a single-view tomographic geometry. Green’s kernel matrices were computed for a number of wavelength combinations using in vivo mouse muscle optical properties (the experimental measurement is addressed in Sec. 4.1) and the firefly luciferase emission spectrum measured at 37°C . For each wavelength, 50 detection points were used. Singular value decomposition was performed on each of the kernel matrices representing the wavelength combinations, and the singular values were used to compute the condition number of the Green’s kernel matrices. Table 1 shows that for the discretized parameterization of source locations in tissue and detection points described before, the poor pose of the Green’s kernel matrix is improved, demonstrated by decreasing condition number with additional wavelength information.

Table 1

Condition numbers for Green’s kernel matrices.

N data50100150200
Wavelengths in Green’s kernel matrix 640nm 620, 640nm 600, 620, 640 nm580, 600, 620, 640 nm
Condition number 1.7×107 2.8×105 4.2×104 1.2×104

2.4.

Optimization for Source Estimation

The source intensities are determined by finding

Eq. 7

minW(ρ̃G̃s)22subjecttos0,
where W is a diagonal matrix of weights of dimension m , and 22 denotes Euclidean length squared. A constrained optimization algorithm28 is implemented to solve Eq. 5, to guarantee physical non-negative values for the estimated source intensities ŝ . The algorithm requires that mn , to ensure that the inverse problem is not underdetermined. The weight Wii is assigned to be the inverse of the instrument noise and statistical noise associated with the photon detection at surface element i wavelength k . Assuming Poisson distributed noise,
Wii1Ncounts+Nrd2,
where Ncounts is the number of photon counts per pixel and Nrd is the read noise of the camera system. The inverse weights, initially in units of counts, are converted to radiance units (photonssecondcm2steradian) , and then to photon density units (photonsmm3) as described in Sec. 2.2.

While high resolution localization is sought, a densely gridded source parametrization contributes to ill-conditioning of the inverse problem. To maintain high rank in the linear problem, the source parameterization is initially defined to be coarse, and an iterative locally adaptive gridding scheme is adopted. In each refinement step, solution voxels below a critical threshold (0.1% of the estimated source maximum) are suppressed in the source parameterization, and solution regions of high intensities are refined by decreasing the grid spacing by one half. The number of source parameters varies between refinements steps. Therefore, we use

χv2=1mnW(ρ̃G̃ŝ)22,
called the reduced χ2 , as the figure of merit to compare the goodness of fit between one refinement step to the next. Adaptively refined meshes have been utilized in fluorescent tomographic studies as well.29

After estimating source locations and intensities by the diffuse luminescence tomographic analysis, solutions can be characterized by calculating the total flux

stot=jŝj,
and the center of mass of the estimated source distribution
rcm=jrjŝjstot.
Source depths are estimated by taking the distance from the source center of mass to the surface, in the z direction.

3.

Experimental Setup

3.1.

Description of the Imaging System

The images in this study were obtained on an IVIS® 200 Imaging System (Xenogen Corporation), as shown in Fig. 3 . The system uses a back-thinned integrating CCD with high quantum efficiency over the spectral range of 400to950nm , which covers the tissue transmission window of 600to950nm . Cooling the camera to less than 90°C reduces the dark current to <100electronsscm2 and the associated noise to near-negligible levels.7 Imaging is performed with a high-sensitivity f/1 lens. The viewing platform for this instrument has mobility along the optical axis, where the fields of view can range from 3.9to26cm . The imaging system is absolutely calibrated to convert electron counts measured on the CCD, to surface radiance emitted at the imaging subject surface in units of photonsseccm2steradian . Conversion coefficients are calibrated for each f-stop and field of view. Multispectral capability in this instrument is possible through six 20-nm -wide bandpass filters spaced every 20nm from 560to660nm , which are set in the 24-position filter wheel. The instrument is equipped with a laser scanner that projects line patterns onto the animal, and this provides images from which the surface topography of the subject can be rendered.

Fig. 3

A schematic diagram of a cross sectional view of the IVIS 200 Imaging System. The laser galvanometer projects lines onto the animal for the structured light imaging feature.

024007_1_031702jbo3.jpg

3.2.

Structured Light Determination of Surface Topography

The surface topography is necessary to define the boundary condition of diffuse photon propagation at the tissue-air boundary. To determine the surface topography of the subject animal, a structured light technique has been implemented. Using a computer-controlled laser galvanometer, a series of parallel lines are projected down on the subject at a 20.6° angle. A structured light image of the illuminated subject is then acquired with the CCD, as shown in Fig. 4a . The animal's height can be determined by analyzing the bending or displacement of the lines as they pass over the animal. To determine surface topography from the structured light image, we utilize a two-step method that involves a 2-D Fourier transform30 combined with a phase-unwrapping algorithm.31

Fig. 4

Structured light analysis to determine surface topography: (a) single-view structured light image; (b) phase shift determined from fast Fourier transform (FFT), and (c) 3-D surface from phase unwrapping.

024007_1_031702jbo4.jpg

In Fourier space, a prominent peak near the reference frequency is broadened, due to the shift in phase of the lines when passing over the animal subject. The reference frequency is the line frequency at the center of the field of view on a flat surface. The constant and low-order terms, due to the positivity of image intensities and brightness of the animal relative to the black stage, also have significant amplitude in Fourier space. To compute the phase shifts relative to the reference frequency, the low-order terms and high-order harmonic terms must be filtered out using a bandpass filter. A cosine-tapered bandpass filter is applied to the Fourier spectrum to filter out the low-order terms, and high-order harmonics. The broadened peak selected out by the Fourier bandpass filter is then shifted to the Fourier space origin, and the inverse Fourier transform is applied to give sF(x,y)=b(x,y)exp[iφ(x,y)] , where the index F denotes the filtered image.

The phase

φ(x,y)=tan1{I[sF(x,y)]R[sF(x,y)]}
at each pixel location can then be calculated in the range φ(x,y)=(2π,2π] , as shown in Fig. 4b. These phase values must be unwrapped across the whole image to obtain the absolute phase. The unwrapping of phase across an image involves integration of phase values over a pathway through the image. As residues often occur as a result of image noise, a well-defined path avoiding residues is essential to correctly determining the surface heights. We have employed a quality-guided path phase-unwrapping technique, which first assigns a phase image subsection quality, and then directs the phase-unwrapping path based on the phase image subsection quality. The quality measure utilized here is the phase-derivative variance.31 Once the absolute phase is calculated for the image, the height can be determined approximately by the relation

Eq. 8

h(x,y)Lϕ(x,y)2πf0D,
where ϕ(x,y) is the absolute phase of the animal subject. The exact expression for Eq. 8 requires some correction terms due to the finite distance between the lens/projector and the animal subject.

The final unwrapped surface topography is shown in Fig. 4c. The line spacing and unwrapping algorithms give a surface accuracy of 0.5mm for smooth surfaces. Only the top surface of the animal facing the imaging system is rendered correctly with this technique. To construct a closed surface, the edge of the mouse surface is cropped vertically to the imaging platform and is closed horizontally across the bottom.

4.

Experimental Results

4.1.

Phantom Mouse Experiment

To assess the validity of the IVIS200 instrumentation and DLIT reconstruction algorithm, tests have been performed with a mouse-shaped homogeneous tissue phantom. To produce this phantom, a nude mouse was sacrificed, frozen, and the surface topography was scanned by a commercial machine vision scanner. The result was a high-resolution 3-D computer model from which a rubber mold was made. The rubber mold was then used to cast a phantom mouse with dimensions of approximately 3to4cm in width (side to side), 9cm in length (nose to feet), and 2.2cm in height (dorsal to ventral). The composition of the phantom is similar to that described by Vernon, 32 using a polyester resin as the host material, TiO2 for scattering, and Disperse Red for absorption. The optical properties were independently measured from a rectangular slab that was made from the same batch of material as the plastic mouse. The absorption μa and reduced scattering μs coefficients were measured using a double integrating sphere apparatus along with the inverse adding doubling method.33, 34 Curve fits to the intensity profile using a slab diffusion model with the partial current boundary condition were used as an alternate method to confirm the optical property values.7, 27 Two 200-μm fiber optics coupled to green LEDs were embedded into the phantom at two different known locations, labeled source A and source B. Depth measurements are from the fiber tip to the phantom surface along the z direction with respect to when the phantom is in a prone or supine position. Before insertion, the power emitted from each fiber optic was measured with a Newport Optical Power Meter (model 840), and converted to units of photons/sec for the wavelength of the spectral peak (567nm) . LED spectral profiles emitted from the fibers were characterized in an Ocean Optics USB2000 Fiber Optic Spectrometer. The profiles of the optical properties and spectra are shown in Fig. 5 . Each LED is independently controlled with a switch to allow for a combination of light-emitting source configurations. A schematic drawing of the phantom mouse is shown in Fig. 6 . An image sequence for tomographic analysis involves acquiring a photographic image (grayscale, illuminated reference image), a structured light image, and several luminescent images through different filters. These images were acquired for the mouse-shaped tissue phantom with the different source configurations and phantom positions. Luminescent data were acquired through four 20-nm -wide bandpass filters centered on 560, 580, 600, and 620nm . The phantom mouse was imaged in two positions: 1. with the dorsal side facing the camera, and 2. with the ventral side facing the camera. Three source configurations were imaged: source A, source B, and source A and B. Raw images of surface radiance for case 1, source A are shown in Fig. 7 .

Fig. 5

(a) Optical properties of the mouse-shaped tissue phantom. The dotted line is μs , the dashed line is μa , and the solid line is μeff . (b) The combined normalized wavelength spectrum of the two green LEDs coupled to fiber optics, which are embedded in the phantom.

024007_1_031702jbo5.jpg

Fig. 6

Schematic drawing of the mouse tissue phantom. The phantom mouse is made of polymer resin with added scatterer and absorbing dye to simulate optical properties of mouse muscle. Two LEDs are coupled to embedded fiber optics (lmax=570nm) , whose tips are labeled as source A and source B in the drawing. The distance of source A to the dorsal side of the phantom is lA=9.4±0.5mm , and the distance of source B to the ventral side of the phantom is lB=4.1±0.5mm .

024007_1_031702jbo6.jpg

Fig. 7

Surface radiance images of the phantom mouse dorsal side with embedded fiber optic source A measured through 20-nm -wide bandpass filters, centered at 560, 580, 600, and 620nm .

024007_1_031702jbo7.jpg

Structured light analysis was performed to render the portion of the phantom surface visible to the camera lens. Data acquired through each filter were selected by setting a dynamic range value for each filter image. The dynamic range is defined here to be the luminescent image maximum signal divided by an a priori threshold signal, below which is considered to be due to camera noise or animal background.35 The signal above the threshold given by the dynamic range is converted to photon density on the surface vertices, as described by Eq. 3. A fixed number of surface data elements for each filtered image are chosen, to equalize data weighting between the bandpass filtered data.

4.1.1.

Single-wavelength versus multiple-wavelength data

The resolution gained by incorporating wavelength-dependent Greens functions is demonstrated by testing combinations of filtered data used in the source estimation. In Fig. 8a , the source intensity distribution when source A was powered on was determined by data from a single wavelength filter centered on 620nm . 400 data points were included in the data vector. The broad source distribution reveals the low resolvability afforded by the model used for single wavelength data. In an inversion including data from the 580-nm centered wavelength band, where the optical properties are distinctly different from those at 620nm , the source localization becomes much tighter, as seen in Fig. 8b. The data vector consisted of 200 data points from each filtered dataset, for a total of 400 data points. The source intensity is also better resolved when using multiple spectra in the data, as shown in Table 2 .

Fig. 8

Diffuse luminescent tomographic estimates of LED source A for (a) single wavelength data at 620-nm filter, and (b) dual wavelength data at 580 and 620nm . The lateral views of each case are shown in (c) and (d), respectively. Source voxels are displayed on a red-black color scale in units of photons/sec. Voxels <2% of the max are suppressed.

024007_1_031702jbo8.jpg

Table 2

Tomographic estimates of the flux for Figs. 8a and 8b.

Source flux( ×1010 photons/sec)Relative error
Physical 7.7±0.4
620nm 4.148%
580nm+620nm 7.53%

4.1.2.

Four wavelength filtered data

Incorporating as much information as possible, image data of the phantom mouse from the four bandpass filters centered at 560, 580, 600, and 620nm were used in the inversion, and results are tabulated in Table 3 for different source and viewing configurations. For each filtered image, 200 data points were selected, giving a total of 800 data points. The source estimations were tightly localized, similar to the distribution in Figs. 8b and 8d. Accuracy in depth is within a millimeter, and intensities were estimated within 19% relative error for the LED sources. Voxel sizes were typically refined to side length of 1mm .

Table 3

DLIT estimates of fiber optic sources in the mouse-shaped tissue phantom. Data included luminescent images collected through the 560-, 580-, 600-, and 620-nm filters.

ExperimentMeasurementSource ASource B
Depth[mm]Flux(photons/sec)Depth[mm]Flux(photons/sec)
DorsalPhysical 9.4±0.5 7.7±0.4 ×1010 17.6±0.5 6.3±0.3 ×1010
Source A onTomographic9.8 7.8×1010
Source B onTomographic16.4 5.1×1010
VentralPhysical 9.9±0.5 7.7±0.4 ×1010 4.1±0.5 6.3±0.3 ×1010
Source A onTomographic10.3 7.1×1010
Source B onTomographic4.5 7.9×1010
Source A and BonTomographic10.3 7.5×1010 4.6 7.6×1010

The measured photon density at the visible surface is simulated well by the photon density forward modeled from the source estimate. The fit for 620-nm data is shown in Fig. 9 , where the photon density on the surface is displayed along the top row as a pseudocolor map overlaying the surface rendered in shaded white. Measured photon density data used in the fit lie within the pseudocolor map area shown on the left of the top row in Fig. 9, the simulated photon density overlaying the surface is shown in the middle, and the relative error map is shown on the right. Measured photon density values below the level of image noise were not used in the fit and were located at the surface area, which appears white in the left and middle maps. Line profiles drawn through the photon density maps are plotted and compared, displayed on the bottom row of Fig. 9.

Fig. 9

Top row: comparison of the measured (left) and simulated (middle), and percentage error (right) images of photon density at 620nm from the dorsal source A reconstruction of Table 3. Profiles on the bottom row show the photon density fit of the simulated photon density (dashed) to the measured photon density (solid) for a horizontal and vertical profile through mouse phantom, as indicated by the red cursor lines on the top row.

024007_1_031702jbo9.jpg

4.2.

In Vivo Studies

In vivo tissue is optically heterogeneous, and blood oxygenation levels cast additional fluctuations into the optical properties of the tissue. To assess the use of the homogeneous tissue approximation in the DLIT algorithm within heterogeneous tissue, we performed experiments with calibrated luminescent sources in vivo, where source locations and strengths are known. Tomographic results are also presented from metastasis model experiments in which tumor cells labeled with firefly luciferase are introduced into in vivo subjects.

4.2.1.

Calibrated luminescent sources

To benchmark the performance of the homogeneous tissue model adopted by the reconstruction algorithm, tritium-filled glass beads with interior phosphor coating were used as calibrated luminescent sources in vivo.5 Tritium is a β -emitter, and the absorption of the β by the phosphor induces photon emission with invariable power over the time scales of the experiments. The glass beads are 1mm in diameter and 3mm in length. Luminescent bead power was measured by placing the beads on black felt as an absorbing background and imaging them in the IVIS 200 Imaging System. Radiance measurements were integrated over the bead interior surface area and 2π steradians to convert to absolute flux. The bead spectrum was characterized on the Ocean Optics spectrometer.

In each of the following experiments, a 25-to30-g male nude mouse was anesthetized with a 4:1 concentration of ketamine and xylazine at a 15μl per 10-g body weight dose rate by interperitoneal injection. A 100-min period state of unconsciousness is expected. Once anesthesia was determined to have taken effect, luminescent beads were surgically implanted in specific locations. After implantation, surgical wounds were sutured and the animal was imaged in the camera system. The sutured incisions were located so that surface radiance from the beads could be imaged through undamaged skin. The animal was imaged in the camera using five bandpass filters centered at 560, 580, 600, 620, and 640nm . Before consciousness could be regained, the animals were euthanized using compressed carbon dioxide. Post-expiration, the bead locations were surgically determined and intact locations were gauged with calipers. In vivo optical properties for muscle and the normalized spectrum of the luminescent beads were used in the DLIT reconstruction. In vivo optical properties of a number of mouse organs were measured with the method described in Bevilacqua 36

Luminescent beads were implanted in contact with the left kidney of a male nude mouse, with one at the cranial end and another at the caudal end of the kidney. The DLIT solution using data from the 560-, 580-, and 640-nm bandpass filters produced the lowest reduced χ2 . The physical and tomographic measurements are outlined in Table 4 . Source depths and intensities are very well estimated, with less than 6% error in source intensity and 0.5-mm error in depth. Figure 10a shows the unfiltered raw data image. An anatomical mouse reference atlas showing brain, spinal column, and kidneys is displayed with the DLIT estimation inside the surface mesh in Fig. 10b. The sources estimated by DLIT are located at the cranial and caudal end of the reference atlas kidneys, suggesting tomographic agreement with the implantation sites.

Fig. 10

(a) Unfiltered image of luminescent beads implanted next to the left kidney. (b) Source reconstruction and organ atlas of mouse kidneys, brain, and spinal column in the surface mesh. Sources above 5% of the maximum intensity are displayed as voxels in red-black fade color.

024007_1_031702jbo10.jpg

Table 4

DLIT estimates of luminescent beads implantation near the left kidney, dorsal view.

MeasurementCranial beadCaudal bead
Depth(mm)Flux(ph/sec)Depth(mm)Flux(ph/sec)
Physical 3.4±0.5 3.5×1010 5.7±0.5 3.5×1010
DLIT3.2 3.7×1010 5.9 3.8×1010

In a second experiment, two luminescent beads were implanted in close proximity, between the right and left trapezius muscles of a male nude mouse. One bead was placed cranially, and the second was placed close to the first thoracic vertebra. Physical and tomographic measurements are summarized in Table 5 and Fig. 11 . For these shallow sources, the center of mass depths were estimated with submillimeter accuracy, and the sources are well localized. The distance between the estimated source center of masses was 7.1mm . Again, the source intensity estimation is in very good agreement with the measured power of the beads.

Fig. 11

(a) Unfiltered image of two luminescent beads implanted between the trapezius muscles. (b) Source reconstruction of luminescent beads implanted in-between the trapezius muscles. Sources above 2% of the maximum intensity are displayed as voxels in red-black fade color, and the mouse reference atlas organ brain and spinal column are displayed within the surface mesh.

024007_1_031702jbo11.jpg

Table 5

DLIT estimates of luminescent beads implanted between the trapezius muscles, dorsal view.

MeasurementCranial beadCervical bead
Depth(mm)Flux(ph/sec)Depth(mm)Flux(ph/sec)
Physical 1.5±0.5 3.5×1010 2.3±0.5 3.5×1010
DLIT1.3 4.0×1010 3.0 3.0×1010

4.2.2.

Tumor models

Monitoring metastatic disease in animals with DLIT was investigated by introducing 2×106 tumorigenic PC-3M-luc cells to a 30-gram male nude mouse by injection into the left ventricle of the heart. The cells were given time to colonize in vivo and develop into metastases. On day 29 postinjection of the PC-3M-luc cells, the animal was administered a dose of 15μl per 10-g body weight of luciferin by intraperitoneal injection. The luciferin kinetic profiles have been measured for this cell line and show that luciferin is taken up for approximately 10min postinjection, and a plateau is sustained for a 25-to30-min period postinjection. Anesthesia was induced with an isofluorane gas mixture. The animal was imaged dorsally in the camera 10min after injection of luciferin to guarantee sufficient signal. The ventral view was imaged 17min after the initial luciferin injection. Four bandpass filtered images were acquired at 580-, 600-, 620-, and 640-nm center frequencies, and the bandpass filtered image sequence was preceded and followed by an image acquired without a bandpass filter to evaluate the stability of luciferin kinetics during imaging. The imaging acquisition time of each filtered sequence totaled 3min .

The data were analyzed with DLIT using optical properties of in vivo muscle and the spectrum of firefly luciferase. The dorsal view image data acquired in the absence of a bandpass filter at the start of imaging were shown in Fig. 12a , and the DLIT source estimation results are shown in Fig. 12b. The same are shown for the ventral view in Fig. 13 . Source depth and intensity estimates are reported in Table 6 . From the dorsal view data, sources in the left femur and in the right peritoneal cavity are estimated by the tomographic algorithm. In the raw data, a bright spot radiating from the left forelimb skin is observed. On the ventral view, a source in the cardiac region is estimated by DLIT, in addition to the source in the right peritoneal cavity. Signal from the femur is not detected on the ventral view raw data image.

Table 6

DLIT estimates of PC-3M-luc lesions, dorsal and ventral views.

MeasurementLeft femurLeft cardiacPeritoneal cavity
Depth(mm)Flux(ph/sec)Depth(mm)Flux(ph/sec)Depth(mm)Flux(ph/sec)
DLIT dorsal1.9 1.0×107 6.2 3.6×107 10.1 3.9×108
DLIT ventral2.3 3.6×108 3.7 3.9×108

Fig. 12

(a) Dorsal view of unfiltered image of mouse injected with PC-3M-luc cells, 29days postinjection. (b) and (c) Top and perspective view of source reconstruction with atlas brain and spinal column are displayed. Voxels <1% of the maximum are suppressed.

024007_1_031702jbo12.jpg

From inspection of the analyzed dorsal data, we associate the source estimated in the animal's left forelimb with light emanating from the cardiac region, as seen from the ventral side of the animal in Fig. 13a. It is conceivable that light from the cardiac source propagated and leaked out toward the skin under the left forelimb. The limited perspective for this source on the dorsal side results in tomographic localization of the source next to the left forelimb and estimates the flux to be ten-fold lower compared to the ventral side estimate. In Fig. 13c, the mouse heart atlas organ overlay indicates that the reconstructed source is in the cardiac region. We interpret that some PC-3M-luc cells introduced to the mouse via intracardiac remained colonized there.

Fig. 13

(a) Ventral view of unfiltered image of mouse injected with PC-3M-luc cells, 29days postinjection. (b) Source reconstruction with atlas brain and spinal column are displayed. Voxels <5% of the maximum are suppressed. (c) is the same as (b), except the heart atlas organ is displayed as well.

024007_1_031702jbo13.jpg

Despite the optically heterogeneous intestinal organs in which the peritoneal source may be embedded, the estimated source intensity is in agreement between the dorsal and ventral views. The location of the peritoneal source is estimated to be 3.8mm deep from the ventral surface and 10.1mm from the dorsal surface. The depths from the opposing surfaces add to 14mm . Measurements of the thickness of this mouse from lateral view images gives 14to15mm at the approximate location of the source, and suggests colocalization by DLIT from the dorsal and ventral sides.

5.

Discussion

A new technique called diffuse luminescence imaging tomography (DLIT) has been developed for obtaining a 3-D reconstruction of the bioluminescent source distribution inside a complex object such as an animal subject. The technique is based on the analysis of multiple images of the light emission from a single perspective of the object acquired through a number of bandpass filters. The image data are combined with an independent measurement of the surface topography to produce a high-resolution map of the photon density at the surface. The advantage of using imaging to characterize the surface light emission is that it allows the complete characterization of subjects having complex surfaces, such as a mouse. This technique can yield in real time detailed information about the strength and position of the internal light sources.

The reconstruction algorithm consists of finding an approximate solution to a system of linear equations that relate the source strength at each point inside the object to the photon density at the surface. We have presented a simple form for the Green’s function, called the tangential plane (TP) approximation, based on the diffusion model for photon transport. The TP approximation takes into account the effective reflectivity of the surface and can be expressed in a very computationally efficient manner. Light sources are estimated by minimizing forward modeled photon density misfit to the data using a non-negative least-squares algorithm.

As a demonstration of the DLIT technique, we have determined the 3-D source distribution from measured single-perspective, multiple-bandpass spectra images of a phantom mouse and live mice with one or two calibrated internal light sources. In these test cases, the reconstruction technique was able to locate the multiple sources with a depth accuracy of <1mm , and 1 to 20% relative error in intensity, for sources inside the complex object surface. An in vivo mouse tumor model experiment was conducted in which the mouse was imaged viewing its dorsal and ventral side. A tumor producing prominent signal on both sides was identified at depth in the peritoneal cavity by DLIT, and the light intensity estimates from the dorsal and ventral side were in agreement. Source strength and location were retrieved with high accuracy in the in vivo bead studies for a number of implantation sites, and we find the agreement in source strength of the deep peritoneal tumor source to be encouraging. The reconstruction algorithm is very fast, requiring about a minute on a 3.4-GHz PC for iterations on 800 data elements and adaptively refined solution space of up to 200 elements.

Our results to date suggest that the homogeneous approximation appropriately determines light source location and intensity for sources in a number of regions in the living mouse. However, for photon propagating through more pronounced heterogeneous paths, through the liver and lungs, for example, the tomographic model may require a more sophisticated description of photon transport and in vivo tissue optical property measurements for more accurate source estimations in living animal models. Methods to extend the DLIT algorithm to include heterogeneous tissue properties are facilitated by the digital mouse atlas and are under investigation.

The tomographic results presented here are in the context of a single-view imaging system, which has the advantage of high throughput, but the disadvantage in that not all sides of the mouse can be imaged simultaneously. The DLIT algorithm can be applied equally well to multiview optical imaging systems, and this will be discussed in a future publication.

Acknowledgments

The authors would like to acknowledge John Hunter for designing the tumor model experiments, and Haroon Ahsan and Joan Dusich for their technical support with the animal experiments. We would like to extend our gratitude to Dan Stearns for his valuable contributions to the foundation of this work.

References

1. 

C. H. Contag, D. Jenkins, P. R. Contag, and R. S. Negrin, “Use of reporter genes for optical measurements of neoplastic disease in vivo,” Neoplasia, 2 41 –52 (2000). https://doi.org/10.1038/sj.neo.7900079 1522-8002 Google Scholar

2. 

A. Rehemtulla, L. D. Stegman, S. J. Cardozo, S. Gupta, D. E. Hall, C. H. Contag, and B. D. Ross, “Rapid and quantitative assessment of cancer treatment response using in vivo bioluminescence imaging,” Neoplasia, 2 491 –495 (2000). https://doi.org/10.1038/sj/neo/7900121 1522-8002 Google Scholar

3. 

K. P. Francis, J. Yu, C. Bellinger-Kawahara, D. Joh, M. J. Hawkinson, G. Xiao, T. F. Purchio, M. G. Caparon, M. Lipsitch, and P. R. Contag, “Visualizing pneumococcal infections in the lungs of live mice using bioluminescent streptococcus pneumoniae transformed with a novel gram-positive lux transposon,” Infect. Immun., 69 3350 –3358 (2001). https://doi.org/10.1128/IAI.69.5.3350-3358.2001 0019-9567 Google Scholar

4. 

S. Gross and D. Piwnica-Worms, “Real-time imaging of ligand-induced IKK activation in intact cells and in living mice,” Nat. Methods, 2 607 –614 (2005). 1548-7091 Google Scholar

5. 

M. Fowler, J. Virostko, Z. Chen, G. Poffenberger, A. Radhika, M. Brissova, M. Shiota, W. E. Nicholson, Y. Shi, B. Hirshberg, D. M. Harlan, E. D. Jansen, and A. C. Powers, “Assessment of pancreatic islet mass after islet transplantation using in vivo bioluminescence imaging,” Transplantation, 79 768 –776 (2005). https://doi.org/10.1097/01.TP.0000152798.03204.5C 0041-1337 Google Scholar

6. 

V. Tuchin, Tissue Optics, SPIE Press, Bellingham, WA (2000). Google Scholar

7. 

B. W. Rice, M. D. Cable, and M. B. Nelson, “In vivo imaging of light emitting probes,” J. Biomed. Opt., 6 (4), 432 –440 (2001). https://doi.org/10.1117/1.1413210 1083-3668 Google Scholar

8. 

H. Zhao, T. C. Doyle, O. Coquoz, F. Kalish, B. W. Rice, and C. H. Contag, “Emission spectra of bioluminescent reporters and interaction with mammalian tissue determine sensitivity of detection in vivo,” J. Biomed. Opt., 10 (4), 41210 (2005). 1083-3668 Google Scholar

9. 

G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys., 31 2289 –2299 (2004). https://doi.org/10.1118/1.1766420 0094-2405 Google Scholar

10. 

X. J. Gu, Q. Z. Zhang, L. Larcom, and H. B. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express, 12 3996 –4000 (2004). https://doi.org/10.1364/OPEX.12.003996 1094-4087 Google Scholar

11. 

W. X. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. H. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “A practical reconstruction method for bioluminescence tomography,” Opt. Express, 13 6756 –6771 (2005). https://doi.org/10.1364/OPEX.13.006756 1094-4087 Google Scholar

12. 

W. X. Cong and G. Wang, “Boundary integral method for bioluminescence tomography,” J. Biomed. Opt., 11 (2), 020503 (2006). https://doi.org/10.1117/1.2191790 1083-3668 Google Scholar

13. 

W. X. Cong, K. Durairaj, L. V. Wang, and G. Wang, “A Born-type approximation method for bioluminescence tomography,” Med. Phys., 33 679 –686 (2006). https://doi.org/10.1118/1.2168293 0094-2405 Google Scholar

14. 

V. Ntziachristos and R. Weissleder, “Charge-coupled-device based scanner for tomography of fluorescent near-infrared probes in turbid media,” Med. Phys., 29 803 –809 (2002). https://doi.org/10.1118/1.1470209 0094-2405 Google Scholar

15. 

V. Ntziachristos, C. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med., 8 757 –760 (2002). 1078-8956 Google Scholar

16. 

V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized Born approximation,” Opt. Lett., 26 893 –895 (2001). 0146-9592 Google Scholar

17. 

C. Kuo, H. Ahsan, J. Hunter, T. Troy, H. Xu, N. Zhang, and B. Rice, “In vivo bioluminescent tomography using multi-spectral and multi-perspective image data,” Biomedical Optics 2006 Technical Digest, (2006) Google Scholar

18. 

G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol., 50 4225 –4241 (2005). https://doi.org/10.1088/0031-9155/50/17/021 0031-9155 Google Scholar

19. 

O. Coquoz, T. L. Troy, D. Jekic-McMullen, and B. W. Rice, “Determination of depth of in vivo bioluminescent signals using spectral imaging technique,” Proc. SPIE, 4967 37 –45 (2003). https://doi.org/10.1117/12.477885 0277-786X Google Scholar

20. 

H. Dehghani, S. C. Davis, S. D. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett., 31 365 –367 (2006). https://doi.org/10.1364/OL.31.000365 0146-9592 Google Scholar

21. 

A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol., 50 5421 –5441 (2005). https://doi.org/10.1088/0031-9155/50/23/001 0031-9155 Google Scholar

22. 

A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York (1978). Google Scholar

23. 

J. Ripoll, V. Ntziachristos, R. Carminati, and M. Nieto-Vesperinas, “Kirchoff approximation for diffusive waves,” Phys. Rev. E, 64 051917 (2001). https://doi.org/10.1103/PhysRevE.64.051917 1063-651X Google Scholar

24. 

M. Keijzer, W. M. Star, and P. R. M. Storchi, “Optical diffusion in layered media,” Appl. Opt., 27 1820 –1824 (1988). 0003-6935 Google Scholar

25. 

A. Lagendijk, R. Vreeker, and P. DeVries, “Influence of internal reflection on diffusive transport in strongly scattering media,” Phys. Lett. A, 136 81 –88 (1989). https://doi.org/10.1016/0375-9601(89)90683-X 0375-9601 Google Scholar

26. 

J. X. Zhu, D. J. Pine, and D. A. Weitz, “Internal reflection of diffusive light in random media,” Phys. Rev. A, 44 3948 –3959 (1991). https://doi.org/10.1103/PhysRevA.44.3948 1050-2947 Google Scholar

27. 

R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion Eq. in radiative transfer,” J. Opt. Soc. Am. A, 11 2727 –2741 (1994). 0740-3232 Google Scholar

28. 

C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, 160 –165 Prentice-Hall, Englewood Cliffs, NJ (1974). Google Scholar

29. 

A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express, 12 5402 –5417 (2004). https://doi.org/10.1364/OPEX.12.005402 1094-4087 Google Scholar

30. 

M. Takeda, H. Ina, and S. Kobayshi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Acoust. Soc. Am., 72 156 –160 (1982). 0001-4966 Google Scholar

31. 

D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping, Theory, Algorithms, and Software, John Wiley and Sons, New York (1998). Google Scholar

32. 

M. L. Vernon, J. Frechette, Y. Painchaud, S. Caron, and P. Beaudry, “Fabrication and characterization of a solid polyurethane phantom for optical imaging through scattering media,” Appl. Opt., 38 4247 –4251 (1999). 0003-6935 Google Scholar

33. 

S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding doubling method,” Appl. Opt., 32 559 –568 (1993). 0003-6935 Google Scholar

34. 

J. W. Pickering, S. A. Prahl, N. van Wieringen, J. F. Beek, J. Sterenborg, and M. J. van Gemert, “Double integrating sphere system for measuring the optical properties of tissue,” Appl. Opt., 32 339 –410 (1993). 0003-6935 Google Scholar

35. 

T. Troy, D. Jekic-McMullen, L. Sambucetti, and B. Rice, “Quantitative comparison of the sensitivity of detection of fluorescent and bioluminescent reporters in animal model,” Mol. Imaging, 3 9 –23 (2004). https://doi.org/10.1162/153535004773861688 1535-3508 Google Scholar

36. 

F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt., 38 4939 –4950 (1999). 0003-6935 Google Scholar
©(2007) Society of Photo-Optical Instrumentation Engineers (SPIE)
Chaincy Kuo, Olivier Coquoz, Tamara L. Troy, Heng Xu, and Bradley W. Rice "Three-dimensional reconstruction of in vivo bioluminescent sources based on multispectral imaging," Journal of Biomedical Optics 12(2), 024007 (1 March 2007). https://doi.org/10.1117/1.2717898
Published: 1 March 2007
Lens.org Logo
CITATIONS
Cited by 173 scholarly publications and 2 patents.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Natural surfaces

Tomography

In vivo imaging

Tissue optics

Bandpass filters

Tissues

3D image processing

Back to Top