1. Introduction

In the last few years, fluorescence tomography (FT) has become a promising technique with many biological and medical applications, for it can produce three-dimensional (3D) visualization, localization, and quantification of fluorescent probes reporting of cancers or tissue abnormalities [1–4]. Since the appearance of FT, a lot of improvements have been obtained in its technological fundamentals and preclinical applications. On one side, many advances have been made in optical probes [5], imaging systems [6–8], and reconstruction algorithms [9–12]. On the other side, FT has developed from single modality to multimodality [13, 14], from static to dynamic imaging [15], and from single wavelength to multi-spectrum [16–18].

Much research has been done on the multispectral fluorescence tomography (MFT), including spectrally resolved emission and excitation. The spectrally resolved emission method uses spectrally-resolved boundary measurements for reconstruction of spatial source distributions, and has been applied in bioluminescence tomography (BLT) [19]. In BLT, multispectral images of a spectrally broad emitting luciferase reporter source are taken on the tissue surface for 3D image reconstruction. The spectrally resolved excitation method uses excitation light of different wavelengths to illuminate fluorophores and generates independent boundary data, i.e. two-dimensional (2D) surface images, at the emission wavelength for tomographic reconstruction. Multispectral excitation-resolved fluorescence tomography (MEFT) has also been reported and confirmed experimentally [16, 17]. The feasibility of excitation-resolved FT of quantum dots has been shown in a simulation study [16]. 3D distribution of fluorescence sources inside tissue can be reconstructed from spectrally varying illumination alone, which has been confirmed by numerical and experimental studies [17].

Different multispectral imaging approaches, including excitation and/or emission resolved methods are compared [18]. It is found that excitation-resolved FT leads to more accurate image reconstructions than multispectral emission detection. Therefore, MEFT is chosen in this study. Compared with common FT reconstructed image excited by a specific excitation light wavelength, MEFT can capture more complete spectrum courses of absorption, distribution, and emission of the fluorophores by adding spectrum as a new dimension. In the dimension of spectrum, the structures containing fluorophores of different concentrations show different spectrum behaviors.

An exposition of separating different structures based on MEFT images comes down to an analysis of the correlative spectrum data obtained over spectrums for fluorophores with different concentrations. As a popular multivariate statistical technique, principal component analysis (PCA) can be used to extract information from the data. In order to achieve these goals, PCA computes new variables called principal components (PCs), which are linear combinations of the original variables [20]. PCA has been widely applied to computer tomography (CT) [21, 22], positron emission tomography (PET) [23, 24] and magnetic resonance imaging (MRI) [25]. In optical imaging, PCA was used to analyze dynamic optical projection data [26]. In that study, optical contrast agents were injected to mice through tail veins. The 2-D fluorescence reflectance data were obtained to dissolve the anatomical structure of various internal organs with PCA. PCA was also used to decomposing the 3D fluorescence tomography images of time series in simulation and phantom studies, but has not been confirmed in vivo [27]. Despite of these successful applications of PCA, no trials have been reported for PCA in 3D MEFT images. As PCA is a data-driven method and 3D MEFT has its own characteristics of poor spatial resolution and diffusive nature, investigation on the application of PCA in 3D MEFT is very important. PCA-based MEFT method provides an attractive approach in studying drug delivery, tumor detection, and treatment monitoring in small animals in vivo.

Challenges remain in resolving functional structures (such as organs) with different kinetics in small animal body in vivo. Each functional structure in the small animal body plays a different role in the course of fluorophore (or drug) circulation, accumulation or metabolism, and so each functional structure shows a specific concentration of fluorophores [28, 29]. In this paper, we propose to apply PCA to a spectral series of 3D MEFT frames in order to separate functional structures containing fluorophores of different concentrations. The proposed method is validated by simulations and experiments with a non-contact, full-angle FT imaging system. In the simulations, a cylinder model containing two tubes with fluorophores of different concentration pairs is used to mimic distributions of fluorophores in the heart and the lungs of a mouse after tail vein injection. Then a digital mouse model with adjacent organs heart and lungs was used for further test. In the experiments, a cylinder with one pair of fluorophore concentrations was used. Indocyanine green (ICG) is employed as the fluorescent dye which is widely used in the near-infrared (NIR) range for small animal imaging. The results of simulations and experiments suggest that the method of performing PCA on MEFT images is capable of detecting and visualizing functional structures with different concentrations.

