1. Introduction

Fluorescence molecular tomography (FMT) is a functional imaging modality which is capable of noninvasively detecting in-vivo biological activities on molecular level. It provides quantitatively accurate three-dimensional (3D) distributions of endogenous or exogenous fluorescent biomarkers in thick and highly scattering media such as biological tissues. Therefore, FMT is a vital tool in preclinical oncology research, drug development and other biological research areas [1–5]. Despite its advantages of non-ionizing radiation and low cost, the emerging clinical applications of FMT have been hampered not only by its low spatial resolution caused by the ill-posed nature of its reconstruction, i.e., inverse problem, but also by the limited penetration of near infrared photons in biological tissues. The ill-posed nature of the reconstruction stems from the strong scattering of near infrared photons in biological tissues. Conventionally, the spatial resolution of FMT could be improved by several methods including utilizing early photons that suffer less from scattering [6–9], applying Lp-norm regularization for sparser reconstruction image [10–12] and so on. Nowadays, FMT in both continuous-wave and time-domain modes could separate small targets with an edge-to-edge distance of about 1 mm within a 30-mm-diameter field of view [9,11], which has already met the needs of most clinical applications. The penetration of near infrared in biological tissues is on the level of several centimeters, which would be the main reason that the applications of FMT are limited to small animals in preclinical studies. Conventional FMT requires photons penetrating the objects for adequate spatial information of the fluorescence distribution. Therefore, the section diameters of what FMT could be applied to is typically smaller than 5 cm, which tallies with the size of small animals such as nude mice or rats for preclinical studies. Few applications of FMT in human subjects aim at fingers or other parts with small section diameters [13].

To overcome the penetration limit of FMT, the most direct approach is to introduce detectors with much lower dark count rates to increase the detection sensitivity, which is unfortunately unavailable at present. On the other hand, in the aspect of diffuse optical tomography (DOT), it has been reported that by utilizing photons with reflective trajectory, there could be applications in parts of human body, such as necks and brains of infants, which near infrared photons cannot penetrate completely [14,15]. For instance, based on sources and detectors with spacings of 20 mm and 30 mm respectively, attached on the transverse plane, the reflected light was detected, and an optical anatomical neck model could be reconstructed for the feasibility study of diffuse optical imaging of malignant lesions in the thyroid gland [14]. Then it is rational to hypothesize that the penetration limit of FMT could be circumvented by a reflective geometry. In this case, reflected photons would provide the information needed to reconstruct the fluorescence distribution in the external region of the subject under scanning, regardless of the section diameter of the subject. Even though the reconstruction of fluorescence distribution in deep regions is still challenging, the distribution in the external region could provide abundant biological information about tissues or organs in this region such as skin, ligament, muscle, or organ exposed during surgery, which is useful in many clinical applications such as surgical navigation [16–18], diagnosis of sports injuries [19] and so on.

There are several studies using FMT in incomplete reflective geometry, which could be applied in objects with larger section diameters, even though the penetration limit would still exit. In a previous study [20], for example, the sources and detectors of an FMT system were placed with an incomplete reflective angle of around 90° and single fluorescent target located 5 mm below the surface of nude mice was successfully reconstructed. Reconstruction of multiple fluorescent targets remains challenging for FMT in reflective geometry, since deeper targets yield much weaker fluorescence, and the reconstruction would be significantly affected by shallower targets. To settle this problem, depth regularization which assigns depth-dependent weights to Lp-norm regularization would be applied to FMT reconstruction, inspired by previous work in DOT [21].

In this paper, a depth-recognizable time-domain fluorescence molecular tomography in reflective geometry (TD-rFMT) is proposed. In the TD-rFMT system, spatial information including depth information is provided not only by the different spacings and positions of the source-detector pairs, but also by the different propagation trajectories of near infrared photons that vary over time. Time derivatives of the data curves collected, i.e., temporal point spread functions, are chosen as reconstruction measurements, which have been proved to improve the quality of reconstruction [9]. Furthermore, with depth regularization applied to the reconstruction, it is proved in both simulations and phantom experiments that TD-rFMT is capable of reconstructing the distribution of multiple targets at diverse depths (i.e., different depths within a large range).

