Depth-recognizable time-domain fluorescence molecular tomography in reflective geometry

: Conventional fluorescence molecular tomography (FMT) reconstruction requires photons penetrating the whole object, which limits its applications to small animals. However, by utilizing reflective photons, fluorescence distribution near the surface could be reconstructed regardless of the object size, which may extend the applications of FMT to surgical navigation and so on. Therefore, time-domain reflective fluorescence molecular tomography (TD-rFMT) is proposed in this paper. The system excites and detects the emission light from the same angle within a field of view of 5 cm. Because the detected intensities of targets depend strongly on the depth, the reconstruction of targets in deep regions would be evidently affected. Therefore, a fluorescence yield reconstruction method with depth regularization and a weighted separation reconstruction strategy for lifetime are developed to enhance the performance for deep targets. Through simulations and phantom experiments, TD-rFMT is proved capable of reconstructing fluorescence distribution within a 2.5-cm depth with accurate reconstructed yield, lifetime, and target position(s).


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][2][3][4][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][7][8][9], applying Lp-norm regularization for sparser reconstruction image [10][11][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][17][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).

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. As shown in Fig. 1(b), the rectangular phantom with a size of 70×40×50 mm 3 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 1 and D 2 (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.

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: where C is a constant depending on the fluorescence molecule. Φ f (r,t) denotes the photon density of the excitation light. Φ e (r,t) denotes the Green's function of the emission light. X(r) denotes the distribution of fluorescence yield and L(r) denotes the distribution of the reciprocal of lifetime, i.e., inverse lifetime, of the fluorescence molecule. The * operator represents convolution of time dimension.
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: where L i denotes the inverse lifetime estimated from the ith TPSF. To ensure the accuracy of L i , only the TPSFs of neighboring source-detector pairs are taken into consideration. And because is the sensitivity function of the ith source-detector pair, L'(r) is the weighted average of the sensitivity functions integrated over time, whose weights are the corresponding L i . In addition, L i is estimated as the solution of the equation When T 2 = 2T 1 , the equation is quadratic, and L i could be easily solved.
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 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 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: where X denotes the fluorescence yield distribution to be reconstructed. Y S denotes the vector of measurements which are the derivatives of TPSF. W S is the coefficient matrix of the inverse problem. And the jth row of W S is given by | t=T j , which is the temporal derivative of the ith sensitivity function at time point T j . It is also named as the sensitivity map and is flattened into a row. For derivative-based reconstruction, a sensitivity map is the probability distribution of the trajectory changes of the photons from an excitation source to a detector at a particular time point T. On the other hand, if the reconstruction is based on the TPSF values, the sensitivity maps would denote the probability distribution of the trajectory of the photons from a source to a detector at a particular time point T. Derivative-based sensitivity maps would increase independence among measurements and enhance spatial resolution and stability of reconstruction.

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: min where X denotes the distribution of fluorescence yield to be reconstructed. W S denotes the coefficient matrix. Y S denotes the vector of measurements. α is the regularization coefficient. D λ is a diagonal matrix. The diagonal elements of D λ are the depth regularization parameters given by D λ , in which d is the depth of the element, d 0 is an offset depth and η is a scaling coefficient. In this paper, d 0 is set as 15 mm and η is set as 0.6 mm −1 . Then, D λ would be normalized by its maximum element.

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 −2 to 1×10 −3 . β 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.

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: where L r (r) denotes the inverse lifetime distribution to be reconstructed. Tar m denotes the region of the mth target. During the lifetime reconstruction of the mth target, only the elements belonging to Tar m are modifiable while the others are fixed as values of L'(r) mentioned above in this paper or values of L r (r) that have been solved. Y V denotes the vector of the TPSF values at specific time points. λ is the regularization parameter whose initial value is 1×10 −10 multiplied by the mean of Y V and would decay over iterations. µ is a smoothing parameter whose value is 1×10 2 . F(L r (r)) is a nonlinear function and its kth element is given by where t k , Φ k f and Φ k e correspond to the kth elements of Y V . C m denotes the weight matrix according to contribution of the mth target in each TPSF. C m is a diagonal matrix, with the kth diagonal element being 1, 10, 50, 2×10 2 or 1×10 3 , depending on whether the contribution rate of the mth target in the TPSF corresponding Tpsf k (t)dt is greater than 0, 0.5, 0.35, 0.6 or 0.9. Moreover, for the steepest descent method, the element of the gradient F(L r (r)) is given by ∇F k r n = X(r n )Φ k f (r n , t) · Φ k e (r n , t) · (1 − L r (r n ))e −L r (r n )t | t=t k .

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 −1 and the reduced scattering coefficient µ s ' is set as 1×10 3 mm −1 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 3 . 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 1 ), the depth of Target 2 (D 2 ), horizontal distance (HD), and the edge-to-edge distance (EED), are shown in Table 1. 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 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 Table 2 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 = ∑︁ r ∈Tar

X(r)L(r)
/︃ ∑︁ r ∈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. 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 homogeneoustarget 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 −1 , while the inverse lifetime of targets with 2 µM ICG/ACH is 1.45 ns −1 , 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.
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. 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

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.