The outline of this paper is as follows. In Section 2, the basic theories of the proposed algorithm are detailed. In Section 3, the methods of simulation and experiment studies are shown. In Section 4, the results are described. Finally, we discuss and conclude the major findings of this study in Section 5.

2. Theory

2.1. Forward and inverse problems

In the FT system, the propagation of near-infrared light in tissues consists of two main parts: 1) excitation of fluorophores (e.g., fluorescent dyes, fluorescent proteins, or quantum dots) using excitation light, 2) subsequent emission of fluorescence light. Different radiative transfer models are developed to obtain an accurate solution. The radiative transfer models illustrate the propagation of light which takes into account the scattering and absorption properties of the light. Based on the theory of diffusion approximation, the course is often modeled by the diffusion equation (DE) when sources and detectors are separated by more than a few mean free path lengths. A coupled diffusion equation is used to illuminate the model of continuous-wave, where Eqs. (1) and (2) describe the propagation of source light and the emission fluorescence light propagation in tissues, respectively [30]:

The $n(r,{\lambda }_e)$ equals to $\eta {\mu }_{af}(r,{\lambda }_e)$, where $\eta$ is the quantum field, and ${\mu }_{af}$ is the absorption coefficient of the fluorophore which decides the amount of absorbed light. Meanwhile, ${\mu }_{af}(r,{\lambda }_e)=2.303\epsilon ({\lambda }_e,C(r))C(r)$, where $\epsilon$ is the extinction coefficient, $C$ is the concentration distribution of the fluorophore. So $n$ is product of $\eta$, $\epsilon$ and $C$. Quantum field $\eta$ is the intrinsic property of the fluorophore and defined as the ratio of the number of photons emitted to the number of photons absorbed, so $\eta$ is often taken as constant for a given fluorophore and environment. The extinction coefficient $\epsilon$ is decided by many factors, such as the kind of fluorophores, solutions of fluorophores, spectrum of illumination light, and the temperature of environment.

Much research has been done on the relationship between extinction coefficient $\epsilon$ and the concentration $C$ of fluorophore, i.e., on the matter to what extent fluorophore obeys the lambert-beer’s law in different solvents. In blood, a strict observance of this law by any dye is not expected, because blood is an inhomogeneous medium containing high concentration of light scattering particles [31]. On the observation of ICG in water or plasma, different results are shown. Some find that the extinction coefficient keeps constant in a small range (up to 2.5mg/L) or not at all [32], either in water or in plasma. Some others find that ICG does not follow the lambert-beer’s law above 15 mg/L in plasma [33]. In common practice, when single spectrum of light is used as the excitation light, $\epsilon$ is often taken as constant which means $n$ is proportional to $C$, when the $C$ is very low. In MEFT, $\epsilon$ varies with spectrum of excitation light ${\lambda }_e$, which means for a given $C$, $n$ has its specific variation trend with spectrum of excitation light. For different $C$, $n$ shows different variation trends with spectrum of excitation light. That is the foundation of PCA in this study.

2.2. Principal component analysis

PCA is a widely used technique for multivariate analysis. The goals of PCA are to (a) extract the most important information from the data, (b) compress the size of the data set by keeping only important information, (c) simplify the description of the data set, and (d) analyze the structure of the observations and the variables [20]. This is achieved by coordinate system translation, which means original variables are transformed to a new set of variables, i.e., the PCs, which are uncorrelated and reordered so that the first few retain most of the variation present in all the original variables.

Computation of the PCs reduces to the solution of an eigenvalue-eigenvector problem to a positive symmetric matrix. In this study, the unknown fluorophore distribution is considered as the variable. To apply PCA to the spatial and spectral MEFT images, fluorophore distribution images at each wavelength are assumed to be samples of variable, and the intensity of each pixel in the images is the observation of the sample. The variables can be denoted as $X=\{X_1,X_2,\mathrm{...},X_s,\mathrm{...},X_M\}$, where $M$ is the number of spectrum used and $X_S$is column vector which is the reconstructed 3D image frame of a given excitation light spectrum $s$. Each image frame has $N=R\times Q\times H$ pixels, where $R$ is the number of rows of each image slice, $Q$ is the number of columns of each image slice and $H$is the number of image slices. Then data of one image frame can be represented as $X_s=(x_{1,s},x_{1,s},\mathrm{...},x_{k,s},\mathrm{...},x_{N,s})$, where $x_{k,s}$ is the density of the pixel $k$ in the image frame at the given spectrum $s$, and the size of the whole 3D image frame $X$ is a $N\times M$ matrix. Then, the covariance or correlation matrix, $L$, is defined as follows,

