Download Report

Principal component analysis of dynamic
fluorescence diffuse optical tomography images
Xin Liu, Daifa Wang, Fei Liu, and Jing Bai*
Department of Biomedical Engineering, School of Medicine, Tsinghua University, Beijing 100084, China
*deabj@tsinghua.edu.cn
Abstract: Challenges remain in resolving drug distributions within small
animals utilizing fluorescence diffuse optical tomography (FDOT). In this
paper, we present a new method for detecting and visualizing organs with
different kinetics utilizing principal component analysis (PCA).
Indocynaine green (ICG) metabolic processes are simulated and imaged
using FDOT. When applied to the time series of generated FDOT images,
PCA provides a set of the principal components (PCs) which can represent
spatial patterns associated with different kinetic behavior. Simulation and
experiment studies are both performed to validate the performance of the
proposed algorithm. The results suggest that we are able to extract and
illustrate changes in ICG kinetic behavior between the heart and the lungs.
©2010 Optical Society of America
OCIS codes: (170.6960) Tomography; (170.3010) Image reconstruction techniques; (170.6280)
Spectroscopy, fluorescence and luminescence; (170.3880) Medical and biological imaging.
References and links
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of
whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005).
V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves
protease activity in vivo,” Nat. Med. 8(7), 757–761 (2002).
J. Haller, D. Hyde, N. Deliolanis, R. de Kleine, M. Niedre, and V. Ntziachristos, “Visualization of pulmonary
inflammation using noninvasive fluorescence molecular imaging,” J. Appl. Physiol. 104(3), 795–802 (2008).
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. U.S.A. 97(6), 2767–2772 (2000).
X. Montet, V. Ntziachristos, J. Grimm, and R. Weissleder, “Tomographic fluorescence mapping of tumor
targets,” Cancer Res. 65(14), 6330–6336 (2005).
A. Corlu, R. Choe, T. Durduran, M. A. Rosen, M. Schweiger, S. R. Arridge, M. D. Schnall, and A. G. Yodh,
“Three-dimensional in vivo fluorescence diffuse optical tomography of breast cancer in humans,” Opt. Express
15(11), 6696–6716 (2007).
V. Ntziachristos, E. A. Schellenberger, J. Ripoll, D. Yessayan, E. Graves, A. Bogdanov, Jr., L. Josephson, and R.
Weissleder, “Visualization of antitumor treatment by means of fluorescence molecular tomography with an
annexin V-Cy5.5 conjugate,” Proc. Natl. Acad. Sci. U.S.A. 101(33), 12294–12299 (2004).
R. Weissleder, and V. Ntziachristos, “Shedding light onto live molecular targets,” Nat. Med. 9(1), 123–128
(2003).
X. Michalet, F. F. Pinaud, L. A. Bentolila, J. M. Tsay, S. Doose, J. J. Li, G. Sundaresan, A. M. Wu, S. S.
Gambhir, and S. Weiss, “Quantum dots for live cells, in vivo imaging, and diagnostics,” Science 307(5709),
538–544 (2005).
G. Hu, J. Yao, and J. Bai, “Full-angle optical imaging of near-infrared fluorescent probes implanted in small
animals,” Prog. Nat. Sci. 18(6), 707–711 (2008).
S. V. Patwardhan, S. R. Bloch, S. A. Achilefu, and J. P. Culver, “Time-dependent whole-body fluorescence
tomography of probe bio-distributions in mice,” Opt. Express 13(7), 2564–2577 (2005).
E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution fluorescence molecular
imaging system for small animal imaging,” Med. Phys. 30(5), 901–911 (2003).
A. Soubret, J. Ripoll, and V. Ntziachristos, “Accuracy of fluorescent tomography in the presence of
heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imaging 24(10), 1377–1386 (2005).
X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography
using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007).
A. Joshi, W. Bangerth, and E. M. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence
optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004).
M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of
light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995).
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6300
17. E. M. C. Hillman, and A. Moore, “All-optical anatomical co-registration for molecular imaging of small animals
using dynamic contrast,” Nat. Photonics 1(9), 526–530 (2007).
18. S. Kawata, K. Sasaki, and S. Minami, “Component analysis of spatial and spectral patterns in multispectral
images. I. Basis,” J. Opt. Soc. Am. A 4(11), 2101–2106 (1987).
19. P. Razifar, H. Engler, G. Blomquist, A. Ringheim, S. Estrada, B. Långström, and M. Bergström, “Principal
component analysis with pre-normalization improves the signal-to-noise ratio and image quality in positron
emission tomography studies of amyloid deposits in Alzheimer’s disease,” Phys. Med. Biol. 54(11), 3595–3612
(2009).
20. F. Pedersen, M. Bergströme, E. Bengtsson, and B. Långström, “Principal component analysis of dynamic
positron emission tomography images,” Eur. J. Nucl. Med. 21(12), 1285–1292 (1994).
21. K. J. Friston, C. D. Frith, P. F. Liddle, and R. S. J. Frackowiak, “Functional connectivity: the principalcomponent analysis of large (PET) data sets,” J. Cereb. Blood Flow Metab. 13(1), 5–14 (1993).
22. A. H. Andersen, D. M. Gash, and M. J. Avison, “Principal component analysis of the dynamic response
measured by fMRI: a generalized linear systems framework,” Magn. Reson. Imaging 17(6), 795–815 (1999).
23. C. G. Thomas, R. A. Harshman, and R. S. Menon, “Noise reduction in BOLD-based fMRI using component
analysis,” Neuroimage 17(3), 1521–1537 (2002).
24. B. Dogdas, D. Stout, A. F. Chatziioannou, and R. M. Leahy, “Digimouse: a 3D whole body mouse atlas from CT
and cryosection data,” Phys. Med. Biol. 52(3), 577–587 (2007).
25. D. Stout, P. Chow, R. Silverman, R. M. Leahy, X. Lewis, S. Gambhir, and A. Chatziioannou, “Creating a whole
body digital mouse atlas with PET, CT and cryosection images,” Mol. Imaging Biol. 4, S27 (2002).
26. A. Kak, and M. Slaney, Computerized Tomographic Imaging (New York: IEEE Press, 1987), ch. 7.
27. K. Esbensen, and P. Geladi, “Strategy of multivariate image analysis (MIA),” Chemom. Intell. Lab. Syst. 7(1-2),
67–86 (1989).
28. 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(17), 4225–
4241 (2005).
29. D. Wang, X. Liu, and J. Bai, “Analysis of fast full angle fluorescence diffuse optical tomography with beamforming illumination,” Opt. Express 17(24), 21376–21395 (2009).
30. V. Saxena, M. Sadoqi, and J. Shao, “Polymeric nanoparticulate delivery system for Indocyanine green:
biodistribution in healthy mice,” Int. J. Pharm. 308(1-2), 200–204 (2006).
31. H. Meyer, A. Garofalakis, G. Zacharakis, S. Psycharakis, C. Mamalaki, D. Kioussis, E. N. Economou, V.
Ntziachristos, and J. Ripoll, “Noncontact optical imaging in mice with full angular coverage and automatic
surface extraction,” Appl. Opt. 46(17), 3617–3627 (2007).
32. D. Wang, X. Liu, Y. Chen, and J. Bai, “A novel finite-element-based algorithm for fluorescence molecular
tomography of heterogeneous media,” IEEE Trans. Inf. Technol. Biomed. 13(5), 766–773 (2009).
1. Introduction
Fluorescence diffuse optical tomography (FDOT) is an optical imaging method that allows
non-invasive, quantitative, three-dimensional (3-D) imaging of biological and biochemical
processes in living small animals. At present, FDOT has been successfully applied in
enzymes [1,2], inflammation [3], oncology [4–6], and therapy [7], etc.
With the advances in optical probes [8,9], imaging systems [8,10,11], and reconstruction
algorithms [12–16], tomographic imaging of dynamic biological activities in small animal
body in vivo is now possible. Compared with static FDOT image reconstructed at a specific
time point, the time series of tomographic imaging impart the ability to capture the complete
dynamic course of absorption, distribution, and elimination of the fluorophores, contrast
agents, or drug, by adding time as a new dimension. We call this time series dynamic
fluorescence diffuse optical tomography (D-FDOT). By allowing longitudinal time-lapse
visualization of metabolic processes of contrast agents, D-FDOT provides an attractive
approach in studying drug delivery, tumor detection, and treatment monitoring in small
animals in vivo.
However, challenges remain in resolving organs and functional structures with different
kinetics in small animal body in vivo. The main problem is that the delivery route of drug
(non-specific contrast agents) is complex. They may distribute all over the body through
circulatory system, instead of accumulating at some specific place, which will result in poor
imaging contrast and specificity. Together, the highly diffusive nature of photon migration in
tissues and ill-posed characteristic of the inverse problem are responsible for the poor image
quality and low spatial resolution. Hence, it is difficult to detect and visualize changes in
kinetic behavior between major organs of small animals directly from FDOT images obtained
at a specific time point. Here, instead of focusing on static tomographic image, we look at the
dynamics of drug during FDOT imaging.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6301
Each organ in the body plays a different role in the course of drug circulating,
accumulating or metabolizing, and so each organ will exhibit a distinctive time course in its
optical signal [17]. An exposition of resolving different functional structures based on DFDOT images reduces to an examination of its correlation structure. Here, the correlation
feature refers to the correlative kinetic behavior obtained over time in the same organ.
Principal component analysis (PCA) is very suitable for this examination. PCA provides the
ability of extracting the important features of the correlation matrix in terms of principal
components or eigenvectors [18]. Each principal component (PC) represents a major trend of
kinetic behavior within which there are high inter-correlations. At present, PCA has been
widely applied in analyzing dynamic positron emission tomography (PET) [19–21] and
functional magnetic resonance imaging (fMRI) data [22,23]. In 2007, Hillman and Moore
[17] firstly demonstrated the capability of applying PCA in analyzing dynamic optical
projections. Using a CCD camera to image the dynamics metabolic processes of optical
contrast agents injected through tail veins, they mainly focused on decomposing the 2-D
fluorescence reflectance imaging and extracted anatomical information of various internal
organs utilizing PCA. Despite of these successful applications of PCA, no trial has been
reported for PCA in 3-D D-FDOT images. As PCA is a data-driven method and 3-D D-FDOT
has its own characteristics of poor spatial resolution and diffusive nature, investigation on the
applying PCA-based methods in 3-D D-FDOT is very important.
In this paper, we propose to apply a PCA-based method on a time series of 3-D FDOT
images to resolve organs and functional structures with different kinetic patterns. The contrast
agent we employed, indocyanine green (ICG), is a widely used optical fluorophores in the
near-infrared (NIR) range for small animal imaging [11], as well as for human breast cancer
imaging [4]. The proposed method was validated by simulation and experiment studies on a
non-contact, full angle FDOT imaging system [10]. In simulation studies, we simulated a case
where ICG distributed through the heart and the lungs of mouse after tail vein injection, based
on a 3-D digital mouse model [24,25]. In experiment studies, we repeated the simulated
conditions by two tubes filled with different concentrations of ICG at different time points.
The results suggest that we are able to detect and visualize regions with different kinetics
based on D-FDOT images reconstructed much better.
The outline of this paper is as follows. In section 2, the methods used are detailed. In
section 3, the simulation and experiment are shown. In section 4, the results are described.
Finally, we discuss and conclude the major findings of this study in section 5.
2. Methods
2.1. Forward and inverse problems
In highly scattering tissue medium, the photon migration in biological tissues can be modeled
using the coupled diffusion equation with the Robin-type boundary condition [16]. In FDOT,
the Green's function G (rs , r ) due to a continuous wave (CW) point excitation source
δ (r − rs ) , describing the photon density at position r, can be obtained as follows,
∇ ⋅ [ D(r )∇G (rs , r )] − µa (r )G (rs , r ) = −δ (r − rs )