2. Setup and method

2.1 System setup

A fiber-coupled time-domain fluorescence molecular tomography system in reflective geometry is built by integrating a TCSPC module (SPC-150, Becker & Hickl GmbH, Germany), a PMT (PML-16-C, Becker & Hickl GmbH, Germany) and a femto-second laser generator (Spectra-Physics, Newport Corporation, Canada) which is set to work at 780-nm wavelength (80 MHz, 100 fs pulse-width), as shown in Fig. 1(a). The experimental system also consists of a rectangular phantom, a 16×2 optical switch for fibers, and a group of filters, including an achromatic doublet (AC254-030-B, Thorlabs, Newton, NJ) and two bandpass fluorescence filters with center wavelength of 840 nm (ff01-840/12-25, Semrock, Rochester, NY; XBPA840, Asahi Spectra, Torrance, CA) to eliminate the excitation light. The temporal resolution of the TCSPC module is about 12.5 ps, and the time window is 12.5 ns to ensure that the majority of the temporal point spread function (TPSF) could be obtained.

 

Fig. 1. Schematic diagram of the TD-rFMT system and phantom experimental setup. (a) TD-rFMT system. (b) Transverse view of the rectangular phantom. (c) The locations and FOV of excitation sources and detectors.

Download Full Size | PDF

As shown in Fig. 1(b), the rectangular phantom with a size of 70×40×50 mm³ is filled with 1% intralipid and two cylindrical targets with a diameter of 5 mm are positioned inside the phantom, while the horizontal distance (HD) between the centers of the targets and the depths of the targets D₁ and D₂ (i.e., the distances from the surface to the centers of the targets) vary according to different experimental setups. As shown in Fig. 1(c), on the surface of the phantom, 11 fibers are placed uniformly with an adjacent distance of 5 mm in a field of view (FOV) of 50 mm for both detection and excitation. These fibers are on the plane with the same height, which is named the excitation plane. During measurements, for each projection, the excitation laser is guided into the phantom through one of the fibers while the other 10 fibers serve as detectors. Thus, there are 10 probe points in each projection. Because the excitation laser is switched among all the fibers producing 11 projections, a total of 110 TPSFs are acquired, from which the distribution of the targets would be reconstructed.

2.2 Inverse problem

For the inverse problem of TD-rFMT, the forward model was calculated by telegrapher equation based on the finite element method (FEM) [7]. And so were all the simulations and phantom experiments in this paper. The relation between the distribution of fluorescence and the TPSF could be expressed as the following physical equation:

To reconstruct the distribution of fluorescence, it requires the distribution of inverse lifetime which is usually unknown. In this paper, we propose a tail-fit method to obtain an estimation of the inverse lifetime distribution L’(r). It would be proved in simulations and phantom experiments in this paper that with L’(r), reconstruction of the distribution of fluorescence yield X(r) is achievable, while the location and value of the reconstructed targets would not be affected. The tail-fit method is formulated as follow:

To illustrate the estimation of L_i, assumptions and approximations are made. Assuming that for the ith TPSF, the distribution of inverse lifetime is uniform with an effective value of L_i, the TPSF could be expressed as $T{\rm{p}}sf(t) = C \cdot A(t) \ast {L_i}{e^{ - {L_i}t}}$, in which $A(t) = \int\limits_r {X(r)\Phi _f^i({r,t} )\ast \Phi _e^i({r,t} )dr}$. In addition, assuming that the tail of the TPSF is dominated by the logarithm tail of A(t) and ${L_i}{e^{ - {L_i}t}}$, the TPSF could be further simplified as $T{\rm{p}}sf(t) = C{e^{ - Kt}} \ast L{e^{ - {L_i}t}}$, where K is the logarithm tail of A(t). So $\ln ({Tpsf(t)} ){\rm{ = }}\ln ({{e^{ - Kt}} - {e^{ - {L_i}t}}} )+ \ln \left( {\frac{{{L_i}}}{{{L_i} - K}}} \right) + C^{\prime}$. Since X(r) is unknown, K could be approximated by the logarithm tail of A’(t) which is K’ mentioned above. Then the equation could be obtained.