The PCs are orthogonal meaning that each component is uncorrelated to each other and the variance of $P_j$ is equal to the eigenvalue corresponding to $E_j$. While transforming $X_0$ to $P$, the total variance keeps the same (equals to the sum of eigenvalues), which means that the PCA does not change the total energy. All PCs can then be arranged according to the variance. Each PC is further utilized to create an image, from which the main components of the samples may be obtained. In this study, in order to prevent the situation that variables with high variances dominate the PCs, correlation matrix instead of covariance matrix is adopted.

To sum up, the flowchart of the proposed algorithm is shown in Fig. 1.

 

Fig. 1 The flowchart of the proposed algorithm.

Download Full Size | PDF

3. Method

3.1. Simulation studies: cylinder model

3.1.1. Setup

The algorithm was first tested on a cylinder model set on a rotating stage. The cylinder model was set on a rotating stage. The rotational axis of the model was defined as the Z-axis with the bottom plane set as Z = 0 cm. The fluorophores were excited by a point source at the height of Z = 1.0 cm. Fluorescence images of 360° full view were collected with 12 projections at every 30° at each spectrum [34].

The model was composed of a cylinder containing two spheres, as shown in Fig. 2. The radius of the cylinder was 1.5 cm and the height was 2.0 cm. Two spheres s1 and s2 with the same radius of 0.2 cm were placed at the same height of Z = 1.0 cm and have an edge-to-edge distance of 0.2 cm. The cylinder was supposed to be filled with 1% intralipid with homogenous optical properties. The concentrations of two spheres were set as that in heart and lungs in mice [28, 29]. The concentration in the sphere s2 (lungs) was set as 32.5$\mu \text{M}$ (controlled concentration), and the concentration in the sphere s1 (heart) was set to range from 40% (13$\mu \text{M}$) to 200% (65$\mu \text{M}$) of the controlled concentration, with an increment of 20% (6.5 $\mu \text{M}$). There were totally nine simulation pairs of concentration combination. Here $\mu \text{M}$stands for $\mu$mol/L.

 

Fig. 2 The cylinder model used in the simulations. A cylinder of 3.0 cm diameter is filled with 1% intralipid. The spheres s1 (magenta) and s2 (green) contained in the cylinder stand for the heart and lungs, respectively, and are assumed to have evenly distributed fluorophores. The two spheres are placed at the same height with an edge-to-edge distance of 0.2 cm.

Download Full Size | PDF

3.1.2. MEFT model

For each concentration pair, the excitation light of different wavelengths was used to illuminate the media. Two excitation wavelengths, 740 nm and 780 nm, and one fluorescence light wavelength of 820 nm were employed. When the wavelength was changed, the optical parameters of the background media varied accordingly.

The procedure of the simulations was described as follows. Firstly, a series of fluorescent yields (concentrations of fluorophores $C$, extinction coefficients $\epsilon$ [33] and quantum yield $\eta$ [35]) was set to the spheres s1 and s2 at two wavelengths, as shown in Fig. 3. Secondly, the optical parameters of 1% intralipid were set at two wavelengths according to [36]. Finally, the tomographic images of fluorophore distribution were reconstructed at each spectrum for each concentration pair. And then these tomographic images were assembled to form the MEFT image frames.

 

Fig. 3 The fluorescence yield of ICG with different concentrations under different excitation light. The horizontal axis represents the concentration of fluorophores in the unit of $\mu M$and the vertical axis represents the corresponding fluorescence yield of the excitation light wavelengths 740 nm (blue) and 780 nm (green).

Download Full Size | PDF

3.1.3. Reconstruction of synthetic data

