Dynamic fluorescence molecular tomography (DFMT) is a rapidly developing technique that can provide three-dimensional (3-D) metabolic information for biological studies and drug development.^1–3 Images of the pharmacokinetic parameters (also known as parametric images) can reflect the absorption and excretion processes of the fluorescence agents in biological tissues. Conventionally, the parametric images are obtained by curve fitting a sequence of FMT images into a compartmental model.^1,3 Due to the illposedness and ill-condition of the FMT problems, the reconstruction quality of the parametric images is limited. Although many methods have been proposed to improve the quality of the reconstruction results,^2,4 challenges in distinguishing different metabolic processes of the fluorescent targets, especially when the distances between the targets are very small, still remain.

Owing to the advantage of obtaining more information in multiple spectrums, the multispectral imaging technology has been applied to bioluminescence tomography⁵ and FMT.⁶ It has been proved that the excitation-resolved multispectral imaging leads to higher accuracy of the reconstructed image than the emission-resolved multispectral detection.⁶ Therefore, the excitation-resolved multispectral imaging is applied in this study.

In this letter, an excitation-resolved multispectral DFMT (MDFMT) method, based on the fact that the fluorescent targets with different concentrations show different variations in the excitation spectral domain, is proposed. In our previous studies, the spatial locations and structures of fluorescent targets with different concentrations in multispectrum FMT images are effectively unmixed by using multivariate image analysis methods, such as principal component analysis⁷ and independent component analysis (ICA).⁸ In this study, we further recover the fluorescent yield of each independent component and investigate the individual metabolic process of each fluorescent target. The results of simulations and phantom experiments show that the proposed method can effectively improve the reconstruction accuracy of the parametric images.

To investigate the metabolic processes of the fluorescent targets, an excitation-resolved MDFMT system⁹ is used to obtain the fluorescence datasets. During the MDFMT experiments, the fluorescent measurements are acquired at $K$ different time points, and the concentrations of the fluorescent targets vary with time. At each time point, the imaged object is sequentially illuminated using an excitation light with $S$ different wavelengths. At each wavelength, the imaged object is rotated for 360 deg in about 1 min with an angular increment of 15 deg, and 24 projections of the fluorescent measurements are acquired with a charge-coupled device camera. Therefore, a total of $P=K\times S$ groups of the fluorescent measurements are collected for the entire MDFMT experiments.

In the FMT forward problem, the photon transportation in biological tissues can be described by the diffusion approximation. After discretizing the imaged object into $N$ voxels, the forward model can be described with a linear system:⁹

where ${\mathbf{W}}_s(s=\mathrm{1,2},\dots ,S)$ is the wavelength-dependent weight matrix, which is generated according to the optical parameters at different wavelengths. ${\mathbf{u}}_s^k(k=\mathrm{1,2},\dots ,K)$ is the vector of the fluorescent measurements for the $s$’th excitation wavelength at the $k$’th time point, and ${\mathbf{x}}_s^k$ is the vector of the corresponding fluorescent yield (also known as FMT image). The fluorescent yields at different time points and excitation wavelengths can be obtained with traditional FMT reconstruction methods, and the Tikhonov regularization is selected in this study.¹⁰

Taking into account multiple excitation spectrums, the fluorescent yields can be defined as follows:^7,8

where $\eta$ stands for the quantum yield, which is often assumed to be constant for a given fluorescent agent and environment, and is set to be 0.12 in this study.⁸ $C$ and $ϵ$ are the concentration and the extinction coefficient of the fluorescent target, respectively. In this study, indocyanine green (ICG) is used as the fluorescent agent. According to previous studies,¹¹ the extinction coefficients of ICG are affected by both concentrations and wavelengths because it tends to aggregate in water at high concentrations. Therefore, in excitation-resolved MDFMT, $ϵ$ varies with ${\lambda }_e$, which means that for different $C$, $x$ shows different variation trends along with ${\lambda }_e$. The fluorescent targets with different concentrations can be considered independent signal sources, which can be decomposed by multivariate image analysis methods.^7,8

After FMT reconstruction, a total of $K$ groups of reconstruction results are obtained (corresponding to $K$ different time points), and each group contains $S$ FMT images (corresponding to $S$ excitation wavelengths). To resolve the different components inside the imaged object, the FMT images of each group are decomposed by ICA, which has the general form:¹²

