A fast 3-D optical imaging method with use of exogenous fluorescence agent is proposed and demonstrated by simulation in a model tissue. After administration of fluorescent agent, ultrashort near-infrared laser pulses are used to illuminate the tissue and excite fluorescence emission. The transient fluorescence signals are detected on the tissue boundaries and employed to reconstruct a 3-D image of relative fluorescence emission distribution inside the tissue. A region with greater fluorescence emission represents a diseased tissue if the fluorescent agent has a close affinity with the disease. We successfully demonstrated the feasibility of this method in the imaging of a small cubic tumor embedded in a cubical tissue phantom with a preassigned uptake distribution of fluorescent indocyanine green dye. The image reconstruction does not involve any inverse optimization. It took less than 5 minutes in a general PC for the two model imaging problems.
©2004 Optical Society of America
In recent years, studies on near-infrared (NIR) optical imaging as a potential biomedical diagnostic modality have attracted increasing attention [1–7]. Unlike conventional optical imaging that uses endogenous contrast of tissue optical properties between normal and diseased tissues to detect and/or diagnose diseases, optical fluorescence imaging utilizes exogenous fluorescent agents to enhance and/or characterize optical contrast between normal and diseased tissues. The use of fluorescence techniques has appealing features like selective absorption spectrum of contrast agents, distinct fluorescence spectrum and lifetime, cellular and molecular fluorescence probes, etc. These traits can be used to provide tissue functional information, cellular and molecular information, as well as structural information.
Two common diffuse optical imaging methods are the time-domain and frequency-domain photon migration measurements. The most challenging task in diffuse optical imaging is the fact that NIR light is highly scattered in human tissues. Unlike X-ray and PET that depend on the ballistic transmitted component of radiation emission, optical imaging has to deal with highly diffused and back-scattered light, and usually involves an inverse problem in image reconstruction, in which tissue optical properties from a given experimental setup and a given set of measurements are determined through an inverse optimization procedure using predicted boundary measurements in forward modeling of photon migration in the given tissue geometry with known optical tissue parameters.
The majority researchers used iterative methods to solve an inverse problem through a convergent algorithm of appropriate objective functions in order to find distributions of tissue optical parameters. For example, O’Leary et al.  formulated an inverse problem for fluorescence lifetime tomography using diffuse photon density waves. Ntziachristos et al.  presented quantitative optical imaging of human breast in vivo by using NIR diffuse optical tomography (DOT) after the administration of indocyanine green (ICG) for contrast enhancement and examined the optical imaging concurrently with magnetic resonance imaging (MRI). Recently Godavarty et al.  developed three-dimensional (3-D) fluorescence tomographic imaging in the frequency-domain using NIR contrast agents in a large tissue-mimic phantom.
Diffusion model based on a diffuse approximation equation was adopted for predicting photon migration in tissues (the so-called forward modeling) in the majority of optical imaging studies. For instance, Wu et al.  investigated fluorescence tomographic imaging using early-arriving photons and Laplace transforms of an analytical diffusion approximation of photon migration. However, the diffusion model is only an approximation to the rigorous radiative transfer equation (RTE) and its applicability is limited to scattering-dominant thick media. Biologic tissues are heterogeneous, and contain some regions that are highly absorbing and/or lowly scattering. For more accurate forward modeling, Chang et al.  utilized Fourier transform to convert time-resolved RTEs to the frequency domain for describing photon migration in highly scattering media. Klose and Hielscher  also considered fluorescence image reconstruction using the stationary RTE as the forward model. Guo and his co-authors solved the transient RTE in time-domain for ultrafast laser transport in biologic tissues using the discrete-ordinates method (DOM) [8, 9] and the radiation element method . An RTE-based forward modeling may consume several hours in computation using a general PC. Even with the use of a simplified diffusion model, a forward modeling may still take dozens of minutes for 3-D heterogeneous tissue geometries.
In the same time, an inverse problem in image reconstruction needs dozens of to hundreds of times of forward modeling. Given the complexity of both the distributions of tissue optical properties and a forward modeling approach, diffuse optical imaging may require several hours to several days to reconstruct an image. Further, an inverse problem is always ill-posed and may lead to divergent results.
In this treatise, a fast and stable optical imaging method is proposed to detect tumor based on relative fluorescence concentration contrast between cancerous tumor and surrounding healthy tissue, rather than the inverse solution of the distribution of tissue optical parameters. Variation in fluorescence concentration distribution due to preferential uptake or washout in diseased tissues leads to detection of diseases. Here we first simulate temporal fluorescence signals using the discrete-ordinates method for a model fluorescence imaging system. The simulated fluorescence data in detectors located on the tissue boundaries are then used for image reconstruction in which a relative fluorescence emission source distribution is established by allocating the detected fluorescence signals to tissue internal voxels. The area with greater fluorescence emission is considered to be cancerous region because the fluorescent agent has a close affinity with cancer. Finally, a reconstructed 3-D image and tomographic images will be presented and compared with the actual tumor shape, size, and location. In actual experimental study in the future, the simulation of fluorescence signals is not needed. Since we only reconstruct the fluorescence emission distribution, the many difficulties associated with the complexities of tissue optical properties, tissue structures, forward modeling, and inverse optimization, are all avoided in the proposed imaging method. Thus, the image reconstruction is expected to be very fast and stable. We will demonstrate our method in the imaging of a small tumor embedded in a model tissue.
2. Simulation and imaging analysis
2.1 Simulation of fluorescence signals
The tissue model for 3-D fluorescence imaging is shown in Fig. 1 where a cubical tissue phantom of side length of 20mm is embedded with a tumor–like object of 4×4×4mm3 in size at the tissue center. After administration of fluorescent agents, exciting laser pulses are incident from the center of one surface of the tissue phantom (e.g., from the plane at X=0) and the consequent fluorescence signals are received by detector arrays located on the four side surfaces of the tissue.
The fluorescence signals are simulated using the discrete-ordinates method for the coupled radiation transfer equations that govern photon transport of both the excitation laser and consequent fluorescence in the model tissue:
where c is the speed of light in the tissue, Il the light intensity in a discrete direction sl with three direction cosines ξ, η and µ, and σ e the extinction coefficient that is the sum of absorption coefficient σa and scattering coefficient σs. The subscripts 1 and 2 in the symbols denote the excitation laser and fluorescence light, respectively. The radiation source terms due to scattering can be formulated in the discrete-ordinates format as
in which Φ il represents the scattering phase function, and ω the scattering albedo defined as ω=σ s /σ e . A quadrature set of n discrete ordinates with appropriate angular weights wi (i=1, 2, …, n) is used. The choice of quadrature scheme is arbitrary, and n=N(N+2) in a S N DOM scheme.
is the source contribution of the collimated laser irradiation and can be written as
where ξc represents the collimated laser incident direction, δ the Dirac delta function, and I 0 the laser beam intensity irradiated to the surface of the tissue.
is the excited fluorescence emission  and can be formulated as
where Cv, αF , and τ are the molarity, quantum yield, and lifetime of fluorescence agent, respectively, and ν is the frequency of radiation.
2.2 Image reconstruction
A typical time-resolved fluorescence signal from the present simulation results is shown in the solid curve in Fig. 2. This signal can be divided into two parts. The rising part to the left of the maximum fluorescence intensity I max depends on early arriving photons of ballistic and snake components, and provides timing information, while the falling part after the I max is affected by fluorescence lifetime. Of course, the whole curve is also strongly influenced by the medium property through which the fluorescent light passed. Let T’i,j,k be the time measured at half I max position in the rising part and T’mi,j,k be the time measured at the peak intensity point. The indices of ‘i’, ‘j’ and ‘k’ indicate the sequence numbers of a detector in the X, Y and Z directions, respectively. Both T’ i,j,k and T’mi,j,k include photon flight time of the excitation laser from incident position to tumor position and the consequent fluorescent photon flight time from the tumor (the fluorescent agent has a close affinity with cancer) to detector (i,j,k). Therefore, the fluorescence flight times from the tumor to detector (i,j,k) at half and full I max, respectively, are corrected as follows:
in which the subscript 0 represents the detector located at the laser incident position. The times T’0 and T’m0 equal to twice of the photon flight time from the laser incident position to the tumor position, respectively measured at half and full I max.
The majority of fluorescent photons cannot travel straightly to the tissue boundary because of multiple scattering events associated with highly scattering turbid media [5–10]. It is nearly impossible to directly measure the photon flight time along a straight line between the tumor center and a detector except for the case that the tissue is extremely thin which might have no real meaning for noninvasive tumor detection. Nevertheless, we need to estimate the distance between the tumor center and a detector. We hypothesize that the flight time Ti,j,k of early arriving fluorescence photons (the rising part in the fluorescence signal curve) is governed by a Gaussian probability distribution:
where P is the probability, and σ2 is the variance of the distribution that can be expressed by
In order to see whether the rising curve in Fig. 2 can be assumed as a Gaussian profile, we have drawn a dotted curve that is symmetric to the rising curve against the I max line. Clearly these two curves form a distribution that is very close to a Gaussian profile. We have examined this assumption using the numerical data of this study and previous experimental data of laser transport in highly scattering media . The outcome was satisfactory.
Now the fluorescence flight time along a straight line between the tumor center and a detector is estimated as
The average distance between the tumor and a detector located at (i, j, k) is then calculated by
The equation for a spherical surface centered at detector (i, j, k) with a radius of Ri,j,k is
in which (Di , Dj , Dk ) are the coordinates of the detector location. We define a voxel space line as X=Vi ’ and Y=Vj ’, and solve for the intersection points between the spherical surface and the space line.
When -(Vi′ ’-Di )2-(Vj′ ’-Dj )2≥0, the intersection points are
Otherwise, there is no intersection point. (Vi’,’,Vj’ ’,Vk ’) is the voxel and/or pixel grid position for image reconstruction.
If the intersections are located in the tissue geometry, the value of I max detected by this detector is then allocated to the voxel (Vi’ ’,Vj’ ’,Vk ’) and its adjacent eight voxels. Repeating the same procedure for every detector, we can get a relative intensity of fluorescence emission from each voxel contributed from all detectors:
Here, E(Vi ’,Vj ’,Vk ’) is the relative fluorescence intensity of voxel (Vi ’,Vj ’,Vk ’), ID (Di , Dj , Dk ; Vi ’,Vj ’,Vk ’) is the intensity contributed from detector (Di , Dj , Dk ) to voxel (Vi ’,Vj ’,Vk ’), and I max(Di , Dj , Dk ) is the maximum fluorescence intensity at detector (i, j, k). The computational grid is N×M×K.
In the simulation model, we considered pulsed excitation light of pulse width 1 ps and wavelength at 780 nm. ICG dye was used as the fluorescent agent. Fluorescence signals were received for wavelength greater than 800 nm. The absorption coefficient and reduced scattering coefficient of the model tissue were 0.0023mm-1 and 1.018mm-1 against the excitation laser, and 0.0031mm-1 and 0.864mm-1 against the fluorescence, respectively. The absorption coefficient of the ICG solution was 0.03mm-1 against the excitation laser (assuming Cv0=1µM), and 0.0025mm-1 against the fluorescence. Those property values were used by Godavarty et al. . We considered an uptake distribution of the fluorescent dye in the model tissue. The dimensionless molarity Cv/Cv0 of the ICG solution was between 0.9 and 1.0 in the tumor region, while it was randomly distributed in the surrounding healthy tissue with values between 0.01 and 0.06. The uptake ratio varied between 15:1 and 100:1. The coupled radiation transfer equations were solved with the S 10 DOM method . A computational grid of 25×25×25 was used in the simulation.
Figure 3 exhibits the reconstructed 3-D tumor image for the imaging model shown in Fig. 1. A filter technique based on the gray-level histogram concept  was used to enhance the image contrast and quality. The blue block in Fig. 3 has greater fluorescence emission value and stands for the embedded cubical tumor. In order to see clearly the location of the reconstructed tumor image, the 3-D image is projected onto three planes normal to the X, Y, and Z directions, respectively. The three projections are shown in Figs. 4a, 4b, and 4c. The squares enclosed with black dash-dot lines show the actual tumor location, size, and shape. The red areas enclosed with bright green lines are the location, size, and shape of the detected tumor. It is seen that the images closely match the actual embedded tumor. Slight displacements of tumor image locations are observed. This is caused by the estimation of the fluorescence flight time along a straight path in order to get the distance information. It is very difficult to precisely measure this timing information because of the highly scattering phenomenon of NIR light in biologic tissues.
Tomographic images at any plane of any orientation are easily obtainable by slicing the 3-D image. Figures 5a, 5b, and 5c display the tomographic images along the X-, Y-, and Z-axis, respectively. The eight slices in each figure are selected in the respective axis locations from 7.5mm to 12.5mm, where the detected tumor shape, size, and location at different axis positions in the model tissue are illustrated clearly. The tomographic images sliced around the tumor center capture the actual tumor very well. However, deformations are found for tomographic images sliced at the edge of the tumor. Because the cubical tumor has an abrupt or discontinuous change at its boundaries, it is a formidable task to reconstruct a sharp image for any imaging method. Also such deformations are partially attributable to the varied fluorescence solution uptake. We have assumed a steady change of fluorescent agent uptake from 0.90 at the edge of the tumor to 1.0 at the tumor center and a random distribution of agent uptake in the background. Considering the facts that the tumor image was retrieved from a certain depth inside a highly scattering dense tissue and the embedded model tumor has a very sharp shape, however, the deformations and slight displacements are reasonably acceptable. As a matter of fact, the present images have excellent quality and high contrast as compared with most published optical imaging results in the literature.
In order to validate our method, we further imaged a small tumor embedded at an off-center position of the tissue. Figure 6 displays the reconstructed 3-D tumor image. The tumor was centered at the position of X=15 mm, Y=15 mm, and Z=15 mm. The size of the tumor was reduced to 2.4×2.4×2.4mm3. Other conditions remain the same like the preceding model. The detected tumor in Fig. 6 closely matches the preassigned tumor location and size conditions.
A unique feature of the present method is that it is extremely computationally efficient in the image reconstruction. When a DELL PC with one CPU of 1.7 GHz was used, it took less than five minutes to reconstruct the 3-D image for the demonstrated model imaging problems. This makes the method very attractive to be deployed for clinical application for early detection of cancer. One limitation is that the use of the fluorescence flight time in the calculation of the distance between a tumor and a detector may restrict the estimation of tumor center location to presence of a single tumor. We are developing new means for estimating centers of multiple tumors.
It is worth noticing that the fluorescence signals used in the present article were obtained through numerical simulation based on forward modeling of radiation transfer. In actual experimental studies, the forward modeling is not required. We only need to measure the time-resolved fluorescence signals in the detector arrays surrounding the four side surfaces of the tissue. In this modeling study, an array of 8×8 detectors was placed on each of the side surfaces. In actual measurements, we may reduce the number of detectors without sacrificing the accuracy. For instance, an array of 4×8 detectors may be deployed (each side of the tissue has 8 detectors and all the detectors are connected with a 32-channel streak camera). The array is then moved vertically (e.g., along the Z-direction in Fig. 1) to get fluorescence signals on the whole side surfaces. When an ultrafast laser with repetition rate of 76 MHz is used for excitation, the time gap between two successive short pulses is only 13 ns. The whole measurement can be completed within a very short time period.
We developed a novel 3-D optical fluorescence imaging method for early detection of small tumor. The image reconstruction in this method does not involve any inverse optimization. It is very efficient and can be completed in minutes by a PC. We demonstrated the feasibility of this method in the imaging of a small tumor embedded in a cubical model tissue. The fluorescence signals were obtained via simulation of the coupled transient RTEs that govern the photon migration of both the excitation laser and consequent fluorescence in the tissue. The simulation was set to mimic actual experiment. The uptake of fluorescent agent was considered to be non-uniformly distributed in the tissue and tumor. The reconstructed 3-D image closely matches the actual tumor in terms of shape, size, and location. Numerous tomographic images are obtainable by slicing the 3-D image. The detected images are of high contrast and high accuracy.
Z. Guo Acknowledges the partial support of this research from the New Jersey Space Grant Consortium, the Charles and Johanna Busch Memorial Fund managed at Rutgers University, and the NSF Grant (CTS-0318001).
References and Links
1. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence Lifetime Imaging in Turbid Media,” Opt. Lett. 21, 158–160 (1996). [CrossRef]
2. V. Ntziachristos, A. G. Yodh, M. Schnall, and B. Chance, “Concurrent MRI and Diffuse Optical Tomography of Breast After Indocyanine Green Enhancement,” Proc. Natl. Acad. Sci. USA 97, 2767–2772 (2000). [CrossRef] [PubMed]
3. R. Roy and E. M. Sevick-Muraca, “Truncated Newton’s Optimization Scheme for Absorption and Fluorescence Optical Tomography: Part 1 Theory and Formulation,” Opt. Express 4, 353–371 (1999). [CrossRef] [PubMed]
4. A. Godavarty, M. J. Eppstein, C. Zhang, S. Theru, A. B. Thompson, M. Gurfinkel, and E. M. Sevick-Muraca, “Fluorescence-Enhanced Optical Imaging in Large Tissue Volumes Using a Gain-Modulated ICCD Camera,” Phys. Med. Biol. 48, 1701–1720 (2003). [CrossRef] [PubMed]
5. J. Wu, L. Perelman, R. R. Dasari, and M. S. Feld, “Fluorescence Tomographic Imaging in Turbid Media Using Early-Arriving Photons and Laplace Transforms,” Proc. Natl. Acad. Sci. USA 94, 8783–8788 (1997). [CrossRef] [PubMed]
8. Z. Guo and S. Kumar, “Discrete Ordinates Solution of Short Pulse Laser Transport in Two-Dimensional Turbid Media,” Appl. Opt. 40, 3156–3163 (2001). [CrossRef]
10. Z. Guo and S. Kumar, “Radiation Element Method for Transient Hyperbolic Radiative Transfer in Plane-Parallel Inhomogeneous Media,” Numerical Heat Transfer B 39, 371–387 (2001). [CrossRef]
11. B. Valeur, Molecular Fluorescence: Principles and Applications (Wiley-VCH, Weinheim, New York, 2002).
12. K. R. Castleman, Digital Image Processing (Prentice Hall, New Jersey, USA, 1996).
13. S. K. Wan, Z. Guo, S. Kumar, J. Aber, and B. A. Garetz, “Noninvasive Detection of Inhomogeneities in Turbid Media with Time-Resolved Log-Slope Analysis,” JQSRT 84, 493–500 (2004). [CrossRef]