∂G (rs , r )

 2qD(r ) ∂n + G (rs , r ) = 0

r ∈Ω
r ∈ ∂Ω
,
(1)
where Ω is the domain of the imaged object and ∂Ω is the boundary; µa (r ) is the absorption
coefficient and D(r ) = 1 / (3µ s' (r )) is the diffusion coefficient of the tissue with the reduced
scattering coefficient µ s' (r ) at position r; n denotes the outward normal vector to the
boundary and q is a constant depending upon the optical reflective index mismatch at the
boundary. A collimated point source spot can be modeled as an isotropic source δ (r − rs ) ,
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6302
where rs is the point one transport mean free path ltr = 1 / µ s' (r ) into the medium from the
illumination spot [16]. The diffusion equation can be solved using the Galerkin finite element
method to obtain the Green's functions [14,15]. After that, the inverse problem is generated
based on Normalized Born approximation [13], which reduces the influence of mouse tissue
heterogeneity. The normalized ratio between the measured fluorescence signal Φ m (rd ) and
the corresponding excitation signal Φ x (rd ) at detector point rd is given as follows,
G (rd , rp )G (rs , rp )n(rp )
Φ m (rd )
drp ,
= Θ ∫V
G (rs , rd )
Φ x (rd )
(2)
where rp is the point inside the volume V considered for reconstruction; n(rp ) denotes the
fluorescence yield at point rp , which is directly proportional to the fluorophores's
concentration; G (rd , rp ) denotes the Green's function value at rp due to a point source
δ (r − rd ) at rd ; Strictly speaking, G (rs , rp ) describes the light propagation at the excitation
wavelength and G (rd , rp ) corresponds to light propagation at the emission wavelength. The
excitation and emission wavelength are close to each other. In addition, the Normalized Born
method can largely reduce the effects of incorrect assumption of optical properties. For the
moment, we assume that the Green's functions are similar at these two wavelengths for
simplification. Θ is a calibration factor which accounts for the unknown gain and attenuation
factors of the system, such as excitation light power. To solve the inverse problem, Eq. (2) is
discretized on a 3-D grid of size N x , N y , N z . For each single source-detector pair (rs1 , rd 1 )
and each point (rp ) p =1... N
xN y Nz
on the grid, Eq. (2) is described as follows,
 n(r1 ) 


 . 
Φ m (rd 1 )
Θ∆V
=
[G (rd 1 , r1 )G (r1 , rs1 ),..., G (rdN , rN )G (rN , rsN )] ⋅  .  ,
Φ x (rd 1 ) G (rs1 , rd 1 )


 . 
 n( r ) 
 N 
(3)
where N is the total number of the descretized voxels. For M source-detector pairs, a linear
system is generated as,
 Φ m (rd 1 ) / Φ x (rd 1 )   W11

  .
.

 

= .
.

 
.

  .
Φ m (rdM ) / Φ x (rdM )  WM 1
.
.
.
.
.
.
W1N   n(r1 ) 
.   . 
.  ⋅  . .
 

.
.   . 
. WMN   n(rN ) 
.
(4)
where Wij is an element of the weight matrix W connecting the measurement to the unknown
fluorescence yield. The unknown n(rp ) in each voxel is obtained by solving Eq. (4) using
algebraic reconstruction technique (ART) [26] with non-negative constraints.
2.2. Principal component analysis
In kinetic analysis, the drug distribution within biological body is sequentially imaged at
different time points. These sequential images can be regarded as multivariate images from
which physiological, biochemical and functional information can be derived. Principal
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6303
component analysis (PCA), as a common multivariate image analysis method, provides a
graceful strategy for resolving different kinetic behavior of the signal in each pixel by
remapping the data into a new coordinate system.
In this study, the input data can be denoted as X = { X 1 , X 2 ,..., X t ,..., X N } , where N is the
number of frames and X t is the frame taken at a given time t. Here, image frame denotes the
reconstructed tomographic image from the same slice but from different time points. To apply
PCA to the spatial and temporal FDOT images at the same slice, we assume that the variables
are the frames taken in a time sequence, and a pixel's temporal distribution is an observation
of these variables. If the reconstructed image has M = R × C pixels, then we can represent
one frame data as X tT = ( x1,t , x2,t ,..., xk ,t ,..., xM ,t ) , where xk ,t is the density of the pixel k at the
frame t (T means transposition operation), and the size of the image sequence X is of M × N
matrix. Then, we define the matrix, S , as follows,
S=
1
XT X,
M −1
(5)
where X is a data set with zero mean obtained from the original one by subtracting its mean
value from each frame. Mathematically, principal component expansion is obtained by the
diagonalization of matrix S , which produces a set of eigenvalues, λi , and eigenvectors, Ei .
If we define the set of eigenvectors as E = {E1 , E2 ,..., Ei ,..., EN } , where the i column contains
the elements of the vector EiT = {e1,i , e2,i ,..., e j ,i ,..., EN ,i } , the principal components of the
image sequence can be obtained as a M × N matrix, P , as follows,
P = XE,
(6)
and each principal components (PCs) is given below,
N
Pi = ∑ et ,i X t .
(7)
t
The condition of Cov ( Pj , Pk ) = 0, j ≠ k is required, meaning that each principal
component is functionally uncorrelated to each other. Further each element of principal
component is utilized as a weight factor for creating image, termed PC image [20,27] that
illustrates organs and functional structures with different kinetic behavior.
To sum up, the corresponding flowchart of the proposed algorithm is shown in Fig. 1.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6304
Fig. 1. The flowchart of the proposed algorithm.
3. Materials
In validating the performance of PCA on D-FDOT images, we performed simulation and
experiment studies based on a non-contact, full angle FDOT system. The sketch of the system
is shown in Fig. 2, which is similar to that in [10].
Fig. 2. The schematic of the imaging system.
3.1. Simulation studies
3.1.1. Setup for simulation studies
In simulation studies, a virtual mouse atlas was employed to provide 3-D anatomical
information [24,25]. To simplify the experimental design, we only focused on imaging kinetic
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6305
behavior inside mouse chest region instead of whole body. Therefore, only the mouse torso
from the neck to the base of the lungs was simulated, with a length of 1.6 cm.
The imaged mouse was suspended on a rotating stage, as shown in Fig. 2. The rotational
axis of the mouse was defined as the Z axis with the bottom plane set as z = 0 cm. According
to [28], optical parameters (absorption, scattering) were assigned to the heart and lungs to
simulate photons propagation in biological tissues. The optical properties outside heart and
lungs were regarded as homogeneous. Detailed information about optical properties of
biological tissues in mouse is presented in Table 1. Based on [29], in this simulations,
fluorescence tomography of 360o full view was performed with 18 projections at every 20o
with the excitation source at height 0.8 cm.
Table 1. Optical parameters of biological tissues in mouse at 700-800 nma
µ a (cm -1 )
Tissue type
Heart
Lungs
0.156
0.516
µ 's (cm -1 )
9.0
21.2
a
From Phys. Med. Biol. 50, 4225 (2005).
3.1.2. Dynamic modeling in simulation studies
To optimally evaluate the performance of the proposed method, we simulated the metabolic
processes of ICG within small animals. Firstly, we set a series of fluorescent yields
(concentrations of ICG) to the heart and the lungs at six time points (5 min, 10 min, 15 min,
30 min, 60 min, and 120 min) according to [30]. The anatomical information of the heart and
the lungs is described in Fig. 3(a) and the concentration time course of ICG is plotted in
Fig. 3(b). Secondly, six tomographic images were reconstructed to depict ICG distribution
within the body at corresponding time points. Finally, we assembled these tomographic
images to form D-FDOT images which were used to simulate dynamic course of ICG
kinetics.
Fig. 3. Setup for simulation studies. (a) The mouse geometry model used in simulation studies.
The gray part in (a) depicts the mouse surface with a length of 1.6 cm from the neck to the
base of the lungs. The red part in (a) depicts anatomical information of the heart and the green
part in (a) depicts anatomical information of the lungs. In order to reduce the boundary
artifacts which interfere with the finite element computation, the model is generated by
sampling the original atlas data (intersections of the vertical and horizontal lines) and then
approximating the curves using spline function to form the torso surface. (b) ICG
concentration time course in the heart and the lungs after tail vein injection. The circles in (b)
depict actual concentration value at corresponding time points according to [30]. Different
colors correspond to different time course curves (red:heart; green:lungs).
It should be noted that, the uptake of ICG in organ was assumed to be approximately
uniformly distributed rather than accumulated in limited regions of organ. Herein, we
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6306
performed simulations assuming that the whole region of organ was tagged with ICG. In
addition, considering that ICG was a kind of non-specific fluorescence contrast agent, the
heart and the lungs were labeled simultaneously in this study.
3.1.3. Reconstructions of synthetic data
Reconstructions were performed to obtain tomographic images at different time points. The
excitation and emission data were synthesized using finite element method based on Eq. (2).
The geometrical model in Fig. 3(a) was discretized into 10,987 nodes and 55,398 tetrahedral
elements. The detectors were located on the boundary finite element nodes which were
between 1.4 cm height range and inside 150o FOV corresponding to each point source. The
volume considered for reconstruction was 2.0 cm × 2.5 cm × 1.6 cm and sampled to
32 × 40 × 32 voxels. Only the 16,518 voxels inside the imaged object were considered for
reconstruction and the reconstruction was terminated after 100 ART iterations.
3.1.4. Principal component analysis for synthetic D-FDOT images
To extract ICG kinetic properties in the heart and the lungs from D-FDOT images, PCA was
performed. Here, the input data was six frames reconstructed tomographic images, having
zero mean, from D-FDOT images at the same slice. Each frame with a spatial size of 32 × 40
pixels was reshaped as a column vector containing 1,280 elements. By the transformation, the
set of input data was placed in a matrix with 1,280 rows and 6 columns. After applying
Eq. (6), we obtained six principal components which were functionally uncorrelated to each
other. Each positive and negative element of principal component was respectively utilized as
a weight factor for creating image that illustrated regions with different kinetic behavior using
different colors. Red color depicted negative PC images and green color depicted positive PC
images.
3.2. Experiment studies
3.2.1. Setup for experiment studies
A physical experiment was performed to further validate the performance of the proposed
method. As shown in Fig. 4(a), two transparent glass tubes (0.3 cm diameter, 0.6 cm length)
filled with different concentrations of ICG were immerged in a cylinder phantom. The
phantom was made of a glass cylinder (3.0 cm diameter) filled with 1% intralipid
( µ s' = 10.0 cm -1 ,µa =0.02 cm -1 ). The two tubes were separated with an edge-to-edge distance
of 0.2 cm along Y axis.
In the experiments, the imaged phantom was placed on a rotating stage which allowed
rotation and shift of the target around Z axis for collecting projections evenly distributed over
360o , as shown in Fig. 2. The small light spot from a 250 W Halogen lamp traveled through a
775 ± 23 nm band-pass filter and focused on the surface of imaged object. A 512 × 512
element, -70o C cooled EMCCD array coupled with a Nikkor 60 mm f/2.8D lens was placed
on the opposite side of imaged object, collecting the photons propagating through the imaged
object. When collecting the fluorescence images, a 840 ± 6 nm band-pass filter was placed in
front of EMCCD camera. When collecting the excitation light images, a neutral density filter
of 1% transmittance was used. When collecting the white light images, a white light bulb was
employed to replace the excitation light. The white light images at different projection angles
were acquired for recovering the 3-D geometry of the imaged object, which was necessary for
modeling diffuse light transportation.
3.2.2. Dynamic modeling in experiment studies
To simulate the time active curve of ICG in the heart and the lungs, the concentrations of ICG
in tube 1 and tube 2 at different time points were set according to Fig. 4(b). Tube 1 depicted
ICG metabolic processes in the heart and tube 2 depicted ICG metabolic processes in the
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6307
lungs. For each experiment at specific point (frame), 24 excitation and emission images were
collected at every 15o using EMCCD with the excitation source at height 2.9 cm. The total
power of the point source was about 20 mW. 72 white light images acquired in 5o steps were
back-projected to form a 3-D geometry of phantom [31].
Fig. 4. Setup for experiment studies. (a) Experiment setup. Two glass tubes (diameter of 0.3
cm and height of 0.6 cm) filled with different concentrations of ICG were placed inside a
cylinder phantom (a glass cylinder of 3.0 cm diameter filled with 1% intralipid). The edge to
edge distance along Y axis was 0.2 cm. (b) ICG concentration time course in tube 1 and tube 2.
Tube 1 simulated ICG metabolic processes in the heart and tube 2 simulated ICG metabolic
processes in the lungs. The circles in (b) depicted actual concentration value at corresponding
time points. Different colors corresponded to different time course curves (red:tube 1;
green:tube 2).
3.2.3. Reconstructions of experimental data
The cylinder phantom in Fig. 4(a) was discretized into 17,558 nodes and 78,782 tetrahedral
elements. The detectors were located on the boundary finite element nodes which were
between 3.2 cm height range and inside 150o FOV corresponding to each point source. The
volume considered for reconstruction was 3.0 cm × 3.0 cm × 5.0 cm and sampled to
30 × 30 × 50 voxels. The 27,178 voxels inside the imaged object were used for reconstruction.
Reconstructions were terminated after 100 ART iterations.
3.2.4. Principal component analysis for experimental D-FDOT images
The main description was the same as in section 3.1.4. Briefly, the input data was 6 frames
reconstructed tomographic imaging with a spatial size of 30 × 30 pixels and the input matrix
had 900 rows and 6 columns.
4. Results
4.1. Reconstructions of synthetic data
Figure 5 shows a slice through the chest, including the heart and the lungs, in the mouse after
ICG injection. It is a study of ICG metabolic processes in the heart and the lungs. The six
reconstructed images (5 min, 10 min, 15 min, 30 min, 60 min, and 120 min) from D-FDOT
images at the same height (z = 0.8 cm) are described in Figs. 5(a) to 5(f). The results suggest
that it is quite difficult to resolve the uptake of ICG in the heart and the lungs directly from
static tomographic image.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6308
Fig. 5. Reconstruction of synthetic data from a dynamic study of ICG metabolic processes in
the heart and the lungs. The images are at z = 0.8 cm. (a)-(f) The reconstructed results at 5 min,
10 min, 15 min, 30 min, 60 min, and 120 min. Different colors correspond to actual boundary
of different organs (red:heart; green:lungs; yellow:surface).
4.2. Principal component analysis for synthetic D- FDOT images
As shown in Fig. 6, when PCA was applied to D-FDOT images at specific slice (z = 0.8 cm)
shown in Fig. 5, we obtained six PC images reflecting different kinetic behavior. In Figs. 6(a)
to 6(f), the obtained positive PC images are depicted using green color. In Figs. 6(g) to 6(l),
the obtained negative PC images are depicted using red color. Based on the anatomical
information of the heart and the lungs, we find that a very high uptake of ICG in the lungs is
nicely demonstrated using the positive PC2 image, as shown in Fig. 6(b). Similarly, in
Fig. 6(h), a very high uptake of ICG in the heart is demonstrated using the negative PC2
image. Higher PC images (PC4 to PC6) are similar to the PC3 image, containing no structure.
The results suggest that PCA is a useful tool to detect and visualize organs and functional
structures with different kinetics. In addition, it is possible to give a simple interpretation that
PCA is also seen as data reduction method since only the first two PC images contain
structure.
Fig. 6. The PC images obtained when PCA was applied to D-FDOT images shown in Fig. 5.
(a)-(f) The six positive PC images. (g)-(l) The six negative PC images. In the PC2 image, the
uptake of ICG in the heart (negative) and the lungs (positive) was indicated. Different colors
corresponded to actual boundary of different organs (red:heart; green:lungs; yellow:surface).
To demonstrate the performance of PCA on D-FDOT images at all slices, we describe the
3-D visualization results of the negative and the positive PC2 images, by extracting
isosurfaces from all obtained PC2 images, as shown in Figs. 7(a) and 7(d). Figure 7(a)
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6309
indicates the 3-D region of ICG distribution in the heart. Figure 7(d) indicates the 3-D region
of ICG distribution in the lungs. In the 3-D visualization of the PC2 images, reconstructed
fluorescence signals, which was less than 2% of the maximum signal intensity, were not
considered. These fluorescence signals did not contain information but induced artifacts into
the reconstructed images. A visual comparison of the PC2 images in Figs. 7(a) and 7(d) to
anatomical information in Figs. 7(b) and 7(e) are shown in Figs. 7(c) and 7(f). The red parts
in Figs. 7(c) and 7(f) depict the 3-D anatomical information of the heart and the lungs. The
green parts in Figs. 7(c) and 7(f) depict the 3-D visualization results of the negative and the
positive PC2 images.
Fig. 7. The 3-D visualization results of the PC2 images. (a) and (d) The 3-D visualization
results of the negative and the positive PC2 images obtained when PCA was applied to the DFDOT images at all slices. (a) indicated the uptake of ICG in the heart. (d) indicated the uptake
of ICG in the lungs. (b) and (e) The 3-D anatomical information of the heart and the lungs in
the mouse. (c) and (f) A visual comparison of the PC2 images in (a) and (d) to anatomical
information in (b) and (e). The red parts in (c) and (f) indicated the 3-D anatomical information
of the heart and the lungs. The green parts in (c) and (f) indicated the 3-D visualization results
of the negative and the positive PC2 images.
Figure 8 shows a pictogram of how principal component analysis works with D-FDOT
images shown in Fig. 5. Each frame in Fig. 5, after subtraction of its mean value, was
arranged as a columns vector of input matrix X . After applying Eq. (6), six functionally
uncorrelated principal components were generated and arranged in the matrix P .
Subsequently, each positive and negative element of principal component was respectively
utilized as a weight factor for creating image that illustrated functional structures with
different kinetic behavior. The relation between matrices X and P was given by the matrix
E.
Comparison between 3-D reconstructed results at a specific time point (5 min) and the
merged results of the negative and positive PC2 images, we further demonstrated that organs
and functional structures with different kinetic patterns can be resolved when PCA is applied
to a time series of FDOT images.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6310
Fig. 8. Pictogram of how principal component expansion works with D-FDOT images. In this
plot, we show how the reconstructed results
( X ,X ,...,X ) in Fig. 5 (after subtraction of
1
2
6
the mean value) were transformed into six principal components
(PC1,PC2,...,PC6) that
were presented as images through the matrix E . The red part in the PC2 image indicated the
uptake of ICG in the heart and the green part in the PC2 image indicated the uptake of ICG in
the lungs. The black curves in 3-D view depicted the height of selected dynamic tomographic
images. The 3-D view (left) depicted the 3-D reconstructed results of frame 1 (5 min). The 3-D
view (right) depicted the merged results of the positive and the negative PC2 images. The red
part in 3-D view (right) depicted the 3-D visualization results of the negative PC2 images in
Fig. 7(a). The green part in 3-D view (right) depicted the 3-D visualization results of the
positive PC2 images in Fig. 7(d).
4.3. Reconstructions of experimental data
Figure 9 shows reconstructed results of experimental data from the same slice but from
different frames. The six reconstructed images at the same height (z = 3.0 cm) are shown in
Figs. 9(a) to 9(f). Similar as in section 4.1, it is difficult to resolve the uptake of ICG in two
tubes based on each static tomographic image.
Fig. 9. Reconstruction of experimental data at different frames. The images are at z = 3.0 cm.
(a)-(f) Reconstruction results at frames 1 to 6. The red curve on the cross images depicts the
phantom boundary, and the black circles depict the actual tubes. All images are displayed at
the same range.
4.4. Principal component analysis for experimental D-FDOT images
The results in section 4.2 suggest that function structures with different dynamic behavior can
be resolved utilizing PCA in simulation studies. As expected, in experiment studies, when
PCA was applied to D-FDOT images shown in Fig. 9, we detected and visualized the uptake
of ICG in two tubes with an edge-to-edge distance of 0.2 cm. As shown in Figs. 10(b) and
10(c), the uptakes of ICG in tube 1 and tube 2 are illustrated using the negative and the
positive PC3 image, respectively. The two parts indicate the kinetics of ICG in the heart and
the lungs. When PCA was applied to D-FDOT images at all slices, we generated the 3-D
visualization results of the positive and the negative PC3 images, as shown in Figs. 10(e) and
10(f). In the 3-D visualization, reconstructed fluorescence signals, which was less than 2% of
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6311
the maximum signal intensity, were not considered. Comparing the reconstructed result of
frame 3 in Figs. 10(a) with the PC3 images in Figs. 10(b) and 10(c), we validated the
capability of PCA in illustrating changes in kinetic behavior in experiment studies.
Fig. 10. Comparison of reconstructed result to PCA-based result. (a) The 2-D reconstructed
result of frame 3. The cross section image is at z = 3.0 cm, which is depicted by the red curve
in (d). The red circle in (a) depicts the phantom boundary, and the black circles depict the
actual tubes. (b) and (c) The positive and the negative PC3 images obtained when PCA was
applied to D-FDOT images shown in Fig. 9. The red circles in (b) and (c) depict the phantom
boundary, and the cyan circles in (b) and (c) depict the actual tubes. (d) The 3-D reconstructed
result of frame 3. (e) and (f) The 3-D visualization results of the positive and the negative PC3
images obtained when PCA was applied to the D-FDOT images at all slices.
Figure 11 further shows differences in the uptake of ICG in two tubes using the merged
results of the positive and the negative PC3 image. Figure 11(a) depicts the 2-D merged result
of the positive PC3 image in Fig. 10(b) and the negative PC3 image in Fig. 10(c).
Figures 11(b) and 11(c) depict the 3-D merged results of the positive PC3 images in
Fig. 10(e) and the negative PC3 images in Fig. 10(f), using different views. The red part in
Fig. 11 indicates the uptake of ICG in tube 1 that is used to simulate ICG metabolic processes
in the heart. The green part in Fig. 11 indicates the uptake of ICG in tube 2 that is used to
simulate ICG metabolic processes in the lungs.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6312
Fig. 11. The merged results of the positive and the negative PC3 images. (a) The 2-D merged
results of the positive PC3 image in Fig. 10(b) and the negative PC3 image in Fig. 10(c). The
cross section image is at z = 3.0 cm, which is depicted by the red curves in (b) and (c). The red
circle in (a) depicts the phantom boundary, and the cyan circles depict the actual tubes. (b) and
(c) The 3-D merged results of the positive PC3 images in Fig. 10(e) and the negative PC3
images in Fig. 10(f) using different views. The red parts indicate the uptake of ICG in tube 1
that is used to simulate ICG metabolic processes in the heart. The green parts indicate the
uptake of ICG in tube 2 that is used to simulate ICG metabolic processes in the lungs.
5. Discussion and conclusion
Fluorescence diffuse optical tomography plays an important role in drug research for allowing
non-invasive, quantitative, 3-D imaging of biological and biochemical processes in living
subjects. However, the complex dynamic behavior of drug often makes it difficult to resolve
organs with different kinetics in small animal body. In this paper, we proposed a PCA-based
method for kinetic studies, to aid in the identification of functional structures with different
kinetic patterns, and then evaluated its performance using simulation and experiment studies.
It was clearly seen from the reconstructed images at different time points, ICG kinetic
behavior in the heart and the lungs (simulation studies) and in the two tubes with an edge-toedge distance of 0.2 cm (experiment studies) could not be resolved. This was probably caused
by the diffusive nature of photon migration in biological tissues and the ill-posed
characteristic of the inverse problem. In contrast, when PCA was applied to the time series of
tomographic images, we generated a set of PC images that illustrate organs and functional
structures with different kinetic patterns. For example, in simulation studies, the uptake of
ICG in the heart and the lungs were illustrated using the negative and the positive PC2
images. Similar in experiment studies, we resolved the spatial structure of two tubes
associated with different kinetic behavior using the negative and the positive PC3 images.
Based on the satisfying results, we believe that PCA-based method provides an attractive
approach in kinetic study. Firstly, this method allows to non-invasive measure changes in the
function of organs. The detecting capability is critical to evaluate the effects of a drug on the
organs in vivo. Secondly, this method has the potential in enhancing the resolution and
specificity of FDOT by emphasizing differences in kinetic behavior. Thirdly, utilizing the
technique, additional anatomical information of organs may be achieved without resorting to
other imaging equipment. Finally, the proposed method is totally independent of any kinetic
model and thus not subject to model-based restrictions.
In this article, the ICG concentration was kept constant while collecting one frame of
fluorescence projection images. It is mainly for simplification. However, it is reasonable. In
experiments, the full angle imaging time of collecting one frame typically takes
approximately 1.3 minutes to 2.0 minutes, depending on the sum of rotation time and total
exposure time of images. During this imaging course, the change of ICG concentration can be
seen as linear approximately. Then, the fluorescence projections of the same frame can be
corrected using linear model, which will lead to the corrected projection images similar to the
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6313
projection images used in this article. Of course, in dynamic FDOT, good data correction
methods for fluorescence projections can improve the reconstructed image quality a lot.
Systematic studies in this field will also be investigated in our future work.
All the obtained results are based on a basic assumption that accurate optical properties
are known. Tissue optical properties are important for FDOT reconstruction and wrong
optical parameters may lead to significant errors to the reconstructed results [32]. It will
directly affect the performance of PCA identifying regions with different kinetics. Moreover,
it should be noted that PCA has difficulties in separating the signal from the noise when the
magnitude of the noise is relatively high [19,23]. Ultimately, it is also difficult to illustrate
how to choose the optimal PC images reflecting different kinetics since the obtained PC
images may be influenced by the reconstructed errors, noise, and time courses of drug, etc.
Hence, in in vivo experiments, many practical factors need further consideration. It is out of
the scope of this article, and will be analyzed in our future work.
In conclusion, using simulation and experiment studies, we have demonstrated the
capabilities of PCA in enhancing FDOT performances in imaging kinetics. This analysis can
be used as a starting point for dynamic studies, e.g. tumor detection, disease progression, and
drug delivery, etc. Future works will be focused applying the method in resolving ICG
distribution within small animal in vivo.
Acknowledgments
This work is supported by the National Natural Science Foundation of China under Grant No.
30670577, 60831003, 30930092, 30872633; the Tsinghua-Yue-Yuen Medical Science
Foundation; the National Basic Research Program of China (973) under Grant No.
2006CB705700; the National High-Tech Research and Development Program of China (863)
under Grant No. 2006AA020803.
#121176 - $15.00 USD
(C) 2010 OSA
Received 10 Dec 2009; revised 6 Feb 2010; accepted 10 Mar 2010; published 12 Mar 2010
15 March 2010 / Vol. 18, No. 6 / OPTICS EXPRESS 6314