1. Introduction

The maturation of near infrared (NIR) fluorescence diffuse optical tomography (DOT) [1–5] and the recent advances in red-shifted and NIR molecule-specific fluorescent probes [6–9] have considerably accelerated the development of fluorescence molecular tomography (FMT). As a novel optical imaging modality of visualizing and quantifying the function of genes and expression of enzymes and proteins deep inside intact living organisms (mainly small animals) [10–14], the potential applications of FMT lie in many fields of biomedical research: 1) non-invasively and near-simultaneously monitoring multiple and molecular events in a specific biochemical pathway in-vivo; 2) optimizing drug and gene therapy and evaluating their effects at a molecular and/or cellular level; 3) assessing disease progression at a molecular pathological level, and 4) rapidly, reproducibly and quantitatively observing time-dependent experimental, developmental, environmental, and therapeutic influences on gene production in the same subject [15, 16]. The image reconstruction approaches of FMT/ fluorescence DOT in steady-state and frequency-domain modes have been fairly demonstrated with simulated and experimental data, which are mostly derived from the Born approximation based on either the relatively simple photon diffusion model [1–5, 12–14] or the more physically rigorous radiative transfer equation [17], it is natural and imminent to extend the techniques to time-domain regime [12, 13], where not only the simultaneous recovery of fluorescent yield and lifetime but also the analysis of multiple components could be implemented in a direct way. The two-parameter reconstruction provides information concerning not only the localization but also the local environment (oxygen, [Ca²⁺] and pH etc.) of fluorophores, while the multi-component analysis, where multiply sorts of molecular targets are traced by using different imaging reporter probes (components), potentially enables assessment of the multi-gene controlling mechanism in a disease progression as well as simultaneous observation of the multiple molecular/cellular events in a specific biochemical pathway [15, 16]. As with time-domain DOT, time-domain FMT can also be performed in either the full time-resolved scheme [18–20] or the featured-data one [18, 21, 22], with the former considerably improving the quantification and spatial resolution of the reconstruction at large computation cost as well as unappealingly high noise-sensitivity [19], while the latter being fast and robust but with degraded image quality due to the loss of the information and the correlation among the date-types [19, 21, 22].

The featured-data algorithms in time-domain DOT that are based on either Mellin or Laplace transform of the time-dependent diffusion model have been proposed by many researchers [18, 21, 22, 25], and their effectiveness as well as robustness to both the systematic errors and measuring noises successfully validated with numerical, in vitro and in vivo experiments [23–26]. Among these, the one that is based on the real-domain Laplace-transformed diffusion model, referred to as the generalized pulse spectrum technique (GPST), has been shown to possess many advantages over those that employs the Mellin-transformd forward model in flexibility, stability, robustness and computational efficiency [25, 27]. The extension of the methodology to time-domain FMT is apparently straight.

In the following sections, we firstly develop a linear GPST scheme of time-domain FMT for simultaneously recovering the fluorescent yield and lifetime distributions in the one- and two-component cases, based on the finite-element method (FEM) solution to the Lapalce-transformed coupled diffusion equations. Then, the proposed algorithms are validated by simulated data for the numerical phantoms that are designed to fairly demonstrate the faith and efficacy of the methodology. Finally, the crucial issues related to the practice of the algorithm are discussed and conclusions presented.

2. Methodology

Analogous to DOT, the task of image reconstruction in FMT is also expressed as an inverse issue for a given photon-migration model. We describe in this section the FEM frameworks for both the forward and inverse models of time-domain FMT, aiming at simultaneous reconstruction of fluorescent yields and lifetimes of multiple components from the time-resolved measurements.

2.1 Forward Model

In the field of diffuse light imaging, the diffusion equation, i.e. the P ₁ approximation to the radiative transfer equation, has been prevalently employed as the photon-migration model, which is mathematically tractable either analytically or numerically [18]. For time-domain FMT, where an ultra-short pulsed laser at the excitation wavelength is launched into tissue and the resulting time-dependent fluorescent emission on the boundary is measured by a time-resolved system, the forward model that governs the relationship between the excitation and emission light propagations in fluorescent turbid medium is thus given by the time-dependent coupled diffusion equations. To apply the GPST, the time-dependent forward model is converted into the complex frequency domain by using Laplace transform