With the estimated lifetime distribution, the inverse problem of TD-rFMT could be formulated based on the first-order temporal derivatives of the TPSF at specified time points, while choosing derivative as data was proved to enable better reconstruction in previous work [9]. The inverse problem is shown as follows:

2.3 Reconstruction method with depth regularization

In previous work of DOT in reflective geometry [21], reconstruction sensitivity in deep regions is extremely low, resulting in false depth of deep targets in the reconstruction results when there are targets in shallow regions simultaneously. The same defect exists in reflective FMT. Depth-dependent regularization assigns smaller regularization parameters with larger depths and could improve the reconstruction uniformity in DOT [21]. Such a method is expected to improve the reconstruction performance of TD-rFMT as well.

In this paper, a depth-regularized Tikhonov-regularization-based projecting sparsity pursuit method (PrSP-Tk-D) was proposed as the reconstruction method with depth regularization to solve the inverse problem of TD-rFMT. The Tikhonov-regularization-based projecting sparsity pursuit method proposed in our previous work [22] has been proved to be a fast and stable reconstruction method. It is not sensitive to the regularization coefficient and is capable of achieving high spatial resolution and high image quality. The PrSP-Tk-D method is described as follows:

Algorithm 1. Depth-regularized Tikhonov-regularization-based projecting sparsity pursuit method

The process of the PrSP-Tk-D method is shown in Algorithm 1. α is the regularization coefficient and varies from 1×10⁻² to 1×10⁻³. β_ac is the acceleration coefficient and is set as 0.6, controlling the step size of acceleration projection [22]. γ_ac is the slackness factor typically set as 0.2. And ɛ_sparse ranging from 0.5 to 2, is the sparsity error factor controlling the sparsity of the solution. Algorithm 1 is very similar to the algorithm previously proposed in [22], with the main difference being the depth regularization parameters as mentioned above.

2.4 Weighted separation reconstruction of lifetime

With the fluorescence yield distribution reconstructed by the PrSP-Tk-D method, the distribution of the reciprocal of fluorescence lifetime could be directly reconstructed based on a time-domain strategy [23]. However, similar to fluorescence yield reconstruction, lifetime reconstruction in the deep regions has a low sensitivity. The reconstructed value of inverse lifetime in the deep regions would be much lower than expected if any target exists in the shallow regions. A weighted separation reconstruction strategy was proposed for accurate inverse lifetime reconstruction in the deep regions. Inverse lifetime reconstruction of targets from the reconstructed yield distribution would be carried out separately from shallow to deep regions. Moreover, for each target, the weight of each TPSF would be resigned according to the contribution of the target in the TPSF and the absolute value of fluorescence yield would be adjusted according to the resigned weights. Therefore, based on a previously reported L1 regularization projected steepest descent algorithm (L1PSD) [23], an L1 regularization weighted separation reconstruction algorithm (L1WSR) was proposed. The L1WSR algorithm is expressed as solving the following optimization problem:

3. Simulations and experiments

Simulations and phantom experiments in this paper were all based on the TD-rFMT system mentioned above. There are 10 detectors for each projection while there are 11 projections in total. In the calculation of simulations and forward model, the absorption coefficient μ_a is set as 3 mm⁻¹ and the reduced scattering coefficient μ_s’ is set as 1×10³ mm⁻¹ while the sampling temporal resolution is set as 25 ps, twice as the temporal resolution of the PMT. The forward model is discretized into mesh grid with resolution of 0.1×0.1×0.2 mm³. For the simulations, Gaussian noise with specific amount of energy was added to approximate the signal-to-noise ratio (SNR) in the phantom experiments. The SNRs in the simulations typically ranged from 10 dB to 60 dB for most of the TPSFs, although they could be less than 0 dB for some of the TPSFs.