A zero-mean Gaussian noise was added to the fluorescence data. The signal-to-noise ratio (SNR) was 100 (20 dB), chosen based on the performance of currently available systems [37]. The excitation and emission data were synthesized using finite element method. The geometry in Fig. 2 was discretized into 3131 nodes and 15561 tetrahedral elements in finite element method software MultiPhysics 3.5a (COMSOL, Burlington, MA, USA). The detectors were located on the boundary nodes between 0.6 cm and 1.4 cm height range and inside 140° FOV corresponding to each point source. A volume of 3.0 × 3.0 × 1.0 cm³ (between 0.5 cm and 1.5 cm in height) was considered and sampled to 31 × 31 × 11 voxels for reconstruction. The conventional algebraic reconstruction technique (ART) was used to solve the derived linear equation and was terminated after 50 iterations [38].

3.1.4. Principal component analysis for synthetic MEFT images

There were nine concentration pairs in total. For each concentration pair, the input data were two reconstructed image frames. Each frame with 31 × 31 × 11 voxels was reshaped as a column vector containing 10,571 elements. By the transformation, the set of input data was placed in a matrix $X$ with 10,571 rows and 2 columns. Before PCA, $X$ was normalized to $X_0$by subtracting each column by the column mean and then dividing each column by its standard deviation. PCA was applied to $X_0$ and two PCs were obtained and arranged by the variance in descending order. The positive components and negative components were displayed in different colors. Each positive and negative element of PC illustrates structures with different concentrations [15, 27].

3.2. Simulation studies: digital mouse model

3.2.1. Setup