where $S$, $N$, and $P$ are the number of excitation spectra, the number of discrete voxels, and the number of ICA decomposed components, respectively. $\mathbf{X}={({\mathbf{x}}_1,\dots {\mathbf{x}}_i,\dots {\mathbf{x}}_S)}^T$ is the observation matrix in which ${\mathbf{x}}_i$ is the $N$-dimensional FMT image obtained at the $i$’th excitation wavelength. $\mathbf{I}={({\mathbf{IC}}_1,\dots ,{\mathbf{IC}}_j,\dots ,{\mathbf{IC}}_P)}^T$ is the independent component ($\mathbf{IC}$) matrix, and ${\mathbf{IC}}_j$ is an $N$-dimensional vector that contains the spatial structural information of the $j$’th component. $\mathbf{E}=({\mathbf{SC}}_1,\dots ,{\mathbf{SC}}_j,\dots ,{\mathbf{SC}}_P)$ is the spectrum course ($\mathbf{SC}$) matrix. The $j$’th column ${\mathbf{SC}}_j$ is an $S$-dimensional vector that reflects the absorption ability of ${\mathbf{IC}}_j$ to the excitation light. Generally, the dimensions of $\mathbf{X}$ and $\mathbf{I}$ are assumed to be equal (i.e., $P=S$). Therefore, in this study, Eq. (3) at a given time point $k$ can be expanded as follows:⁸

Assuming that the target appears in the $d$’th IC, the true fluorescent yield of this target can be recovered by setting the other ICs to be zero and multiplying by SCs as follows:

where ${\widehat{\mathbf{x}}}_s^k$ represents the recovered fluorescent yield for the $s$’th excitation wavelength and at the $k$’th time point; $\mathbf{Z}$ denotes the $N$-dimensional zero vector.

Finally, the pharmacokinetic parameters of each fluorescent target are obtained by fitting the fluorescent yields into a two compartmental model. Let ${\phi }_s^j={[A_s^j,B_s^j,{\alpha }_s^j,{\beta }_s^j]}^T$ denote the pharmacokinetic parameter vector of the $j$’th voxel at the $s$’th spectrum and $t_k$ represent the $k$’th time point; then, the corresponding parameters can be calculated as follows:^1,3

where $A$ and $B$ are the zero-time intercepts, which denote the initial concentration of the fluorescent agents. $\alpha$ and $\beta ({\min }^{-1})$ are the metabolic rates, which reflect the absorption and excretion processes of the fluorescent agents in organs and biological tissues.

The actual fluorescent targets can be determined based on two criteria. First, the objects should repeatedly appear in the ICs at different time points. Second, the metabolic properties of the objects should have physiological significance, which reflects the absorption and excretion processes of the fluorescent agents in biological tissues. The objects that do not satisfy these criteria should be considered noise components.

Numerical simulations and phantom experiments are carried out to estimate the performance of the proposed MDFMT method. To investigate the influences of the excitation wavelengths, several wavelength combinations with different numbers of spectra are tested in both the simulations and phantom experiments. The results demonstrate that the different metabolic properties can be effectively separated as long as the number of spectra $S$ is larger than the number of fluorescent targets. Moreover, when $S$ increases, the accuracy of the parametric images improves while the temporal resolution decreases. Therefore, $S$ should be properly selected to balance the accuracy and temporal resolution of the parametric images. In this study, four wavelengths, 700, 720, 740, and 760 nm, are selected in the excitation, and the fluorescent wavelength is set to be 840 nm in both the simulations and phantom experiments.

A 3-D digital mouse atlas is employed in the numerical simulations. As depicted in Fig. 1(b), the mouse torso with a total height of 3 cm containing the heart (red) and lungs (green) is selected as the target region. The region outside the organs is regarded as adipose. As listed in Table 1,¹³ different absorption coefficients ${\mu }_{\mathrm{a}}$ and reduced scattering coefficients ${\mu }_{\mathrm{s}}^{\prime }$ are assigned to corresponding regions to construct a heterogeneous model. To mimic the metabolic processes of ICG in the two organs, 50 groups ($K=50$) of different fluorescent concentrations varying with time are set according to the metabolic curves shown in Fig. 1(a). The curves are obtained with Eq. (6), and the pharmacokinetic parameters are listed in Table 2.^14,15 The fluorescent yields ${\mathbf{x}}_s^k$ are calculated based on Eq. (2) in which the extinction coefficients $ϵ$ at different wavelengths and concentrations are set according to the parameters reported by Landsman et al.¹¹ The corresponding boundary measurements ${\mathbf{u}}_s^k$ are calculated according to Eq. (1). A zero-mean white Gaussian noise with a signal-to-noise ratio of 40 dB is added to ${\mathbf{u}}_s^k$ to generate noisy measurements.

Table 1

Optical parameters at different wavelengths in the simulations and phantom experiments.

Table 2

Setups of the pharmacokinetic parameters for different regions in the simulations and phantom experiments.