To evaluate the performance of the proposed PrSP-Tk-D and L1WSR methods, 9 sets of simulations were carried out, including 5 sets of homogeneous-target simulations and 4 sets of heterogeneous-target simulations. In the simulations, there are two targets with a diameter of 5 mm and a height of 10 mm placed symmetrically about the excitation plane, whose values of fluorescence yield are the same. For homogeneous-target simulations marked as simulations 1–5, the lifetimes of the targets are both set as 1 ns. Meanwhile, for heterogeneous-target simulations marked as simulations 6–9, the lifetime of Target 1 is set as 1 ns while the lifetime of Target 2 is set as 0.6 ns. The position setups of simulations 1–9, including the depth of Target 1 (D₁), the depth of Target 2 (D₂), horizontal distance (HD), and the edge-to-edge distance (EED), are shown in Table 1.

Table 1. Position setups of simulations.

View Table | View all tables in this article

The reconstruction was performed using the PrSP-Tk-D and L1WSR methods. As shown in Fig. 2, the fluorescence yield distributions for all the simulations are successfully reconstructed, regardless of targets at different depths in Fig. 2(d) and 2(h), and even the deep targets with HD = 0 mm are recovered as displayed in Fig. 2(e) and 2(i). And it can be found that the spatial resolution decreases as the depth increases. As can be seen from Table 2, the positioning errors (PE) of the targets in the heterogeneous-target simulations 6–9 are mostly smaller than 1 mm, indicating that the tail-fit estimation method for L’(r) is practical in achieving accurate fluorescence yield distribution reconstruction even the lifetimes of targets are diverse (i.e., with a large range). The PE is given by PE = x_recon4 − x_true where x_true and ${x_{recon}} = {{\sum\limits_{r \in Tar} {X(r)x(r)} } \mathord{\left/ {\vphantom {{\sum\limits_{r \in Tar} {X(r)x(r)} } {\sum\limits_{r \in Tar} {X(r)} }}} \right.} {\sum\limits_{r \in Tar} {X(r)} }}$ are the horizontal coordinate or depth of the center of the true target and the value-weighted center (i.e., centroid) of the reconstructed target. In addition, it can be found from Table 2 that the PEs of targets in the homogeneous-target simulations are also mostly smaller than 1 mm while the PEs in depth (DE) of Target 2 are slightly better than those in the heterogeneous-target simulations. It is because larger inverse lifetime would result in larger deviation by the tail-fit method. For heterogeneous Target 2 with inverse lifetime of 1.67 ns⁻¹, the deviation could be about 0.4 ns⁻¹ while the deviation for Target 1 is smaller than 0.1 ns⁻¹. Meanwhile, it can be seen that although depth regularization improves the performance, diverse depths still result in slightly lower performance in terms of positioning and relative reconstructed values (RV). The RV is given by RV = EV₁:EV₂ in which effective value (EV) is defined as the average yield value of the corresponding target.

 

Fig. 2. Reconstructed yield images on the excitation plane (a)-(i) for simulations 1–9. All the images are normalized by their maximum values.

Download Full Size | PDF

Table 2. Positioning errors (PE), relative values (RV), effective inverse lifetimes (EIL) and the relative errors of EIL (EILE) in the simulations.

View Table | View all tables in this article

