|
1.IntroductionBioluminescent 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 and lower in the red and near-infrared (NIR) part of the spectrum .6 Firefly luciferase has a rather broad emission spectrum ranging from , 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 , 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 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.Theory2.1.Photon Diffusion in TissueLight 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 be much greater than the absorption coefficient , so that the change in the photon density is very small between scattering events. In this case, the photon density in a homogeneous medium is governed by the steady-state diffusion equation: where is the photon rate density and is the diffusion coefficient,The source flux in a volume element is given byThe 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 . 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,a condition that is easily satisfied in most practical cases. A similar technique has been proposed by Ripoll 23The 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]: Here, , where is the perpendicular distance from the source to the tangent plane, and is the distance along the plane from the surface element at to the perpendicular bisector of length , andThe parameter 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 are combined by giving the former more weight for large , and giving the latter more weight for small . This simplified expression has excellent agreement ( 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 MapPhotons emitted from the tissue surface are detected by a charge-coupled device (CCD)-based imaging system, as shown in Fig. 2 . The surface radiance is directly related to the photon density 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 and is where is the speed of light, is the index of refraction of the object medium, 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 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 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 (measured with respect to the normal to the surface element) into the solid angle 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.2.3.Discretizing the Diffusive Transport ProblemParameterization of the source intensity solution space is in the form of cubic voxels with a point light source at the center. We assign the value of the strength (or flux in photons/sec) of the point source inside the ’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 in each voxel and the photon density measured at surface elements described by the Green’s function [Eq. 2]. For observations of photon density comprising column vector , and source points comprising column vector , represents the data kernel matrix of dimension 2.3.1.Wavelength dependence of Green’s functionsIn 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 . In the wavelength range of interest , decreases monotonically with wavelength as both and [and hence ] generally decrease with wavelength in this window. If we measure the surface radiance at two or more wavelengths with significantly different functions, and the wavelength dependence of and 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 where the wavelength dependence in elements of comes into the terms , , and of Eq. 2. is the relative fraction at which the wavelength contributes in the source emission spectrum, and is given bywhere and denote the lower and upper limits of the bandpass filter centered on wavelength , and is the light source emission spectrum. 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 plane. The slab was pixelized into a 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 . 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 1Condition numbers for Green’s kernel matrices.
2.4.Optimization for Source EstimationThe source intensities are determined by finding where is a diagonal matrix of weights of dimension , and 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 , to ensure that the inverse problem is not underdetermined. The weight is assigned to be the inverse of the instrument noise and statistical noise associated with the photon detection at surface element wavelength . Assuming Poisson distributed noise,where is the number of photon counts per pixel and is the read noise of the camera system. The inverse weights, initially in units of counts, are converted to radiance units , and then to photon density units 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 called the reduced , 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.29After estimating source locations and intensities by the diffuse luminescence tomographic analysis, solutions can be characterized by calculating the total flux and the center of mass of the estimated source distributionSource depths are estimated by taking the distance from the source center of mass to the surface, in the direction.3.Experimental Setup3.1.Description of the Imaging SystemThe 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 , which covers the tissue transmission window of . Cooling the camera to less than reduces the dark current to 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 . 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 . Conversion coefficients are calibrated for each f-stop and field of view. Multispectral capability in this instrument is possible through six -wide bandpass filters spaced every from , 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. 3.2.Structured Light Determination of Surface TopographyThe 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 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 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 , where the index denotes the filtered image. The phase at each pixel location can then be calculated in the range , 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 relationwhere 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 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 Results4.1.Phantom Mouse ExperimentTo 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 in width (side to side), in length (nose to feet), and 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, 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 and reduced scattering 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 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 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 . 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 -wide bandpass filters centered on 560, 580, 600, and . 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 . 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 dataThe 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 . 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 centered wavelength band, where the optical properties are distinctly different from those at , 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 . Table 2Tomographic estimates of the flux for Figs. 8a and 8b.
4.1.2.Four wavelength filtered dataIncorporating as much information as possible, image data of the phantom mouse from the four bandpass filters centered at 560, 580, 600, and 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 relative error for the LED sources. Voxel sizes were typically refined to side length of . Table 3DLIT 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.
The measured photon density at the visible surface is simulated well by the photon density forward modeled from the source estimate. The fit for 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. 4.2.In Vivo StudiesIn 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 sourcesTo 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 in diameter and 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 steradians to convert to absolute flux. The bead spectrum was characterized on the Ocean Optics spectrometer. In each of the following experiments, a male nude mouse was anesthetized with a 4:1 concentration of ketamine and xylazine at a per body weight dose rate by interperitoneal injection. A 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 . 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 bandpass filters produced the lowest reduced . 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 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. Table 4DLIT estimates of luminescent beads implantation near the left kidney, dorsal view.
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 . Again, the source intensity estimation is in very good agreement with the measured power of the beads. Table 5DLIT estimates of luminescent beads implanted between the trapezius muscles, dorsal view.
4.2.2.Tumor modelsMonitoring metastatic disease in animals with DLIT was investigated by introducing tumorigenic PC-3M-luc cells to a 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 per 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 postinjection, and a plateau is sustained for a period postinjection. Anesthesia was induced with an isofluorane gas mixture. The animal was imaged dorsally in the camera after injection of luciferin to guarantee sufficient signal. The ventral view was imaged after the initial luciferin injection. Four bandpass filtered images were acquired at 580-, 600-, 620-, and 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 . 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 6DLIT estimates of PC-3M-luc lesions, dorsal and ventral views.
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. 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 deep from the ventral surface and from the dorsal surface. The depths from the opposing surfaces add to . Measurements of the thickness of this mouse from lateral view images gives at the approximate location of the source, and suggests colocalization by DLIT from the dorsal and ventral sides. 5.DiscussionA 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 , 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 PC for iterations on 800 data elements and adaptively refined solution space of up to 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. AcknowledgmentsThe 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. ReferencesC. 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
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
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
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
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
V. Tuchin, Tissue Optics, SPIE Press, Bellingham, WA
(2000). Google Scholar
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York
(1978). Google Scholar
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
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
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
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
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
C. L. Lawson and
R. J. Hanson, Solving Least Squares Problems, 160
–165 Prentice-Hall, Englewood Cliffs, NJ
(1974). Google Scholar
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
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
D. C. Ghiglia and
M. D. Pritt, Two-Dimensional Phase Unwrapping, Theory, Algorithms, and Software, John Wiley and Sons, New York
(1998). Google Scholar
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
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
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
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
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
|