In order to further test the algorithm, numerical simulations were also performed on an irregular model, which employed the data from a 3D digital mouse model [39]. As shown in Fig. 4, The mouse torso containing heart and lungs was chosen as the target region, with a total length of 2.5 cm. Similarly to the cylinder model, the rotational axis of the mouse was defined as the Z axis with the bottom plane set as Z = 0 cm. The mouse was rotated over 360° with 30° increments, and 12 excitation point sources placed at the height of Z = 1.5 cm were used for illumination. In order to validate the accuracy of the reconstructed results and to evaluate the performance of PCA in resolving structures of organs, corresponding XCT data supplying the structure information were downloaded from the website (http://neuroimage.usc.edu/Digimouse.html).

 

Fig. 4 Geometry of the mouse torso region including the main target organs, i.e., the heart (magenta) and lungs (green).

Download Full Size | PDF

3.2.2. MEFT model

In the digital mouse model, the heart and lungs were set as target organs, while other parts in the mouse torso were set as background tissue. To simulate photons propagation, a heterogeneous mouse model was set up. The optical parameters were assigned to corresponding organs according to [40]. To simplify the problem, optical properties outside the organs were regarded as fat. The ICG concentration in the heart was set to 30 $\mu \text{M}$while in the lungs was set to 50 $\mu \text{M}$. As a non-target fluorescence dye, ICG would be retained in background tissue to some extent. So the ICG concentration in the background tissue which was taken as background fluorescence was set to 16 $\mu \text{M}$. In order to construct the MEFT model, the excitation light of two wavelengths 740 nm and 780 nm was used to illuminate the digital mouse model. Similarly to the cylinder model, series of fluorescent yields were set to heart, lungs and background tissue at the two wavelengths.

3.2.3. Reconstruction of synthetic data

The digital mouse model was discretized into 8,685 nodes and 44,060 tetrahedral elements. A volume of 2.1 × 3.4 × 2.5 cm³ (between 0 cm and 2.5 cm in height) was considered for reconstruction and sampled to 22 × 35 × 26 voxels with a size of 0.1 × 0.1 × 0.1 cm³. The detectors were located on the boundary nodes between 0 cm to 2.5 cm height range and inside 140° FOV corresponding to each point source. Reconstructions were terminated after 50 ART iterations.

3.2.4. Principal component analysis on synthetic MEFT images

The PCA procedure was similar to that in section 3.1.4. Briefly, the input data were two reconstructed image frames. Each frame with 22 × 35 × 26 voxels was reshaped as a column vector containing 20,020 elements. The matrix $X$ of input data was assembled with 20,020 rows and 2 columns. Then $X$ was centered and standardized to $X_0$. After performing PCA on $X_0$, the PCs were obtained. The positive component and negative components were displayed in different colors.

3.3. Phantom experiments

3.3.1. Setup

Phantom experiments were performed to further test the performance of the proposed algorithm. Similar to the simulated cylinder model, two transparent glass tubes (0.4 cm in outer diameter and 0.6 cm in length) filled with different concentrations of ICG were immerged in a cylinder phantom. The phantom was made of a glass cylinder (3.0 cm in diameter and 4.6 cm in height) filled with 1% intralipid. The two tubes had an edge-to-edge distance of 0.2 cm.

In the experiments, the phantom was placed on a rotating stage that allowed rotation and shift of the target around Z axis for collecting projections evenly distributed over 360°. The small light spot from a 250W Halogen lamp traveled through a set of band-pass filters and focused on the surface of imaged object. The set of excitation band-pass filters had the passbands of 740 ± 6 nm and 780 ± 6 nm, respectively. A 512 × 512 element, a −70°C cooled charge-coupled device (CCD) camera (iXon DU-897, Andor Technologies, Belfast, Northern Ireland) was placed on the opposite side of the imaged object, collecting the photons propagating through the imaged object. When collecting the fluorescence images, an 820 ± 6 nm band-pass filter was placed in front of the CCD camera. When collecting the excitation light images, a neutral density filter of 1% transmittance was used. A white light bulb was used to replace the excitation light and collect the white light images for recovering the 3D geometry of the imaged object.

3.3.2. MEFT model

Tube 1 depicted the ICG distribution in the heart with a concentration of 30 $\mu \text{M}$, while tube 2 depicted the ICG distribution in the lungs with a concentration of 50 $\mu \text{M}$. For each spectrum, 12 excitation and emission images were collected at every 30° using EMCCD with the excitation source at the height of 2.8 cm. The reduced scattering coefficients and absorption coefficients of 1% intralipid were set according to [36]. The total power of the point source was about 20 mW. 72 white light images acquired in 5° steps were back-projected to form a 3D geometry of phantom [31].

3.3.3. Reconstruction of experimental data

The cylinder phantom was discretized into 17,486 nodes and 88,286 tetrahedral elements. The detectors were located on the boundary nodes between 3.2 cm height range and inside 140° FOV corresponding to each point source. A volume of 3.1 × 3.1 × 2.2 cm³ (between1.8 cm and 4.0 cm) was considered for reconstruction and sampled to 32 × 32 × 23 voxels. Reconstructions were terminated after 50 ART iterations.

3.3.4. Principal component analysis for experimental MEFT images

The PCA procedure was the same as in section 3.1.4. Briefly, the input data were two reconstructed image frames. Each frame with 32 × 32 × 23 voxels was reshaped as a column vector containing 23,552 elements. By transformation, the set of input data was placed in a matrix $X$ with 23,552 rows and 2 columns. $X$ was then centered and normalized to $X_0$. After performing PCA on $X_0$, the positive components and negative components of PCs were displayed in different colors.

4. Results

4.1. Simulation studies: cylinder model

4.1.1. Reconstruction of simulation data

There are totally nine concentration pairs after combination. For each concentration pair, two image frames are obtained, with excitation light wavelengths of 740 nm and 780 nm, respectively. The reconstructed fluorescence slices (Z = 1.0 cm) at both excitation light of all nine concentration pairs are shown in Fig. 5. The two small red circles indicate the true locations of the fluorophores. As shown in the slices, the two fluorophores are difficult to distinguish.

 

Fig. 5 The reconstructed fluorescence slices (Z = 1.0 cm) at the excitation light wavelengths of 740 nm (upper row) and 780 nm (lower row) of all concentration pairs. From left to right, the according concentration pairs are shown below the lower row. All the images are displayed in the same range.

Download Full Size | PDF

The 3D reconstruction results with an excitation wavelength of 780 nm and 32.5-52$\mu \text{M}$ concentration pair are shown in Fig. 6. Figure 6(a) is volume rendering and Fig. 6(b) is surface rendering using an isosurface value equal to 20% of the maximum value of the volume data.

 

Fig. 6 The 3D visualization results of reconstructed fluorescence frame of 32.5-52 $\mu M$concentration pair at the excitation light wavelength of 780 nm. (a) Volume rendering. (b) Surface rendering.

Download Full Size | PDF

The results in Figs. 5 and 6 show that the two tubes almost emerged together and it is quite difficult to resolve the locations and the structures of each tube.

4.1.2. PCA on synthetic MEFT images

After the fluorescence image frames $X_0$ is assembled, the correlation (r) of the nine concentration pair samples is analyzed, as shown in Table 1. High correlation (r>0.994) between each variable is observed, indicating the abundant redundant information in the image frames.

Table 1. The correlation coefficients of the nine concentration pairs

View Table | View all tables in this article

PCA are performed on $X_0$, then a two-by-two coefficients matrix is obtained. Each column of the coefficients matrix contains coefficients for one PC which is the unitized and orthogonal eigenvectors. Each eigenvector consists of two values. The columns are arranged in the order of decreasing component variance. After analysis of the eigenvalues and eigenvectors of nine concentration pairs, we find that the eigenvectors almost have the same values in the precision of 10⁻⁴, although the eigenvalues are different. The corresponding eigenvectors$E$of the eigenvalues (latent roots ${\lambda }^{\text{'}}$) are shown in Table 2. The first row contains values of the eigenvalues ${\lambda }^{\text{'}}$ arranged in a decreasing order, and the corresponding eigenvectors are shown as columns below the eigenvalues. From the eigenvectors, the contributions of each variable to the PC are represented. The eigenvector $E_1$ illuminates that two variables in $X_0$ contribute almost the same to the first PC (PC1) which means that PC1 is approximately the mean value of two variables in $X_0$. The eigenvector $E_2$ illuminates that the second PC (PC2) is the result of variable 2 minus variable 1 in $X_0$, which illustrates the difference of two variables in $X_0$. We mainly focus on PC2 for PC2 illustrate the difference of fluorescence yields between spectrum.

Table 2. The eigenvalues ${\lambda }^{\text{'}}$ (latent roots) and the corresponding eigenvectors

View Table | View all tables in this article

PC2 contains positive and negative values, which are processed separately. The positive and the negative components of PC2 at the height of Z = 1.0 cm are displayed in Fig. 7. The upper row is the positive component depicted in magenta color while the lower row is the negative component depicted in green color. As shown in Fig. 7, the spheres s1 and s2 appear separately in the positive or negative components for every concentration pair except 32.5-32.5 $\mu \text{M}$, which means that the PCs reflect different structures of fluorophores of different concentrations.

 

Fig. 7 The PC2 images obtained from PCA applied to MEFT frames. The upper row is the positive components depicted in magenta color and the lower row is the negative components depicted in green color. The two small red circles show the true location and boundary of the spheres. The corresponding concentration pairs are shown below the lower row. All the images are displayed in the same range.

Download Full Size | PDF

The positive and negative components of PC2 image of 32.5-52 $\mu \text{M}$are also displayed in the same image in Fig. 8(a). As shown in Fig. 8(a), the two spheres are separated. To demonstrate the performance of PCA at all slices, 3D visualization results of the negative and the positive PC2 images for 32.5-52$\mu \text{M}$ concentration pair are obtained by extracting isosurfaces from all obtained PC2 images(with an isosurface value equal to 20% of the maximum value), as shown in Fig. 8(b).

 

Fig. 8 The obtained PC2 image frame of 32.5-52 $\mu \text{M}$ concentration pair. (a) The image slice at Z = 1.0 cm. (b) The 3D surface rendering results (by extracting isosurfaces). The green sphere is sphere s1 and the magenta sphere is sphere s2.

Download Full Size | PDF

The comparison between 3D reconstructed results at a specific wavelength of excitation light in Fig. 6(b) and the 3D results of PC2 in Fig. 8(b) demonstrates that organs and functional structures with different kinetic patterns can be separated when PCA is applied to MEFT image frames.

4.2. Simulation studies: digital mouse model

4.2.1. Reconstructions of simulation data

The reconstruction results of the digital mouse model from the same slice (Z = 1.5 cm) and different excitation wavelengths are shown in Fig. 9. Using the structure information of XCT data of digital mouse, the boundaries of heart and lungs are drawn. It is difficult to separate the structures of heart and lungs containing different concentrations of ICG based on each single spectrum tomographic image.

 

Fig. 9 The reconstructed fluorescence slice images (Z = 1.5 cm)of the digital mouse model at two excitation wavelengths (a) 740 nm and (b) 780 nm. . The boundary of mouse surface is depicted in red, the boundary of heart is depicted in yellow, and the boundary of lungs is depicted in mauve. All the images are displayed in the same range.

Download Full Size | PDF

4.2.2. Principal component analysis for synthetic MEFT images

Similarly to the cylinder model, before PCA is applied to the fluorescence image frames $X_0$, the correlation of the two samples is analyzed. Abundant redundant information is shown between two spectrums for the correlation coefficient is 0.9773. Similar to the cylinder simulation, the second principal component PC2 is mainly focused.

As shown in Fig. 10, PC2 images are displayed in the form of cross-sectional image. The height of image slice is Z = 1.5 cm. The structure information of organs by XCT is drawn in line with different colors in the image slice. The distributions of ICG in the heart and the lungs are respectively illustrated using the positive and the negative PC2 images in Figs. 10(a) and 10(b). The 2D merged result of positive and negative components in PC2 image is displayed in Fig. 10(c). The 3D merged result in Fig. 10(d) is the surface-rendering result of the PC2 images (by extracting isosurfaces of 20% of the maximum value). The results in Fig. 10 suggest that function structures of fluorophores of different concentrations can be resolved utilizing PCA in the digital mouse model.

 

Fig. 10 The PC2 images. (a)The positive component in PC2 is labeled in magenta color and (b) the negative component in PC2is labeled in green color. (c) The 2D merged results of the PC2 image. (d) The 3D merged results of the PC2 images. The boundary of mouse surface is depicted in red, the boundary of heart is depicted in yellow, and the boundary of lungs is depicted in mauve.

Download Full Size | PDF

4.3. Phantom experiments

4.3.1. Reconstructions of phantom data

The reconstruction results of the phantom from the same slice and different excitation wavelengths are shown in Fig. 11. The two reconstructed images at the same height (Z = 3.0 cm) are shown. Similar to in section 4.1, it is difficult to separate the structures of two tubes containing fluorophores of different concentrations based on each single spectrum tomographic image.

 

Fig. 11 The reconstructed fluorescence slice (Z = 3.0 cm) images at two excitation wavelengths (a) 740 nm and (b) 780 nm. The two small red circles show the true location and boundary of two spheres. All the images are displayed in the same range.

Download Full Size | PDF

Figure 12 shows 3D reconstructed results of experimental data at the wavelength of 780 nm. Figure 12(a) is volume rendering and Fig. 12(b) is surface rendering using an isosurface value equal to 20% of the maximum value of the volume data.

 

Fig. 12 The 3D visualization results of reconstructed fluorescence frame at the excitation light wavelength of 780 nm. (a) Volume rendering. (b) Surface rendering.

Download Full Size | PDF

4.3.2. Principal component analysis for phantom MEFT images

Similar to the simulation model, before PCA is applied to the fluorescence image frames $X_0$, the correlation of the two samples is analyzed. The correlation coefficient is 0.9985, which indicates the abundant redundant information.

After performing PCA to $X_0$, a two-by-two coefficients matrix is obtained. Each column of the coefficients matrix is the unitized and orthogonal eigenvector. Similar to the simulation study, the second principal component PC2 is mainly focused.

As expected, when PCA is applied to experimental MEFT images the fluorophores in two tubes in PC2 can be visualized and separated. As shown in Fig. 13, the fluorophores in tube 1 and tube 2 are illustrated using the negative and the positive PC2 images, respectively (upper row). The height of image slice is Z = 3.0 cm. The two small red circles show the true location and boundary of two tubes. The merged results of the positive and the negative P2 images shows the fluorophores in two tubes further (lower row). The 3D merged result is the surface-rendering result of the PC2 images (by extracting isosurfaces of 20% of the maximum value).

 

Fig. 13 The PC2 images.(a)The positive component in PC2 is labeled in magenta color and (b) the negative component in PC2 is labeled in green color. (c) The 2-D merged results of the PC2 image. (d) The 3D merged results of the PC2 images. The magenta parts indicate the fluorophores in tube 1 and the green parts indicate the fluorophores in tube 2.

Download Full Size | PDF

When PCA is applied to MEFT images, the two tubes can be separated, which means that the structures of two fluorehpores of different concentrations can be obtained.

5. Discussion and conclusion

Fluorescence tomography plays an important role in drug research for allowing non-invasive, quantitative, 3D imaging of biological and biochemical processes in living subjects. However, for the nature of the absorption and scattering to light of the media, it is difficult to get a high resolution of media structure. Especially in the in-vivo case, each functional unit (such as organ) has its own structure in which the ratio of HbO₂ to Hb is different and hence the optical parameters are different. So it is difficult to resolve organs with different kinetics in small animals. In this paper, we proposed a PCA-based method to aid in the identification of functional structures, and then evaluated its performance using simulations and phantom experiments.

As shown in the reconstructed images at a given spectrum wavelength, the two tubes with an edge-to-edge distance of 0.2 cm cannot be separated in either simulations or phantom experiments (Figs. 5, 6, 11 and 12) and the adjacent organs (i.e., heart and lungs) cannot be resolved in the digital mouse simulation (Fig. 9). Meanwhile it is difficult to recognize the structure of the object. In MEFT, the situation can be improved, but is still a challenge. In contrast, when PCA is applied to the MEFT image frames, a set of PC images is generated, which are the linear combination of image frames. The main component can be shown in the main PC images. In this study, the components are fluorescence yields that illustrate functional structures with different kinetic patterns. In the simulations of a cylinder model, the spheres on behalf of heart and lungs containing fluorophores of different concentrations are illustrated in the negative and the positive PC2 images. In the simulations of a digital mouse model, the structures of adjacent organs (heart and lungs) are separated in the negative and the positive PC2 images, respectively. Similarly, spatial structures of fluorophores in two tubes are separated in the negative and the positive PC2 images in the phantom experiments, because the concentrations of ICG solutions in two tubes are different. The locations and structures of fluorophores with different concentration distributions can be obtained in both simulations and phantom experiments.

In this study, only two wavelengths were used. With more wavelengths, more accurate results may be obtained. The problem of to what extent of concentration difference the proposed method can work is complex. As shown in Fig. 7, two inclusions with the same concentration 32.5$\mu \text{M}$cannot be resolved in the cylinder model, for fluorescence yields of the two inclusions having the same spectrum variation trends. This is mainly because the data used to proceed is the fluorescence yield that is affected by extinction coefficient $\epsilon$ and the concentration $C$. Extinction coefficient $\epsilon$ varies with many factors and concentration difference is related to not only the ratio of concentration but also the absolute value. As shown in the simulation results, the proposed method performs well for all different concentration pairs.

According to the results, the PCA-based method provides an attractive approach in kinetic study of fluorophores in organs in vivo. Each organ in the body plays a different role in the course of fluorophores circulation, accumulation or metabolism, and so each organ shows different concentrations of fluorophores. Therefore, the fluorescence yield of each organ shows different variation trends with spectrum. By taking advantage of this trend, PCA can resolve different fluorophores in organs.

All the obtained results are based on a basic assumption that accurate optical properties are known. When the method is applied to small animals, tissue optical properties are important for FT reconstruction and wrong optical parameters may lead to significant errors to the reconstructed results [41]. It may directly affect the performance of PCA identifying regions with different structures. Moreover, the uptakes of ICG are considered homogenous within organs for simplification in a normal mouse [28]. While in some particular cases, such as organs with tumors, the uptake of ICG within organs will presumably be heterogeneous and the heterogeneity may affect the analysis results of the proposed approach to some extent. In these cases, the regions inside and outside tumors may be taken as different parts in order to obtain more reasonable result. Ultimately, after performing PCA on MEFT, the principal components only reflect the structure information, so there is no direct relationship between the values of PC2 in the images to the known concentrations. The structure of fluorophores of different concentrations rather than the accurate values of the fluorescence yields in the structures are obtained in this study.

We mainly concerned on PCA to MEFT in this paper. However, fluorescence yields of fluorophores with different concentrations presumably show different variation trends with the emission spectrum in the emission-resolved fluorescence tomography. So this work could work for other kinds of multispectral FT. It is also worth noting that ART was used as reconstruction algorithm for its convenience and popularity in this work. But classic ART takes much time to compute and converge. The convergence speed can be significantly improved by selecting proper projection access order in which data are accessed in ART [42]. Other reconstruction algorithms [43] or package, such as NIRFAST [44], may be referred to in our future work.

In conclusion, the capabilities of the PCA-based method in separating structures containing fluorophores with different concentrations are demonstrated in simulations and phantom experiments. The proposed method not only can improve the detection precision, but also may be used as a base for structure study and dynamic study. Future work will focus on applying the method in separating structures of organs in small animals in vivo.