The reconstructed inverse lifetime images of simulations 1, 3, 4 and 6–8 are shown in Fig. 3. As can be observed in Fig. 3(a)-(c), the inverse lifetime distributions of the homogeneous targets are similar, while the inverse lifetimes of the heterogeneous targets are apparently distinguishable in Fig. 3(d)-(f). It proves that the L1WSR method in TD-rFMT is capable of qualitative inverse lifetime reconstruction. Furthermore, as shown in Table 2, the effective inverse lifetimes (EIL) of all the simulations are close to the true inverse lifetimes, with relative errors of EIL (EILE) smaller than 10% except for the targets with true lifetime of 0.6 ns in the deep region. The EIL is calculated as the value-weighted average of inverse lifetime $EIL = {{\sum\limits_{r \in Tar} {X(r)L(r)} } \mathord{\left/ {\vphantom {{\sum\limits_{r \in Tar} {X(r)L(r)} } {\sum\limits_{r \in Tar} {X(r)} }}} \right.} {\sum\limits_{r \in Tar} {X(r)} }}.$ As can be seen in simulation 5, because the yield distribution is not commendably reconstructed, the inverse lifetime of the target in the deep region has been affected, resulting in larger errors of EIL. And although the weighted separation reconstruction strategy ensures lifetime reconstruction of the homogeneous targets in the deep region, there are still improvements to be made for more accurate reconstruction of quantitative distribution of targets with large inverse lifetime. Therefore, it is proved in the simulations that accurate depth-recognizable fluorescence yield and lifetime images of TD-rFMT could be achieved by the PrSP-Tk-D and L1WSR methods.

 

Fig. 3. Reconstructed images of inverse lifetime on the excitation plane (a)-(c) for homogeneous-target simulations 1, 3, 4 and (d)-(f) for heterogeneous-target simulations 6–8. Effective lifetimes (EL, reciprocal of EIL) of the targets are shown. Zero values were removed.

Download Full Size | PDF

The performance of TD-rFMT and the methods is further evaluated in the phantom experiments. Hence, 8 sets of phantom experiments were carried out, which included 4 sets of homogeneous-target experiments and 4 sets of heterogeneous-target experiments. In the experiments, two tubes with a diameter of 5 mm were buried inside the rectangular phantom and placed symmetrically about the excitation plane. In the homogeneous-target experiments numbered as 1–4, the tubes filled with 10 μM Indocyanine Green/dimethylsulphoxide (ICG/DMSO) with liquid height of 20 mm served as homogeneous targets. On the other hand, in the heterogeneous-target experiments numbered as 5–8, one tube filled with 2 μM ICG/alcohol (ACH) was marked as Target 1 and the other filled with 10 μM ICG/DMSO was marked as Target 2. According to our measurements, the inverse lifetime of targets with 10 μM ICG/DMSO is 0.95 ns⁻¹, while the inverse lifetime of targets with 2 μM ICG/ACH is 1.45 ns⁻¹, and the fluorescence yield ratio of target with 2 μM ICG/ACH to target with 10 μM ICG/DMSO is 0.54:1. The position setups of phantom experiments 1–8 are shown in Table 3. Besides the lifetime and yield of Target 1, the position of targets in heterogeneous-target experiments 5–8 were the same as those in homogeneous-target experiments 1–4, respectively.

Table 3. Position setups of phantom experiments.

View Table | View all tables in this article

As shown in Fig. 4, the reconstruction results of the phantom experiments agree with the results of the simulations. As shown in Table 4, the targets are reconstructed successfully with most HD errors smaller than 1 mm and positioning errors in depth (DE) smaller than 1 mm. In addition, for heterogeneous targets with diverse yield values, the relative values are close to the true ratio 0.54:1, which further proves the capability of the PrSP-Tk-D method in accurate yield reconstruction. However, in the experiments with diverse depths, albeit with depth regularization, larger PEs still occur, especially for targets in the deep region. And as can be seen in Fig. 4(e)-(d), the depths of targets with larger inverse lifetime appear to be shallower than expected, compared with the simulation results. And for deep targets close to the shallow targets in the horizontal direction, as shown in Fig. 4(d) and 4(h), the reconstructed yield is not accurate enough in terms of either shape or PE. It may be on account of the shadowing effects. In the deep region right below the shallow targets, the reconstruction sensitivity would be much lower than that of other deep regions. This is because in any measured TPSF, the intensity contribution of the shallow targets would be dominant compared with that of the targets in the deep region below it.

 