The 3-D FMT reconstruction results at the time point of $t_k=10\min$ and the excitation wavelength of 700 nm are shown in Fig. 1(c). It is very difficult to distinguish the locations and the structures of the two organs. Four ICs are obtained with the decomposition of ICA, and the results are depicted in Figs. 1(e)–1(h). IC1 and IC2 are noise components, while IC3 [Fig. 1(g)] and IC4 [Fig. 1(h)] correspond to the lungs and heart, respectively. The two organs are effectively extracted from the FMT images, and the metabolic processes of these two organs can be separately analyzed. By combining IC3 and IC4, the merged 3-D result is shown in Fig. 1(d), which is very similar to the true setups, as depicted in Fig. 1(b).

Figure 2 shows the cross-sectional parametric images obtained with the proposed MDFMT and the traditional DFMT methods. The images correspond to the slice indicated by the blue dashed line in Fig. 1(d). The MDFMT results are obtained at the excitation wavelength of 760 nm, and similar results can be obtained at other excitation wavelengths. As depicted in Figs. 2(d1)–2(d4), the two organs can be clearly distinguished in the parametric images obtained with the MDFMT method. By comparison, as shown in Figs. 2(e1)–2(e4), it is difficult to recognize the different metabolic processes in the DFMT results.

The normalized root mean square error ($\mathrm{NRMSE}={\Vert \widehat{e}-e\Vert }_2/{\Vert e\Vert }_2$, where $e$ and $\widehat{e}$ stand for the true and reconstructed parametric images, respectively) is used to quantify the relative errors between the true and the reconstructed images. A smaller NRMSE indicates higher accuracy of the reconstructed image.⁴ As listed in Table 3, the NRMSEs of the four parametric images obtained with the MDFMT method are smaller than those with the DFMT methods, which quantitatively confirms that the proposed MDFMT method obtains higher reconstruction accuracy than the traditional DFMT method.

Table 3

NRMSEs of the parametric images obtained with the two methods in the simulations and phantom experiments.

The phantom experiments are carried out on a cylinder with a height of 4 cm and a radius of 1.5 cm, as depicted in Fig. 3(b). The imaged phantom is filled with 1% intralipid, and the optical parameters at different wavelengths are listed in Table 1.^8,16 Two tubes with an edge-to-edge distance (EED) of 0.1 cm are placed inside the cylinder to be the fluorescent targets. The height and radius of the tubes are 0.5 and 0.2 cm, respectively. Ten groups of ICG with different concentrations are filled in the tubes to mimic the metabolic processes of two organs. The concentrations are set according to the metabolic curves at 10 different time points ($K=10$), as indicated by the red stars and green circles in Fig 3(a). The pharmacokinetic parameters of the two curves are listed in Table 2.

The 3-D reconstruction results are shown in Fig. 3. The four ICs decomposed with ICA are depicted in Figs. 3(e)–3(h) in which IC3 [Fig. 3(g)] and IC1 [Fig. 3(e)] correspond to tubes 1 and 2, respectively, and the other two components are noises. Figure 3(d) shows the merged results of IC1 and IC3. The structures and spatial locations of the two tubes in the merged results are very similar to the true setups [Fig. 3(b)]. By contrast, as shown in Fig. 3(c), the two tubes cannot be distinguished in the FMT reconstruction results.

Figure 4 shows the cross-sectional parametric images obtained with the MDFMT and DFMT methods, and the slice is indicated by the red solid line in Fig. 3(d). The MDFMT results are obtained at the excitation wavelength of 760 nm, and similar results can be obtained at other excitation wavelengths. Due to the small EED of the fluorescent targets (0.1 cm), the two tubes are mixed together and cannot be distinguished in the parametric images obtained with the traditional DFMT method, as depicted in Figs. 4(e1)–4(e4). By comparison, as shown in Figs. 4(d1)–4(d4), the different metabolic processes of the two tubes can be clearly observed in the reconstruction results obtained with the MDFMT method. It can be seen from the results that the MDFMT method obtains higher image quality than the DFMT method. The findings are quantitatively proved by the NRMSEs listed in Table 3 in which the MDFMT method obtains smaller NRMSEs for all the parametric images than the DFMT method.

In conclusion, an excitation-resolved MDFMT method is proposed to improve the accuracy of the parametric images. Different fluorescent targets in the imaged object are decomposed using ICA, and the metabolic processes of the targets are separately analyzed. As depicted in Figs. 2 and 4, the metabolic characteristics of different fluorescent targets can be clearly observed in the parametric images obtained with the MDFMT method but cannot be distinguished in the DFMT results. Moreover, as listed in Table 3, the proposed MDFMT method obtains smaller NRMSEs in both the simulations and phantom experiments, which quantitatively confirms that the MDFMT method can achieve higher accuracy of the parametric image than the traditional DFMT method. It is worth mentioning that the fluorescent yields, not the fluorescent concentrations, are fitted to the two compartmental models to obtain the pharmacokinetic parameters because the extinction coefficient $ϵ$ in Eq. (2) is unknown. This may affect the accuracy of the parametric images and will be further studied in our future work.