where subscripts x and m denote the excitation and emission wavelengths, respectively; Φ_v(r,r _s,p) = ${\int }_0^{+\propto }$ Φ_v(r,r _s,t)e ^-pt dt (v∈[x,m]) is the Laplace transform of the time-dependent photon density Φ_v(r,r _s,t) with a complex transform-factor p ; the optical parameters involved are the absorption coefficient μ_av(r), the reduced scattering coefficient μ′_sv(r) and the diffusion coefficient D_v (r,t) = c/[3μ′_sv(r)]; the fluorescence parameters are the fluorescent yield ημ(r) and lifetime τ(r). These quantities, in general, are functions of the position vector r. We employ uniquely the Robin boundary condition for the above equations

where K = (1 + R_f )/(1 - R_f ) and R_f ≈ -1.4399n ^-2 + 0.7099n ^-1 + 0.6681 + 0.0636n is the internal reflection coefficient at the air-tissue boundary, with n being the relative refractive index of tissue to air [28]. The measurable flux, i.e., the data-type at the boundary site ξ_d (d = 1,2,…,D)) and for the source site ζ_s (s = 1,2,…,S), can be calculated by the Fick’s law with consideration of Eq.(2)

Eq. (1) can be efficiently solved with the Galerkin FEM by expanding Φ_v(r,r _s,p) into finite elements: Φ_v(r,r _s,p) ≈ ${\mathrm{\Sigma }}_{n=1}^N$ Φ_v(n,p) u_n (r) = Φ _v(p)^T u(r) with u(r) = [u ₁(r), u ₂(r), ⋯, u_N (r)^T and Φ _v(p) = [Φ_v(1,p),Φ_v(2,p),⋯,Φ_v(N,p)]^T being the shape functions and the Laplace-transformed photon density at the N nodes of the FEM mesh, respectively, resulting in following matrix equation

The elements of the coefficient matrices A _v, B and Q _v are given below

where i and j index the N nodes of the mesh; ημ_af(j) and τ(j) the fluorescent yield and lifetime at the j-th meshing node, respectively.

2.2 Inverse model

According to Eq. (1) and by applying the Fick’s Law at emission wavelength, we have an integral equation as follows

where Γ_m(ξ_d, ζ_s,p) is the Laplace transform of the transient emission flux measured at boundary site ξ_d and for excitation site ζ_s, and G_m (ξ_d, r, p) the flux at ξ_d for a source at r, predicted by the diffusion equation at the emission wavelength

Let x(r,p)≈ ${\mathrm{\Sigma }}_{n=1}^N$ x_n (p)u_n (r) = x ^T(p)u(r) , where x(p) = [x ₁(p), x ₂(p),…, x_N (p)]^T , Eq.(6) can be discretized into a matrix equation

where

and

with the element given by

where Ω_n numerates all the elements that are joined at the n-th node; Ḡ^(Ω_m) _m(ξ_d,p) and $\overline{\Phi }$ ^(Ω_m) _x(ξ_d,p) are the mean values of G_m (ξ_d,r,p) and Φ_x(r,ζ_s,p) at the nodes of the element Ω_n, respectively. By the sense of Eq. (8), the inverse procedure is based on the mesh nodes, which usually represents a significant reduction in the number of the unknowns, as compared with those that are based on the elements. This will be further explained in the next section.

For solving Eq. (8) that is in general of large-scale, underdetermined and ill-posed, a Kacmarcz method, commonly known as the Algebraic Reconstruction Technique (ART) might be very efficient [18, 22, 29]

where Γ ⁽ⁱ⁾(p) is the i-th element of Γ(p) and W ⁽ⁱ⁾(p) the i-th row of W(p), with (i mod j) denoting the modulus after i divided by j ; M is referred to as the ART-circulation number. The initial point x ₀(p) is set to what reflects the morphology of the background fluorescence as closely as possible. In principle, the ART sequentially projects a solution estimate onto the hyperplanes defined by the individual rows of the linear system. The relaxation parameter λ has a significant influence on the reconstruction and has been proven in a range of [0, 2] to make the algorithm converge to a point on the intersection of the governing equations that is nearest to the initial point [29]. The regularization strategy in the ART-based algorithms is accomplished by limiting the number of the iterations, whose choice is task-dependent and mandatory in presence of noise. The primary advantage of this implicit inversion method over the other schemes for solving underdetermined linear systems, such as the Levenberg-Marquardt, the truncated Singular Value Decomposition and the Conjugate-Gradient algorithms, is its near independence of memory occupation, since the rows of the weight matrix are successively employed during the solution process. However, the accuracy of the ART strongly relies on the initial guess of the unknowns. Sometimes, the prior information on the background as well as the targets is necessarily required to attain an image reconstruction of high-quality. To suppress the artifacts arising in the above ART inversion, a median filter operating on the adjacent nodes (i.e., those that belong to the same element as the output node) is employed at the end of each circulation of the ART inversion.

2.2.1 One-component case

In the normal one-component case, one kind of biochemical molecule of interest is targeted by its specific probe, and a set of two unknown distributions (one component), i.e. the fluorescent yield ημ_af(r) and lifetime τ(r) of the probe, is to be recovered. This can be explicitly done from the images of x(r,p ₁) and x(r,p ₂) by employing a pair of transform-factors: p ₁ and p ₂, in the Laplace transforms

2.2.2 Multi-component case

By the multi-component case we mean that multiple kinds of molecular targets under investigation are labeled by their respective probes of distinct fluorescent properties either in yield or lifetime (with an assumption that the probes have the same excitation and emission wavelengths as each others, and there are no interactions among them), and therefore multiple independent sets of the molecule-specific fluorescence distributions (multiple components), saying μ_afi(r) and τ_i (r), i = 1, 2,…, each related to one probe, need to be simultaneously recovered from the measured time-resolved data, which reflects the mixture of the yield-weighted decay characteristics of the multiple fluorescent probes. Consequently, the distribution x(r, p) in Eq. (6) changes to the following form

where N_c is the number of the components under investigation. The recovery of the 2N_c unknown distributions in Eq. (14) can in principle be achieved by using at least 2N_c transform-factors in the Laplace transforms: p ₁, p ₂,…, p2N_c . For the two-component case, this leads to a set of four joint equations as follows

where x(r,p_i ) (i = 1,2,…,2N_c ) is obtained from solving Eq. (8) for each of the transform-factors. Let

, Eq. (15) can be transformed into the following matrix equation

with

By solving Eq. (16) for α(r), we can further obtain the four unknown distributions of the fluorescence parameters, i.e., ημ_af1(r), τ₁(r), ημ_af2(r) and τ₂(r), in terms of Eq. (16)

For the case of N_c > 2, a similar procedure to the two-component case can be followed, which is mathematically more complex. It should be noted that within the above FEM framework Eq. (17) and Eq. (21) need to be solved successively for each node of the mesh, given the values of x(r) at the positions of the discrete nodes.

Although, according to Eq. (17), 2N_c , real transform-factors are theoretically enough for the reconstruction of the 2N_c , independent unknowns in a N_c -component case, this may bring about considerable errors in the images due to the evident underestimation of x(r, p) in Eq. (8). To overcome the adversity, more than 2N_c real transform-factors are in general employed, making Eq. (17) both overestimated and ill-posed. We therefore introduce an extended Levenberg-Marquardt method for the reliable solution, which is expressed as an optimization problem with the Tikhonov-Miller regularization as the follows [30]

where α _B(r) is the value of α(r) calculated for the known background fluorescent properties at each node. The regularization matrix R(r)is chosen in such a way that R ^T R (RR ^T) and M ^TM (MM ^T) commute each other, meaning that the two matrices have the same complete set of eigenvectors and the reversely-ordered eigenvalues, i.e., $\mathbf{M}=\sqrt{{\mu }_k}{\mathbf{u}}_k{\mathbf{v}}_k^T\mathrm{and}\mathbf{R}={\sum }_{k=1}^{R_M}\sqrt{{\mu }_{R_M-k}}{\mathbf{\mu }}_k{\mathbf{v}}_k^T,$ where R_M is the rank of M, u _k and v _k are the eigenvectors of MM ^T and M ^T M belonging to the nonzero eigenvalue μ_k. The solution to the above optimization problem can be readily found and leads to the matrix equation below

The fact that the quantities in the above equation are all the functions of position vector r implies that they needs to be solved for each FEM node, as indicated before. Mathematically, Eq. (22) or Eq. (23) seeks the filtered least-squares solution to Eq. (17) on the basis of its smoothness around the baseline (the known background properties). The smoothing weight is controlled by the regularizing factor β, whose value is task-dependent. Again, the median filtering is applied to α(r) images as well as the resultant ημ^af(r) and τ(r) ones, respectively, for a further suppression of the artifacts.

 

Fig. 1. FEM mesh and optode configuration employed in the study.

Download Full Size | PDF

Table 1. Optical, fluorescent and geometrical parameters of Phantom 1

View Table | View all tables in this article

Table 2. Optical, fluorescent and geometrical parameters of Phantom 2

View Table | View all tables in this article

3. Validations

Although, for computational simplicity, the validation of the proposed algorithm is performed on a two-dimensional circular domain with a diameter of R = 25 mm using simulated data, the principle of the methodology is applicable to three-dimensional models of the realistic geometry and the results represent an ideal experimental situation where the effects of systematic error on the reconstruction are ignored. To apply FEM, the domain is divided into 3750 triangles that join at 1951 nodes. 16 coaxial source-detector optodes (S = D = 16), are assumed at equal spacing around the annulus, of which the 16 detectors collect the exiting photons in parallel as the 16 sources illuminate the surface successively. This leads to a total of 256 time-resolved measurements. In our algorithm only measurements from the 9 detectors opposite to the illuminating source, i.e. a set of 144 time-resolved data, are employed for the image reconstruction. It was found in our previous studies that such a configuration of data set can effectively reduce the difference in the order of the data magnitude and therefore significantly improve the image quality [19, 25–27]. Again, we emphasize herein that, the inverse model in the proposed algorithm is node-oriented where the fluorescence parameters to be recovered are given at the nodes, and in the Galerkin-FEM forward calculation those at each element are calculated as the average of those at all the nodes consisting of each element. By this strategy, both the ill-posedness and computational complexity of the inverse model are considerably lessened as the result of reduction in the number of the unknowns, e.g., from 3750 to 1951 in our case, as compared to a full element-based scheme. Figure 1 shows the FEM mesh of the domain and the deployment of the optodes used throughout the paper.

 

Fig. 2. (a) Original (top) and reconstructed (bottom) images of fluorescent yield and lifetime for Phantom 1, and (b) their profiles along X- and Y-axes.

Download Full Size | PDF

 

Fig. 3. (a) Original (top) and reconstructed (bottom) images of fluorescent yield and lifetime for Phantom 2, and (b) original (red) and reconstructed (green) profiles along X-axis.

Download Full Size | PDF

 

Fig. 4. (a) Schematic of the phantom used for the evaluating spatial resolution of the algorithm, (b) reconstructed images for CCS=13 mm (top), 15 mm (middle) and 17 mm (bottom), respectively, and (c) profiles of the reconstructed yield and lifetime images along the X-axis.

Download Full Size | PDF

3.1 One-component case

In this case, we uniquely employ a pair of real transform-factors: p ^1,2 = ∓0.1P, where P = 1/[1|μ${(}_{\mathit{\text{ax}}}^{B)}$ c + 1/${\mathrm{\mu }}_{\mathit{\text{am}}}^{\left(B\right)}$ c + τ^(B)], for the Laplace transforms of the time-domain forward and inverse models, and a circulation number M = 20 as well as a relaxation parameter λ = 0.5 for the ART solution to Eq. (8), the derived linear system. The validation is firstly enforced by two numerical phantoms, for which the noiseless sets of the data-type Γ(p) are generated from the forward model with the same mesh as that to be used in the inverse model, to evaluate the intrinsic performance of the algorithm. In Phantom 1, four fluorescent circular targets with the same radius but different fluorescent properties are included to investigate the ability of the algorithm to discern the difference in the fluorescent properties of the targets, while Phantom 2 embeds two disks with same fluorescence parameters but different radiuses in order to probe the performance of the algorithm to reconstruct the target size. The optical properties of the targets in the two phantoms are the same as those of the background, which are set to μ_ax,m = 0.035mm^-1 and μ′_sx,m = 1.0mm^-1 for both the excitation and emission wavelengths. These values are in the range of the optical properties for in vivo muscle [31]. All the optical and fluorescent properties for Phantom 1 and 2 are listed in Table 1 and 2, respectively. The original and reconstructed images and their profiles along the axes for both the phantoms are shown in Fig. 2 and 3, respectively, where all the differences in the yield, lifetime and size among the targets are reasonably disclosed and further improvements in the spatial resolution and quantitative accuracy are desired. It is worthy to note in Fig. 3 that the reconstruction of the large sized target is much more quantitatively faithful than that of the small-sized one.

Secondly, we investigate the spatial resolution of the algorithm using a phantom with the same domain and background optical/fluorescent properties as the above examples. Two fluorescent disks with the same radius of 3 mm and optical/fluorescent properties of μ_ax,m =0.035 mm^-1 , μ′_sx,m =1.0 mm^-1, ημ_af =0.005 mm^-1 and τ = 1000 ps , but a variable center-to-center separation (CCS), are placed along the X-axis, as shown in Fig. 4(a). Figure 4(b) presents the reconstructed images and Fig. 4(c) their profiles along the X-axis, for three different values of the CCS: 13 mm, 15mm and 17 mm. It is seen from the results that the two disks can still be resolved as the CCS is less than 13 mm (the valleys drops from the peaks by about 9% and 5% for the yield and lifetime, respectively), and the recovery of the fluorescent yield appears to have better resolution but poorer quantitativeness (only about 45% of the true value is reconstructed by the peaks) than the lifetime in terms of the profiles.

 

Fig. 5. Investigation on the noise robustness of the algorithm by imaging the same phantom as in Fig. 4, with the CCS of the two target disks equals to 17 mm, for a varying SNR of 35 db (Top), 40 dB (Middle) and 45 dB (Bottom). (a) Reconstructed yield and lifetime images, and (b) their profiles along the X-axis.

Download Full Size | PDF

Finally, we demonstrate the robustness of the algorithm to noise by using noisy data with different signal-to-noise ratio (SNR). This time we use the same phantom as described above (contrasts 5:1 and 10:1 for the fluorescent yield and lifetime, respectively), but fix the CCS of the two targets at 17 mm. Originally, the noise is embedded in the measured time-resolved flux and is composed of many types. But here we model the noise directly in each of the featured data-types, i.e., the Laplace transforms of the time-resolved data, as an additive Gaussian random variable with a standard deviation proportional to the data-type: σ(ξ_d, ζ_s, p) = Γ_m(ξ_d, ζ_s, p)10^-χ/20, where d = 1,2,…,D , s = 1,2,…,S , and χ is the SNR in decibels. Figure 5 illustrates the reconstructed images of the phantom a varying SNR: χ = 35 dB , 40 dB and 45 dB . The results reveal that the noise-robustness of the algorithm is moderate and the reconstruction of the fluorescent yield is much more insensitive to the noise but less accurate in the quantitativeness than that of the lifetime.

3.2 Two-component case

Simultaneous reconstruction of two sets of fluorescent yield and lifetime images in the two-component case is in principle a direct extension to the one-component case but involves more complex mathematics, as aforementioned. Here we only present preliminary results for imaging a simple two-component phantom. The phantom has the same domain and background/target optical properties as in Fig. 4 and embeds two circular targets with the same radius of 4 mm and a fixed CCS of 25 mm. The background and targets all contain two kinds of fluorescent probes, as listed in Table 3. A set of 11 real transform-factors from -0.25P to +0.25P with a step of 0.05P is used for the Laplace transforms. This makes Eq. (17) overdetermined and ill-posed, with the condition number of its normal equation being 10⁵ in the order of magnitude, for which the least-square solution can be approximately found with the Levenberg-Marquardt scheme. The reconstruction is performed with Gaussian noisy data of χ = 45 dB and a fixed regularizing factor β = 10^-3. Figure 6 shows the original and reconstructed images for the two-component phantom, as well as their profiles along X-axis. It is observed that the reconstruction of the fluorescent yield correctly reflects the difference between the two components but both are underestimated, while the resultant images of the lifetime are somewhat distorted with the first component underestimated and the second overestimated.

Table 3. Fluorescent parameters in the two-component phantom

View Table | View all tables in this article

4. Discussions and conclusions

In general sense, the Laplace transform of the time-resolved data in the above GPST-FMT algorithm can be implemented in the complex-domain. In our simulations only real transform-factors are employed. Compared to the methods using the Fourier transform, which is the imaginary-domain case of the Laplace transform and can be directly provided by a frequency-domain system, those that are based on the real-domain Laplace transform are less computationally intensive since no complex-number arithmetic is involved in the inversion procedure. Unlike the Fourier transform whose intensity and phase concern the information on the fluorescent yield and lifetime, respectively, the Laplace transform utilizes the information embedded in the full time-resolved data in an implicit way: the Laplace transform using a positive or negative transform-factor can be physically interpreted as an exponential time weighting that favors “early-time” or “late-time” features of the time-resolved curves that are jointly influenced by the fluorescent yield and the lifetime, and thus a pair of positive and/or negative frequencies would effectively distinguish between them. Originated from this idea, for the multi-component case with N_c targets, a proper choice of at least N_c transform-factor pairs would selectively pick up the signal weights that favor the different fluorescent emissions, and lead to a reliable separation among the fluorophores. A major disadvantages of the real-domain Lapalace transform is that it is now unclear if the information embedded in the time-resolved data can be extracted effectively, as with the Fourier transform, which mathematically complete along the imaginary domain.

 

Fig. 6. (a) Original and (b) reconstructed images of the two-component phantom with the fluorescence parameters listed in Table 3: The top row is for the first component and the bottom row for the second component in each sub-figure. The original (red) and reconstructed (green) profiles along X-axis are shown in (c) and (d) for Component 1 and Component 2 respectively.

Download Full Size | PDF

The foundation of the above linear scheme of time-domain FMT lies on the assumption that the influence of the fluorophore absorption on the optical properties at both the excitation and emission wavelength is negligible, implicating that the targets to be imaged are of either low-contrast or small-size so that the evolvement of their fluorescent parameters during the reconstruction will not result in an evident change in the weighting matrix of Eq. (8). For accurately recovering targets of large-size and high-contrast, or for precisely quantifying fluorophores of low-concentration that is usual case in the diagnosis and pathogenesis research of early tumor, the influence of the probe absorption on the weighting matrix should be necessarily taken into account, even though some authors have demonstrated that the FMT image reconstruction might be insensitive to the background heterogeneity to some extend [32]. This eventually leads to a nonlinear inverse issue that needs to be solved with an iterative scheme.

It is rather difficult to characterize the noise in the measurement since it can come from many sources. Generally, two types of noise sources are dominant in the measured signals. The first is due to thermal noise in the detector electronics. This type of noise is independent of the measured signal and well modeled as a sequence of independent identically distributed Gaussian random variables. The second one, known as shot noise, originates from the quantum fluctuations of the detected photons, and is typically modeled as Poisson distributed in the fluence rate. This model can be approximated with independent, but not identically distributed Gaussian random variables, with an expectation that fluence rate will need to be reasonably high to exceed the thermal noise floor of the detectors. Since the real Lapalce transform of the original time-resolved signal can be in broad sense regarded as a windowed intensity, we thus, for the most simplification, use the approximated shot noise model aforementioned. Because the noise model does not have a constant variance across the source-detector pairs, a whitening procedure might be incorporated for more robust reconstruction. This is accomplished by weighting Eq. (8) with the inverse covariance matirx

where C ^-1 = diag(σ(ξ^d,ζ_s,p)) is the inverse of the covariance matrix of the data-type. Eq. (24) requires, however, that the covariance matrix is either known as a prior knowledge or estimated as auxiliary parameters within an inversion procedure [5]. It should be noted that the assumption of the constant SNR throughout the detectors is somewhat inconsistent with the real scenario where the SNR varies proportionally with the square root of the fluence rate present at each detector [5]. Nevertheless, the shot-noise model with a constant SNR (most favorably set to the minimum) helps reasonably identify the lower limit of the noise-robustness of the algorithm since a time-resolved multi-channel system usually try to equalize the measured intensities through attenuator adjustment. According to the status of the time-resolved instrumentation, in particular a time-correlated single photon counting technique, a noise-robustness of SNR around 40 dB would be favored in practice.

The above examples of the image reconstruction are based on the absolute values of the Laplace-transformed time-varying flux. The evaluations show that the noise-robustness of this absolute imaging scheme is moderate. This indicates some potential difficulties in achieving a high-sensitive FMT for small animal, where SNR might be rather low due to the strong background fluorescence, the inevitable systematic noises as well as probably inaccurate estimation of the background optical and fluorescent properties. Fortunately, the differential imaging scheme is available in the most practice of FMT, where the measurements before and after intravascular injection of the fluorescent probes can be obtained and the difference between them is used to form the data-types that is to be used for the reconstruction, i.e.,

where ${\mathrm{\Gamma }}_m^{\left(B\right)}$(ξ_d, ζ_s, p) and ${\mathrm{\Gamma }}_m^{\left(A\right)}$(ξ_d, ζ_s, p) are the data-types before and after probe injection, respectively; ${\mathrm{\Gamma }}_m^{\left(M\right)}$(ξ_d, ζ_s, p) the model prediction in terms of the background parameters. Differential imaging has been proved to be able to cancel many systematic and measuring noises and greatly improve the image quality [23–27]. It is also important to note that with the differential imaging scheme the three-dimensional FMT task can be performed in a layered imaging mode where the source-detector pairs are deployed to obtain the image at each of the multiple planes, by use of a two-dimensional inverse model. The validity of this method has been justified in DOT and FMT [10, 26].

In the two-component simulation, the regularized Levenberg-Marquardt method, i.e., Eq. (23), has been promisingly applied for generating the reliable α -images from the overdetermined but ill-conditioned linear system that are derived from the x -images at a number of different tramsform-factors, with a constant regularizing-factor β. Furthermore, it has been shown that an L-curve criterion for determining β can better trade off between the data misfit ||x(r) - M(r)α(r)|| and the solution-to-baseline seminorm ||R(r)[α(r) - α_B(r)]||, but at a considerably increasing cost of computation. By comparison, we found that this regularized overestimation scheme can produce much higher quality of images than the simple matrix inversion using only four transform-factors. It is also pointed out here that, although we confine the fluorescent heterogeneities for the two components to the same regions for simplicity, it is unnecessary in the methodology for the fluorescence regions of the multiple probes to be the same as each other. The allowance for the maximal freedom of the probe distributions agrees with the fact that the multiple probes specifically bound to their respective biochemical molecules have little probability of the same fluorescent emission behavior, due to the intrinsic differences among the distributions of the molecule concentrations and among the efficacies of the fluorescent probes.

The Laplace transformed time-dependent diffusion model with a real transform-factor p, i.e., Eq. (1), can mathematically interpreted as the time-independent coupled diffusion equations with absorption coefficients of μ_av + p/c (v ∈ {x, m}). Therefore, together with considering the exponentially decaying response of the fluorescent emission, the choice of the transform-factor p should meet the condition that p ≥ max{- μ_av c, - 1/τ} so that the time-independent diffusion model is physically meaningful as well mathematically stable. Besides this criterion, there is no theoretical guide to the choice of the transform-factor pairs. Therefore the transform-factor pairs used in the above examples are somewhat empirical. Nevertheless, the above numerical validations of the methodology with these empirically chosen transform-factor pairs have achieved our goal with considerable success. The selection of the transform-factor pairs that can optimally characterize the original time-resolved signal will be a challenge in the future investigation.

With the objective of developing a general methodology for the image reconstruction of TD-FMT that can work with a broad range of fluorescent properties of the probes that might be encountered in applications, the fluorescence parameters in the above examples do not all exactly conforms to those of the common fluorescent dyes, such as Indocyanine green (ICG) (with the extinction coefficient of about 0.013 mm^-1μM^-1, the lifetime of about 0.56 ns and the quantum efficiency of 0.016 at the peak excitation wavelength of 778 nm) and Cy5.5 (with the extinction coefficient of about 0.019 mm^-1μM^-1, the lifetime of about 1 ns and the quantum efficiency of 0.23 at the peak excitation wavelength of 670 nm). Nevertheless, with regard to these two dyes, the fluorophore concentrations corresponding to the chosen fluorescent yields of the targets are all on the orders of 1μM and 100 nM, respectively, as preferably expected in vivo.

According to the above examples, the target size and contrast significantly influences the image quality. A formal investigation on these effects is shown in Fig. 7, where an increasing improvement in the quantitativeness is observed with the size increasing and/or the contrast decreasing. This observation conforms to that found in DOT. For the two-component case, because the lifetime values of the two components are reciprocal in Eq. (16), an underestimation of one lifetime value may result in an overestimation of another one, and vice versa. This adverse effect has been clearly seen in Fig. 5, where the lifetime reconstruction for the first component is underestimated while that for the second slightly overestimated. In this sense, an effective enhancement of the solution accuracy of the first-stage equation, i.e., Eq. (17), is critical to the reconstruction fidelity. Normally, this would be attainable through regularizing ill-posedness of the inverse issue and incorporating a priori knowledge. It is also argued that two-component reconstruction is target-dependent and would be improved for targets of large-size and low-contrast.

 

Fig. 7. Effects of target size (a) and contrast (b) on the reconstruction quantitativeness. (a) Target contrast fixed at 3:1 for a varying radius; (b) Target radius fixed at 4 mm for varying contrast. The same phantom as in Fig. 4 is used with the target CCS=25 mm, and the results are shown for the peak values in the reconstructed images.

Download Full Size | PDF

It is clearly seen from the profiles of the reconstructed images, that the centers of the reconstructed targets appear to be deeper and the reconstructed sizes larger than the true ones. The latter is the intrinsic nature of diffuse light imaging modality that stems from ill-posedness of the inverse issue; while the reason for the former observation is complex and can be preliminarily ascribed to the use of the row-fashioned ART inversion. It is expected that both the adversities would be substantially suppressed with adoption of the whole-weight-matrix based inversion scheme and incorporation of the prior-knowledge based constrains.

In summary, we have presented a linear GPST methodology of TD-FMT for simultaneous reconstruction of fluorescent yields and lifetimes of multi-components, from which the algorithms for one- and two-component recovery are exemplified, respectively. The simulative validations of the algorithm for a 2-D domain have been performed for its abilities to discern differences in target size and fluorescence parameters, as well as for its noise-robustness and spatial resolution. The results have proved the effectiveness of the methodology. Also, the crucial issues on the implementation and validation of the algorithm have been discussed. We will explore phantom and in vivo experimental validations of the methodology in the on-going work.