Fig. 4. Reconstructed yield images on the excitation plane (a)-(h) for phantom experiments 1–8. All the images were normalized by their maximum values.

Download Full Size | PDF

Table 4. Positioning errors in depth (DE), horizontal distances (HD), relative values (RV), effective inverse lifetimes (EIL) and the relative errors of EIL (EILE) in the phantom experiments.^a

View Table | View all tables in this article

The reconstruction results of inverse lifetime in the phantom experiments also agree with the simulation results. As shown in Fig. 5, the similar inverse lifetime distribution of homogeneous target as and the apparently distinguishable distribution of heterogeneous targets support the good performance of the L1WSR method in TD-rFMT in the phantom experiments. And the relative errors of effective inverse lifetime in the phantom experiments are mostly smaller than 10%, which further validates the accuracy of lifetime reconstructed by the method, although the shadowing effects mentioned above have harmed the results of phantom experiments 4 and 8. Although the performance is imperfect and the spatial resolution of fluorescence yield needs further improvement, the PrSP-Tk-D and L1WSR methods in TD-rFMT show good performance in the phantom experiments, with accurate yield values, accurate lifetimes, and the depth problem are largely alleviated.

 

Fig. 5. Reconstructed images of inverse lifetime on the excitation plane (a)-(d) for homogeneous-target phantom experiments 1–4 and (e)-(h) for heterogeneous-target phantom experiments 5–8. Effective lifetimes (EL, reciprocal of EIL) of the targets are shown. Zero values were removed.

Download Full Size | PDF

4. Discussion

In conclusion, depth-recognizable time-domain fluorescence molecular tomography in reflective geometry was proposed in this work. TD-rFMT not only could circumvent the penetration limit of FMT, enabling its application in objects with large section diameters, but also could prevent potential excitation light leak during measurements that happens in conventional geometry that requires transmission. In addition, TD-rFMT could be non-contact and does not require the rotation of the object, which will furtherly benefit its applications. The depth-recognizable reconstruction of fluorescence yield distribution and inverse lifetime distribution were achieved by the depth-regularized Tikhonov-regularization-based projecting sparsity pursuit method and the L1 regularization weighted separation reconstruction method proposed along with TD-rFMT in this paper. The reconstruction of fluorescence yield distribution is based on derivative of the TPSF and an estimated lifetime distribution, which provide the information for accurate reconstruction. And it is proved in the simulations and phantom experiments that by using the PrSP-Tk-D method in TD-rFMT, fluorescence yield images could be successfully reconstructed with precise positioning, accurate relative value, and acceptable spatial resolution, no matter when the targets are located in the shallow region, in the deep region or in both regions. Furthermore, with the reconstructed yield distribution, inverse lifetime distribution of each target could be separately reconstructed by the L1WSR method. However, L1WSR method cost a large amount of computation. Computation cost could be reduced by approximating the nonlinear problem with a linearized method [24]. The simulations and phantom experiments also proved that the reconstructed inverse lifetime distributions are accurate. An accurate reconstructed lifetime distribution is valuable because it is absolute quantitative, and lifetime of the fluorescence molecule could be quantitative indicator of coefficients of the biological environment such as pH and temperature [25].

However, the spatial resolution of fluorescence yield images and accuracy of inverse lifetime images are not as good as expected for extreme cases like phantom experiments 4 and 8. In the extreme cases, the target in the deep region was placed right below the shallow target and as mentioned above, the contribution of the deep target in each measurement would be much lower than that of the shallow one. Therefore, further improvements are required. A deep learning method that could enhance the performance of the extreme cases and the spatial resolution in the deep regions would be taken into consideration in the future.