Total light approach of time-domain fluorescence diffuse optical tomography

Andhi Marjono; Akira Yano; Shinpei Okawa; Feng Gao; Yukio Yamada

doi:10.1364/OE.16.015268

1. Introduction

Over the last decade, several researchers have investigated the use of florescent dyes as contrast agents to distinguish diseased tissues from normal ones through non-invasive optical measurements, including fluorescence diffuse optical tomography (FDOT). Recently, FDOT is rapidly evolving as a new modality to monitor molecular events in either humans or small animals [1]. FDOT is an emerging technology to image spatial distributions of fluorophores or other structures in tissues which are irradiated by a near-infrared (NIR) excitation light source and emit fluorescence light.

FDOT can image three parameters of fluorescence, i.e., fluorophore concentration, yield and lifetime, while in this paper we focus our attention only on imaging the fluorophore concentration. Images of the fluorophore concentrations provide valuable information of the location and functional status of the targeted tissue [2].

Many investigations on FDOT have been conducted using continuous wave (CW) and frequency domain (FD) techniques with some numerical simulations as well as phantom and in vivo experiments [3–14]. Another measurement mode is time domain (TD) technique [15–25], where the excitation is performed using ultra-short laser pulses in conjunction with time-resolved detection. TD technique gives better results of image reconstruction than CW technique [26] particularly from the point of view of spatial resolution. TD measurement data with ultra-short laser pulses implicitly contain all modulation frequencies of FD measurement, and therefore provides the richest optical information. As in TD systems for DOT, FDOT can be performed by either the full time-resolved scheme [27–29] or the featured-data one [27, 30–32]. The full time-resolved scheme is believed to provide the best image quality such as the quantification and spatial resolution of the reconstructed images at large computation cost as well as high noise-sensitivity [28], while the featured-data schemes are fast in computation and robust to noise but with degraded image quality due to the loss of the information and the correlation among the data-types [28, 30–32]. In order to obtain the images of distributions of fluorophore concentrations and internal optical properties of tissue with high quality, accurate models of the forward problem to solve the propagation of excitation and emission lights in tissues are very essential, because the errors in the forward solution will be enhanced in the inverse process and will lead to low quality of the reconstructed images. Also, a fast computation of the forward problem is highly desirable for reconstructing the image by the inversion process. The total light approach proposed by Schulz et al. [33] and further developed by Marjono et al. [34,35] has a capability of fast computation by avoiding the emission light computation using convolution calculation during the reconstruction process.

The total light approach can generally be applied to CW, FD and TD schemes. If the total light approach is employed in TD-FDOT, the convolution calculation to compute the emission light which is highly time-consuming in the whole reconstruction process is not necessary. This is possible because the measured emission data are deconvolved to obtain zero-lifetime emission data prior to the reconstruction process, although the lifetime function must be known a priori. By avoiding the computation of the emission light at each iteration process of image reconstruction, the computation time is significantly reduced. In order to realize a fast computation but still with good image quality, the total light approach using the mean time of flight as the data type is employed in this paper. The data-types in TD systems have been discussed by Arridge [27]. The total light approach is similar to double DOT because we have two measurement datasets for the excitation and total light measurements, and it is possible to simultaneously reconstruct the images of the fluorophore concentration and the absorption coefficient of the tissue.

The purpose of this research is to develop a TD-FDOT scheme for simultaneous reconstruction of the fluorophore concentration and absorption coefficient of tissues by use of the total light approach which brings high qualities of the reconstructed images with fast computation. We will start with the fundamental equations for the forward model in Section 2, followed by the inverse model in Section 3. Then simulation models and solving methods are described in Section 4, the results are shown and discussed in Section 5, and finally we summarize our concluding remarks in Section 6.

2. Forward model - formulation using total light

2.1. Photon transport model – formulation using total light

The propagation of both the excitation and emission lights in heterogeneous turbid medium can be modeled by the coupled time-domain photon diffusion equations as follows [36–39]

[- \nabla \cdot D_{x} (r) \nabla + µ_{ax} (r) + \frac{1}{c} \frac{\partial}{\partial t}] Φ_{x} (r, t) = q_{x} {(r, t)}_{,}

[- \nabla \cdot D_{m} (r) \nabla + μ_{am} (r) + \frac{1}{c} \frac{\partial}{\partial t}] Φ_{m} (r, t) = \frac{γ ε N (r)}{τ} \int_{0}^{t} Φ_{x} (r, t^{'}) e^{- \frac{t - t^{'}}{τ}} d t_{,}^{'}

where r [mm] is the position vector, t [ps] time, c [0.225 mm/ps] the speed of light, Φ [W/mm²] the diffuse photon fluence rate in a tissue volume Ω, q [W/mm³] the light source inside the medium, D [mm] the optical diffusion coefficient represented by 1/(3µ _s′) [40] where µ _s ^′ [mm^-1] the reduced scattering coefficient, and µ _a [mm^-1] the absorption coefficient, ε [mm^-1µM^-1] the molar extinction coefficient of the fluorophore, N [µM] the fluorophore concentration inside the object, γ the quantum yield, and τ [ps] the fluorescence lifetime, and the subscripts x and m indicating the excitation and emission wavelengths, respectively. As shown in Eq. (2), the emission light source is given by the convolution of the fluence rate of the excitation light Φ_x(r,t) and the lifetime function $T (t) = \frac{1}{τ} e^{- \frac{t}{τ}} .$ . Fluorophores absorb light at the excitation wavelength with the absorption coefficient of ε N(r) [33, 41]. Assuming that the optical properties of the background tissue at the excitation and emission wavelengths are equal, µ _ax(r)=µ _am(r)=µ _a (r) [18], and that all the fluorescence photons are emitted from fluorophores instantaneously after absorption of the excitation photons, meaning that the lifetime is zero, then Eqs. (1) and (2) for propagation of excitation and emission lights can be modified as the following, [33–35]

[- \nabla \cdot D (r) \nabla + µ_{at} (r) + \frac{1}{c} \frac{\partial}{\partial t}] Φ_{x} (r, t) = q_{x} (r, t)

[- \nabla \cdot + D (r) \nabla + µ_{a} (r) + \frac{1}{c} \frac{\partial}{\partial t}] Φ_{m}^{*} (r, t) = γ ε N (r) Φ_{x} (r, t)

where µ _at(r)=µ _a (r)+ε N(r) and Φ ^* _m(r,t) is the diffuse photon fluence rate for the zero lifetime emission light. Then the photon diffusion equation for the total light, Φ_t(r,t), is obtained as the following,

[- \nabla \cdot D (r) \nabla + µ_{a} (r) + \frac{1}{c} \frac{\partial}{\partial t}] Φ_{t} (r, t) = q_{x} (r, t)

where

Φ_{t} (r, t) = \frac{1}{γ} Φ_{m}^{*} (r, t) + Φ_{x} (r, t)

Equations (3) and (5) can be solved independently for Φ_x(r,t) and Φt(r,t), respectively, and Φ* _m(r,t) is obtained from Eq. (7),

Φ_{m}^{*} (r, t) = γ (Φ_{t} (r, t) - Φ_{x} (r, t)) .

The diffuse photon fluence rate of the emission light with a finite lifetime of fluorescence, Φ_m(r,t), is obtained by convolving Φ* _m(r,t) with the lifetime function, T(t), as expressed by Eq. (8) under the assumption that the fluorescence lifetime, τ, is homogeneous throughout the object, [41]

Φ_{m} (r, t) = \frac{1}{τ} \int_{0}^{t} Φ_{m}^{*} (r, t^{'}) e^{\frac{- (t - t^{'})}{τ}} d t^{'}

The excitation light source of an impulse light is incident on the position at the surface, r ₀, at time zero, and expressed by a delta function, q _x(r,t)=δ(r-r ₀, 0). The strength of the input power of the impulse light source is assumed to be unity (normalized).

Equations (3) and (5) need boundary and initial conditions for being solved. The boundary conditions are derived from the energy balance at the surface, and described by the Robin boundary conditions as Eq. (9), [42,43]

- D (r_{b}) \nabla_{n} Φ_{v} (r_{b}, t) = \frac{1}{2 A} Φ_{v} (r_{b}, t), (ν = x, t)

where r _b is the position of the surface, ∇_n denotes the gradient outward normal to the surface, A=(1+R _f)/(1 - R _f) is a factor for internal reflection, and the subscript v stands for x or t. R _f is the internal reflection coefficient at the surface and is given by Eq. (10) [44]

R_{f} = 1 . 440 n^{- 2} + 0.710 n^{- 1} + 0.668 + 0 . 0636 n,

where n is the ratio of the refractive index of the object to that of the surroundings, and is given as 1.54 in this study. The initial conditions for Φv(r,t) (ν = x, t) is given as follows,

Φ_{v} (r, t) = 0 for t \leq 0, (ν = x, t)

Φ_x(r, t) and Φ_t(r, t) are obtained by solving Eqs. (3) and (5), respectively, under the boundary and initial conditions, Eqs. (9) and (11). Then Φ* _m(r, t) is calculated from Eq. (7), and Φ_m(r, t) is obtained by use of Eq. (8).

Measured intensities of the excitation and emission lights at the surface of the object, Ψ_x(r _b,t) and Ψ_m(r _b,t) [W/mm²], are given by the fluxes of their fluence rates as Eq. (12)

Ψ_{v} (r_{b}, t) = - D (r_{b}) \nabla_{n} Φ_{v} (r_{b}, t), (ν = x, m)

2.2. Generating simulated measurement datasets

Simulated measurement datasets are obtained for the models with the given distributions of optical and fluorescence properties. By denoting ζ _s (s = 1, 2, …, N _s) and ξ _d (d = 1, 2, …, N _d) (N _s and N _d indicating the numbers of source and detectors, respectively) as the source and detector positions at the model surface, respectively, the measurement excitation and emission data are described as Ψ_x(ξ _d, ζ _s, t) and Ψ_m(ξ _d, ζ _s, t), respectively. Ψx(ξ _d, ζ _s, t) is obtained by solving Eq. (3) for Φ_x(r, t) and by substituting Φ_x(r, t) into Eq. (12). Ψ_m(ξ _d, ζ _s, t) is obtained by solving Eq. (5) for Φ_t(r, t), calculating Φ_m(r, t) from Eqs. (7) and (8), and substituting Φ_m(r, t) into Eq. (12).

3. Inverse process

The inverse process for reconstruction of the fluorophore concentration from the simulated measurement datasets of Ψ_x(ξ _d, ζ _s, t) and Ψ_m(ξ _d, ζ _s, t) is explained in the following.

Fig. 1. Normalized deconvolution result of zero-lifetime emission light after five iterations with the measured emission light and lifetime function.

Download Full Size | PDF

3.1. Setup of the measurement datasets for reconstruction

As the first step of the reconstruction process, from a simulated measurement dataset of the excitation and emission fluxes at the boundary, Ψ_x(ξ _d, ζ _s, t) and Ψ_m(ξ _d, ζ _s, t), a simulated measurement dataset of the emission light for zero-lifetime at the boundary, Ψ* _m(ξ _d, ζ _s, t) must be calculated by deconvolving Ψ_m(ξ _d, ζ _s, t) with the lifetime function, T(t), using the iterative technique of Van Cittert’s method given by the following equation [45–47]

Ψ_{m, k}^{*} = Ψ_{m, k - 1}^{*} + (Ψ_{m} - Ψ_{m, k - 1}^{*} * T)

where k represents the iteration number and * denotes convolution. The dependences of Ψ* _m, Ψ_m and T on ξ _d, ζ _s and t are not explicitly written in Eq. (13). The first guess of Ψ* _m is given as Ψ_m. Figure 1 shows the lifetime function, T(t), measured emission light, Ψ_m(t), and zero-lifetime emission light, Ψ* _m(t), which is obtained by deconvolving Ψ_m(t) with T(t) for the source and detector at the positions No. 7 and No. 10, respectively. Ψ* m(t) in Fig. 1 agrees very well with Ψ* _m(t) calculated from Φ* _m(t) of Eq. (7) in the forward calculation process without deconvolution. If Ψ* _m(t) obtained by the forward calculation is plotted in Fig. 1, it is almost perfectly superposed on the curve of Ψ* _m(t) obtained by deconvolution although very small differences between them are found in the period after 1000 s, thus showing that the deconvolution process is highly reliable and performed accurately. The deconvolution calculation for one curve of measured Ψ_m(t) took only 0.062 s for computation. Therefore, the total computation time for deconvolution all the 12×12 curves of Ψ_m(t) took 8.88 s which was necessary only once before the iterative reconstruction process started. We have calculated the fluxes for the period of 5000 ps from the impulse incidence. As the second step, the simulated dataset of the total light, Ψ_t(t), must be calculated by use of a relation similar to Eq. (6) with Φ replaced by Ψ. The computations were performed using DELL XPS 630 Intel Core 2 Duo computer running Windows XP having 2.4 GHz processor 4.00 GB of RAM.

In this study, we employ not the full time-resolved data but the mean time-of-flight (TOF) data for the inversion process. Then the measurement datasets expressed by the vectors M _x(ξ _d, ζ _s) and M _t(ξ _d, ζ _s) are given by Eq. (14),

M_{ν} (ξ_{d}, ς_{s}) = \frac{\int_{0}^{\infty} t Ψ_{ν} (ξ_{d}, ς_{s}, t) d t}{\int_{0}^{\infty} Ψ_{ν} (ξ_{d}, ς_{s}, t) d t} (ν = x, t)

In real measurements, measured instrumental response function (IRF) that can be some hundreds of picosecond must be considered. In order to be consistent with the ideal excitation light source of δ-function, the measured measurement datasets are deconvolved with IRF. When the mean TOF scheme is employed this deconvolution is equivalent to subtracting the mean-time of IRF from the mean TOF.

3.2 Inversion process for reconstruction of fluorophore concentration

Figure 2 shows the reconstruction process in our simulation. Once we obtain the measurement datasets of the mean time-of-flight of the excitation and total lights, M _x(ξ _d, ζ _s) and M _t(ξ _d, ζ _s), we start the inversion process for reconstructing µ _at(r) and µ _a(r) independently.

Fig. 2. Diagram of reconstruction diagram process of fluorescence diffuse optical tomography using total light approach.

Download Full Size | PDF

We allocate the absorption coefficients, µ _at(r) and µ _a(r), are given by allocating the absorption coefficients at all the nodes (from 1 to M) of the finite elements; i.e., µ _a(r)=(µ _a)^T u(r) and µ _at(r)=(µ_at)^T u(r), where µ_a=[µ _a ₁, µ _a ₂, …, µ _aM]^T, µ_at=[µ _at ₁, µ_at ₂, …, µ _atM]^T and u(r)=[u ₁(r), u ₂(r), …, u _M(r)]^T are the vectors of the absorption coefficients and the basis function, respectively. The lifetime, τ, is assumed to be known and unnecessary to be reconstructed. The Jacobian matrixes, J _ν (ν = x and t), for the mean time-of-flight data type is given as follows, [28, 48, 49]

J_{v} (ξ_{d}, ζ_{s}) = \frac{J_{v}^{T^{(1)}} (ξ_{d}, ζ_{s})}{E^{(0)} (ξ_{d}, ζ s)} - M_{v} (ξ_{d}, ζ_{s}) \frac{J_{v}^{T^{(0)}} (ξ_{d}, ζ_{s})}{E^{(0)} (ξ_{d}, ζ_{s})} (ν = x, t)

where

J_{ν}^{T^{(0)}} (ξ_{d}, ζ_{s}) = {[\begin{matrix} - \int_{Ω} E_{ν}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s}) u_{1} (r) d r \\ - \int_{Ω} E_{ν}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s}) u_{2} (r) d r \\ ⋮ \\ - \int_{Ω} {E_{ν}}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s}) u_{M} (r) d r \end{matrix}]}^{T} (ν = x, t)

J_{ν}^{T^{(1)}} (ξ_{d}, ζ_{s}) = {[\begin{matrix} - \int_{Ω} [E_{ν}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(1)} (r, ζ_{s}) + E_{ν}^{(1)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s})] u_{1} (r) d r \\ - \int_{Ω} [E_{ν}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(1)} (r, ζ_{s}) + E_{ν}^{(1)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s})] u_{2} (r) d r \\ ⋮ \\ - \int_{Ω} [E_{ν}^{(0)} (ξ_{d}, r) Ψ_{ν}^{(1)} (r, ζ_{s}) + E_{ν}^{(1)} (ξ_{d}, r) Ψ_{ν}^{(0)} (r, ζ_{s})] u_{M} (r) d r \end{matrix}]}_{(v = x, t)}^{T}

E_{ν}^{(i)} (ξ_{d}, ζ_{s}) = \frac{1}{2 A} Ψ_{ν}^{(i)} (ξ_{d}, ζ_{s}), (i = 0, 1 and ν = x, t)

Ψ_{ν}^{(i)} (ξ_{d}, ζ_{s}) = \int_{0}^{\infty} t^{i} Φ_{v} (ξ_{d}, ζ_{s}, t) d t, (i = 0, 1 and ν = x, t)

The absorption coefficients are updated by the iteration as χ _k+1=χ _k+δ _χk. Here, χ _k denotes the vectors of the absorption coefficients at the k ^th iteration, µ _a or µ _at, at the nodes of the FEM meshes, and their perturbations, δ _χk, are calculated by Algebraic Reconstruction Technique (ART),

δ_{χ_{k}^{j + 1}} = δ_{χ_{k}^{j}} + \frac{[h^{(j)} - J_{ν}^{(j)} (χ_{k}) \cdot δ_{χ_{k}^{j}}]}{∥ J_{ν}^{(j)} {(χ_{k}) ∥}^{2}} {[J_{ν}^{(j)} (χ_{k})]}^{T} (ν = x, t)

where the superscript j for δ _χk is the updating number and the updating continues until j = N _d×N _s. δ _χk satisfies the equation M _ν - F _ν(χ _k)=J _v δ _χ _k where F is the operator for the forward problem solving the photon diffusion equation. h ^(j) is the j ^th element of the vector [M _ν - F _ν (χ _k)], and J _v ^(j)(χ _k) is the j ^th row of the Jacobian matrix J _v(χ _k). We do not employ any relaxation parameter in this study. As shown in Fig. 2, after we reconstruct the distributions of µ _at(r) and µ _a(r), we finally can calculate the fluorophore concentration, N(r)=(1/ε) [µ _at(r) - µ _a(r)], meaning that the optical properties and the fluorophore concentration are simultaneously reconstructed. Unless the total light approach is employed, the convolution process for the emission light expressed by the right hand side of Eq. (2) must be done in every iteration process. This convolution calculation is highly time-consuming and the whole reconstruction process needs long computation time. The total light approach does not need this convolution calculation because the excitation and total lights are independently computed from Eqs. (3) and (5), respectively. Therefore, the computation time for the whole reconstruction process is significantly reduced to about 10% of that using the original Eqs. (1) and (2) for the same number of meshing.

4. Simulation model and solving method

The validation of the algorithm by simulation is performed on a 2-D circular domain. The 2-D geometries of the simulation objects with the diameter of 20 mm are shown in Figs. 3(a)–3(d) in the coordinate system (x, y) with the circle center at (0, 0). There is one circular fluorophore target located on the x-axis with the diameter of 4 mm and the center at (3,0) and (0,0) in Fig. 3(a) and 3(b), respectively. Figure 3(c) contains one circular fluorophore target with the diameter of 6 mm and the center at (3,0). Figure 3(d) contains two circular fluorophore targets located symmetrically with respect to the origin on the x-axis with their diameters of 4 mm and their centers at (±3,0) which means that the center-to-center separation (CCS) is 6.0 mm. Numbers around the circles denote the positions where the 12 optical probes are attached.

Fig. 3. Circular model with the diameter of 20 mm and with (a) one circular target with the diameter of 4 mm and the center at (3,0), (b) one circular target with the diameter of 4 mm and the center at (0,0), (c) one circular target with the diameter of 6 mm and the center at (3,0), (d) two circular targets with their diameters of 4 mm and the CCS of 4 mm, 6 mm, 8 mm or 10 mm.

Download Full Size | PDF

Figures 4(a)–4(d) show the FEM meshings which are used to generate the simulated measurement excitation and emission datasets for the objects shown in Figs. 3(a)–3(d), respectively. Figure 4(e) shows a homogeneous FEM meshing with 1089 nodes and 2048 triangles which is used for the reconstruction process. 12 coaxial source-detector (N _s = N _d = 12) fiber bundles are assumed to be attached on the boundary with an equal spacing.

Fig. 4. (a)–(d): meshings of the objects shown in Fig. 3 (a)–(d) for generating simulated measurement datasets. (e): meshing for image reconstruction by the inversion process.

Download Full Size | PDF

The optical properties and the fluorophore concentrations of the background and targets are listed in Table I. The diffusion coefficient is assumed to be homogeneous and known throughout the objects. We assume to use Indocyanine-green (ICG) as the fluorophore, which has the extinction coefficient ε = 0.013 mm^-1µM^-1, fluorescence lifetime, τ = 0.56 ns, and quantum yield, γ = 0.016 for the excitation/emission wavelengths of 778 nm/830 nm. [50].

The background absorption coefficient, µ _a = 0.0006 mm^-1, is equivalent to that of polyacetal which was the material used for phantoms [42], and this value is much smaller than the absorption coefficient of water at 800 nm which is about 0.0023 mm^-1. All the results shown in the following section were obtained with µ _a = 0.0006 mm^-1 for the background. Although not shown in this paper we also tried background µ _a up to 0.1 mm^-1, and we found that the results with higher background µ _a were still as good as (or even better than) the results presented in this paper with no significant change. The inherent absorption coefficient of the target is assumed to be about ten times the background µ _a. This is because targets are often diseased tissues with higher absorption than background.

The fundamental partial differential equations, Eq. (3) for Φ_x and Eq. (5) for Φ_t were numerically solved by FEM. We have calculated the temporal and spatial profiles for both the excitation and emission lights for the period of 5000 ps from the time of incidence (0 ps) of the impulse light source with a time interval Δt = 10 ps.

Table I. Optical properties and fluorophore concentrations of the medium

View Table

5. Results and discussions

5.1 Generation of measurement datasets

The measurement datasets were generated by solving the forward problem, and as the examples, the calculated spatial profiles of the diffuse photon fluence rates at the time of 750 ps in the model of Fig. 3(a) are shown in Figs. 5(a)–5(d) for the excitation, emission, zero-lifetime emission and total lights, respectively. The light source is irradiated from a source fiber located at the position No. 7 shown in Fig. 3. As for the excitation light shown in Fig. 5, the maximum value of the diffuse photon fluence rate of the excitation light at 750 ps is about 1.18×10^-6 found at around the center area. We assume that there is no fluorophore in the background medium which has an absorption coefficient of 0.0006 mm^-1, and that the target has a fluorophore concentration of 1 µM with the inherent absorption coefficient of 0.0066 mm^-1. Therefore, the total absorption coefficient of the target including the fluorophore is 0.0066 mm^-1+0.013 mm^-1 equal to 0.0196 mm^-1, resulting in the difference in the total absorption coefficient between the target and the background to be 0.0196 mm^-1-0.0006 mm^-1 = 0.0190 mm^-1. Consequently, the target area having a larger absorption coefficient absorbs the excitation light more strongly than the background as shown in Fig. 5(a). And as shown in Fig. 5(b), the maximum value of the diffuse photon fluence rate of the emission light is 6.13×10^-10 at 750 ps which is observed inside the fluorophore target as expected. In Fig. 5(c), the maximum value of the diffuse photon fluence rate of the zero-lifetime emission light is about 3.25×10^-9 at 750 ps which is observed inside the fluorophore target similarly to the emission light shown in Fig. 5(b). We can see that at the time of 750 ps the zero-lifetime emission light has spread more widely into the whole area of the object than the actual emission light. In Fig. 5(d), the maximum value of the diffuse photon fluence rate of the total light is about 1.30×10^-6 at 750 ps which is observed at around the center of the object. The total absorption coefficient of the target, µ _at = µ _a+ε N, for the excitation light in Eq. (3) is larger as much as ε N than the absorption coefficient for the total light in Eq. (5). Therefore, the target absorbs the excitation light more strongly than the total light, and consequently, as shown in Figs. 5(a) and 5(d), the diffuse photon fluence rate of the total light is slightly larger than that of the excitation light. At longer time, the spatial profiles of the diffuse photon fluence rate of the excitation and total lights tend to be circularly symmetric with respect to the center of the circle although their intensities rapidly decrease and finally vanish.

Fig. 5. Spatial profiles of the diffuse photon fluence rates in the model of Fig. 3(a), (a) the excitation light, (b) the emission light, (c) the zero-lifetime emission light, and (d) the total light, at the time of 750 ps from the incidence when irradiated from the source position No. 7 with a unit input power.

Download Full Size | PDF

Figures 6(a)–6(d) show the temporal profiles of the measured fluxes of the excitation, emission, zero-lifetime emission and total lights, respectively. The results are obtained by irradiating the light source of a unit input power from the source fiber placed at the position No. 7 and the fluxes are measured at the detector fibers placed at the positions No. 1, 4 and 10 for the object of Fig. 3(a). The excitation flux reaches the peak at the time of about 350 ps. The intensities of the measured excitation fluxes at the detectors No. 4 and 10 are equal due to the symmetry. The excitation flux measured at the detector No. 1 is smaller than those at the detectors No. 4 and 10 because the detector No. 1 located furthest from the position No. 7 and the excitation light is blocked by the target. As for the emission light, the detectors closer to the target will have larger measured fluxes. In the case of Fig. 6(b), the maximum emission flux of 8.4×10^-11 is observed at the detector No. 1 at the time of about 1000 ps. Then it slowly decreases and almost vanishes at the time longer than 4000 ps. The emission fluxes at the detectors No. 4 and 10 are the same again due to the symmetric arrangement to the source and target. The zero-lifetime emission flux in Fig. 5(c) show narrower temporal profiles than the actual emission light. Its maximum at the detector No. 1 is observed as 1.72×10^-9 at the time of about 500 ps and the zero-lifetime emission vanishes at the time of larger than 2800 ps. The total light, Ψ_t, is defined as the sum of the excitation light, Ψ_x, and the zero-lifetime emission light divided by the quantum yield, Ψ* _m/_γ. As shown in Figs. 6(a) and 6(c), Ψ_x is three orders of magnitude larger than Ψ* _m. Even after Ψ* _m is divided by γ = 0.016, Ψ* _m/_γ is still one order of magnitude smaller than Ψ_x. Therefore, Ψ_t shown in Fig. 6(d) is slightly larger than Ψ_x shown in Fig. 6(a).

Fig. 6. Temporal profiles of the fluxes in the model of Fig. 3(a), (a) the excitation light, (b) the emission light, (c) the zero-lifetime emission light, and (d) the total light, at the detectors No. 1, 4 and 10 for the source position No. 7.

Download Full Size | PDF

5.2 Results of reconstruction

Figure 7 show the results of reconstruction for the model of Fig. 3 (a). Figures 7(a)–7(d) show the reconstructed µ _at = µ _a+ε N, the reconstructed µ _a, the reconstructed N and its profile along the x-axis, respectively, and Fig. 7(e) shows the convergences of the inversion process by plotting the relative errors defined as the normalized differences between the measured and the estimated data of M _x and M _t during 40 iterations. 40 iterations took less than 100 seconds in computation. The convergence of M _x was slightly better than that of M _t. Both M _x and M _t well converged after 40 iterations, and µ _at and µ _a were simultaneously reconstructed. The image of N(r) shown in Fig. 7(c) was obtained by subtracting µ _a image from µ _at image and dividing it by ε. The quantitativeness, which is defined as the ratio of the reconstructed value to the true one of the peak fluorophore concentration, is about 80% as shown in Fig. 7(d). The position of the peak fluorophore concentration in the reconstructed image is shifted about one mm to the center from the true position. The full width at half maximum (FWHM) of the reconstructed fluorophore concentration is about 3.6 mm which is close to the true diameter of the target of 4 mm.

Figures 8(a)–8(d) show the reconstructed µ _at, the reconstructed µ _a, the reconstructed N, its profile along the x-axis and normalized errors during 40 iterations, respectively, for the model of Fig. 3(b). Figure 8(d) shows that the FWHM of the reconstructed profile is about 3.1 mm and the quantitativeness is about 52% but the position of the peak fluorophore concentration in the reconstructed image is accurately at the true position. As shown in Fig. 8(e), the inversion process converged after 40 iterations. By comparing the images in Figs. 7 and 8, the quantitativeness is found to be better for the target closer to the boundary. The target closer to the boundary is more strongly excited by the excitation light and the emitted fluorescent lights are absorbed more weakly by the medium between the target and the boundary. Therefore, the measured intensities of emission lights are stronger and the SNR of the measured data becomes better for the target closer to the boundary, and finally the fluorophore concentrations are recovered more closely to the true values.

Fig. 7. Images of reconstructed (a) µ _at [mm^-1], (b) µ _a, [mm^-1] (c) N [µM] for the model of Fig. 3 (a). (d) profile of the reconstructed N along the x-axis, and (e) graph of convergence during 40 iterations.

Download Full Size | PDF

Fig. 8. Images of reconstructed (a) µ _at, (b) µ _a, (c) N for the model of Fig. 3 (b). (d) profile of the reconstructed N along the x-axis, and (e) graph of convergence during 40 iterations.

Download Full Size | PDF

Fig. 9. Images of reconstructed (a) µ _at, (b) µ _a, (c) N for the model of Fig. 3 (c). (d) profile of N along the x-axis, and (e) graph of convergence during 40 iterations.

Download Full Size | PDF

Figure 9 shows the results for the model of Fig. 3(c). Contrary to the previous results, the convergence of M _x was slightly better than that of M _t while both well converged after 40 iterations. Figure 9(d) shows that the FWHM of the reconstructed profile is about 5 mm. The position of the peak fluorophore concentration in the reconstructed image is slightly shifted to the center, and the quantitativeness is about 97%. Although the shape of the reconstructed target is little bit deformed as seen in Fig. 9 (a) to (c), the quantitativeness is much better than those for the models of Figs. 3(a) and 3(b). This better recovery of the quantitativeness of the fluorophore concentration is mainly due to the fact that the target having a larger area containing the same concentration of the fluorophore emits fluorescence light more strongly. Therefore, the measured fluxes of emission lights are also stronger, the SNR of the measured data becomes better, and finally the fluorophore concentrations are recovered more closely to the true values.

Figure 10 shows the reconstructed image of N and its profile along the x-axis for the model of Fig. 3(d) having two targets with CCS = 6 mm. The quantitativeness of the peak fluorophore concentration is about 80% which is the same as for the model of Fig. 3(a) with single target located at the same position. The two positions of the peak fluorophore concentration are shifted to the center similarly to that for the single target model. The two targets are not perfectly separated but clearly seen to be resolved for CCS = 6 mm.

Fig. 10. (a). reconstructed image of N [µM] and (b) profile of N along the x-axis for the model of Fig. 3 (d) with CCS = 6 mm.

Download Full Size | PDF

We have also done the reconstruction of the model to Fig. 3(d) with CCS = 8 mm as shown in Fig. 11. The quantitativeness is recovered about 98%, this is again caused by the shorter distance between the target and the medium boundary. For a smaller CCS, we have done the reconstruction with the smallest CCS of 4 mm and the quantitativeness is recovered about 70% as shown in Fig. 12.

Fig. 11. (a). reconstructed image of N [µM] and (b) profile of N along the x-axis for the model to Fig. 3(d) with CCS = 8 mm.

Download Full Size | PDF

Fig. 12. (a). reconstructed image of N [µM] and (b) profile of N along the x-axis for the model to Fig. 3(d) with CCS = 4 mm.

Download Full Size | PDF

Fig. 13. The parameter S showing the separation between the two targets in the model of Fig. 3(d) with CCS = 4, 6, 8 and 10 mm.

Download Full Size | PDF

The separation of the two targets can be evaluated by the parameter S as the following

S = \frac{max [N (x, 0)] - N (0, 0)}{max [N (x, 0)] - min [N (x, 0)]}

where N(0,0) is the value at (x,y)=(0,0), max[N(x,0)] and min[N(x,0)] are the maximum and minimum of N on the x-axis. The results of S are 0.12, 0.53, 0.65 and 0.98 for CCS = 4 mm, 6 mm, 8 mm and 10 mm, respectively. The ideal result of S is 1.00 when the two targets are ideally reconstructed and completely separated.

5.3 Reconstruction from noisy data

Finally, we demonstrate the robustness of the algorithm to noise by using noisy data with different signal-to-noise ratio (SNR) for the model of Fig. 3(a).

Fig. 14. Reconstructed images from noisy data for the model of Fig. 3(a). SNR of the noise are (a) 55 dB, (b) 40 dB and (c) 25 dB. (d) Profiles of N along the x-axis.

Download Full Size | PDF

The noise is originally contained in the measured time-resolved fluxes. However, here we model the noise directly in each of the featured measurement data as an additive to the data with the signal-to-noise ratio (SNR, in decibels) which is related to the standard deviation of the Gaussian noise, σ, and the measurement data, M, by the equation ‖σ‖ = 10^-SNR/²⁰ ‖M‖ meaning that the norm of σ is given proportional to the norm of the measurement data M.

Figure 14 shows the reconstructed images of the fluorescence concentration for the model of Fig. 3(a) with SNR = 55 dB, 40 dB, and 25 dB and their profiles along the x-axis. SNR equal to 55 dB, 40 dB, and 25 dB are equivalent to the standard deviation of the noise equal to 0.18%, 1.00% and 5.62% of the values of the measured data. We can see that the quantitativenesses are about 78%, 68%, and 37% for SNR = 55 dB, 40 dB, and 25 dB, respectively, while it is about 80% for noiseless data as shown in Fig. 7(d). For the case of SNR = 25 dB, the artifact in the reconstructed image is noticeable and the quantitativeness is low, but the position of the peak concentration is well recovered to locate the exact position of the target. Although the results are not shown, we have found that reconstruction almost fails for SNR≤20 dB and that it is almost the same as the noiseless case for SNR≥60 dB. These results show that the noise-robustness of the algorithm is moderate.

6. Conclusions

We have presented the formulation and simulation results of the total light approach in time-domain FDOT for reconstruction of the fluorophore concentration in biological media. The validations of the algorithm have been conducted for 2-D models by showing its performances such as quantitativeness, position accuracy, separation and noise-robustness. The effectiveness of using the total light approach has been proven because this approach has reduced the computation time significantly to about 10% of those by the conventional methods and has given excellent reconstructed results of fluorophore concentration showing good quantitativeness, position accuracy, separation and noise-robustness. The total light approach in this study employed the scheme using the mean time-of-flight data which are obtained from a time-domain measurement technique. The total light approach actually can be performed using a full time-resolved data scheme but the computation time will be much longer while the result will be better than those of the featured data-type schemes. As an ongoing and future works, we are exploring real phantom and in vivo experiments for validations of the methodology.

Acknowledgments

This work was supported by Grant-in-Aid for Scientific Research (B) (KAKENHI) No. 19360098 of the Japanese government, and A. M. thanks for the support from Japanese Ministry of Education, Culture, Sports, Science and Technology. F. G. acknowledges the funding supports from the National Natural Science Foundation of China (60678049), the National Basic Research Program of China (2006CB705700).

References and links

1. J. Lewis, S. Achilefu, J. R. Garbow, R. Laforest, and M. J. Welch, “Small animal imaging: current technology and perspectives for oncological imaging,” Eur. J. Cancer 38, 2173–88 (2002). [CrossRef] [PubMed]

2. V. Ntziachristos, C. H. Tung, C. Bremer, and R. Weissleder, “Fluorescence molecular tomography resolves protease activity in vivo,” Nat. Med. 8, 757–760 (2002). [CrossRef] [PubMed]

3. M. Patterson and W. Pogue, “Mathematical model for time-resolved and frequency domain fluorescence spectroscopy in biological tissues,” Appl. Opt 33, 1963–1974 (1992). [CrossRef]

4. X. D. Li, M. A. O’Leary, D. A. Boas, B. Chance, and A. G. Yodh, “Fluorescent diffuse photon density waves in homogeneous and heterogeneous turbid media: analytic solutions and applications,” Appl. Opt 35, 3746–3758 (1996). [CrossRef] [PubMed]

5. D. J. Hawrysz and E. M. Sevick-Muraca, “Development toward diagnostic breast cancer imaging using near-infrared optical measurements and contrast agents,” Neoplasia 2, 388–417 (2000). [CrossRef]

6. V. Ntziachristos and R. Weissleder, “Experimental three-dimensional fluorescence reconstruction of diffuse media by use of a normalized born approximation,” Opt. Lett. 26, 893–895 (2001). [CrossRef]

7. E. M. Sevick-Muraca, J. P. Houston, and M. Gurfinkel, “Fluorescence-enhanced, near infrared diagnostic imaging with contrast agent,” Curr. Opin. Chem. Biol. 6, 642–650 (2002). [CrossRef] [PubMed]

8. E. E. Graves, J. Ripoll, R. Weissleder, and V. Ntziachristos, “A submillimeter resolution for small animal imaging,” Med. Phys. 30, 901–911 (2003). [CrossRef] [PubMed]

9. A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt. 42, 3081–3094 (2003). [CrossRef] [PubMed]

10. T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. 17, 545–580 (2003). [CrossRef] [PubMed]

11. A. B. Milstein, S. Oh, K. J. Webb, C. A. Bouman, Q. Zhang, D. A. Boas, and R. P. Millane, “Fluorescence optical diffusion tomography,” Appl. Opt 42, 3081–3094 (2003). [CrossRef] [PubMed]

12. A. B. Milstein, J. J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffusion tomography using multiple-frequency data,” J. Opt. Soc. Am. A 21, 1035–1049 (2004). [CrossRef]

13. S. R. Cherry, “In vivo molecular and genomic imaging: new challenges for imaging physics,” Phys. Med. Biol. 49, R13–48 (2004). [CrossRef] [PubMed]

14. V. Ntziachiristos, J. Ripoll, L. H. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotech. 23, 313–320 (2005). [CrossRef]

15. J. Wu, L. Perelman, R. R. Dasari, and M. S. Feld, “Fluorescence tomographic imaging in turbid media using early-arriving photons and Laplace transforms,” Proc. Natl. Acad. Sci. USA 94, 8783–8788 (1997). [CrossRef] [PubMed]

16. K. Chen, L. T. Perelman, Q. G. Zhang, R. R. Dasari, and M. S. Feld, “Optical computed tomography in a turbid medium using early arriving photons,” J. Biomed. Opt. 5, 144–154 (2000). [CrossRef] [PubMed]

17. D. Hattery, V. Chernomordik, M. Loew, I. Gannot, and A. Gandjbakhche, “Analytical solutions for time-resolved fluorescence lifetime imaging in a turbid medium such as tissue,” J. Opt. Soc. Am. 18, 1523–1530 (2001). [CrossRef]

18. M. Sadoqi, P. Riseborough, and S. Kumar, “Analytical models for time resolved fluorescence spectroscopy in tissues,” Phys. Med. Biol. 46, 2725–2743 (2001). [CrossRef] [PubMed]

19. D. Hall, G. Ma, F. Lesage, and Y. Wang, “Simple time-domain optical method for estimating the depth and concentration of a fluorescent inclusion in a turbid medium,” Opt. Lett. 29, 2258–2260 (2004). [CrossRef] [PubMed]

20. G. M. Turner, G Zacharakis, A. Sourbet, J. Ripoll, and V. Ntziachristos, “Complete-angle projection diffuse optical tomography by use of early photons,” Opt. Lett. 30, 409–411 (2005). [CrossRef] [PubMed]

21. G. Ma, N. Mincu, F. Lesage, P. Gallant, and L. McIntosh, “System irf impact on fluorescence lifetime fitting in turbid medium,” Proc. SPIE 5699, 263–273 (2005). [CrossRef]

22. S. Bloch, F. Lesage, L. Mackintosh, A. Gandjbakche, K. Liang, and S Achilefu, “Whole-body fluorescence lifetime imaging of a tumor-targeted near-infrared molecular probe in mice,” J. Biomed. Opt. 10, 054003 (2005). [CrossRef] [PubMed]

23. A. T. N. Kumar, J. Skoch, B. J. Bacskai, D. A. Boas, and A. K. Dunn, “Fluorescence lifetime-based tomography for turbid media,” Opt. Lett. 30, 3347–3349 (2005). [CrossRef]

24. X. Lam, F. Lesage, and X. Intes, “Time domain fluorescent diffuse optical tomography: analytical expressions,” Opt. Express 13, 2263–2275 (2005). [CrossRef] [PubMed]

25. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14, 7109–7124 (2006). [CrossRef] [PubMed]

26. S. Keren, O. Gheysens, C. S. Levin, and S. S. Gambhir, “A comparison between a time domain and continuous wave small animal optical imaging system,” IEEE Trans Med Imaging 27, 58–63 (2008). [CrossRef] [PubMed]

27. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–93 (1999). [CrossRef]

28. F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. 41, 778–791 (2002). [CrossRef] [PubMed]

29. R. Model, M. Orlt, and M. Walzel, “Reconstruction algorithm for near-infrared imaging in turbid media by means of time-domain data,” J. Opt. Soc. Am. A 14, 313–323 (1997). [CrossRef]

30. M. Schweiger and S. R. Arridge, “Application of temporal filters to time resolved data in optical tomography,” Phys. Med. Biol. 44, 1699–1717 (1999). [CrossRef] [PubMed]

31. H. Zhao, F. Gao, Y. Tanikawa, K. Homma, and Y. Yamada, “Time-resolved optical tomographic imaging for the provision of both anatomical and functional information about biological tissue,” Appl. Opt. 43, 1905–1916 (2005). [CrossRef]

32. F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” Appl. Opt 39, 5898–5910 (2001). [CrossRef]

33. R. Schulz, J. Peter, W. Semmler, and W. Bangerth, “Independent modeling of fluorescence excitation and emission with the finite element method,” in OSA Biomedical Optics Topical Meeting, Technical Digest, ThF24, OSA (2004).

34. A. Marjono, S. Okawa, F. Gao, and Y. Yamada, “Light Propagation for Time-Domain Fluorescence Diffuse Optical Tomography by Convolution Using Lifetime Function,” Optical Review 14, No. 3, 131–138 (2007). [CrossRef]

35. A. Marjono, Y. Akira, S. Okawa, F. Gao, and Y. Yamada, “Full time-resolved fluorescence diffuse optical tomography using total light approach,” in OSA Biomedical Optics Topical Meeting, Technical Digest, BMD33, OSA (2008).

36. Y. Yamada , “Light-tissue interaction and optical imaging in biomedicine,” in Annual Review of Heat Transfer, C. L. Tien , ed., (Begell House, 1995), Vol. 6, pp. 1–59.

37. M. Patterson and B. W. Pogue, “Mathematical model for time-resolved and frequency-domain fluorescence spectroscopy in biological tissues,” Appl. Opt. 33, 1963–1974 (1994). [CrossRef] [PubMed]

38. L. Hutchinson, R. Lakowicz, and M. Sevick-Muraca, “Fluorescence lifetime based sensing in tissues: a computational study,” Biophys. J. 68, 1574–1582 (1995). [CrossRef] [PubMed]

39. S. Muraca and L. Burch, “Origins of phosphorescence signals reemitted from tissues,” Opt. Lett. 19, 1928–1930 (1994). [CrossRef]

40. K. Furutsu and Y. Yamada, “Diffusion approximation for a dissipative random medium and the application,” Phys. Rev. E 50, 3634–3640 (1994). [CrossRef]

41. T. Farrel and M. Patterson , “Diffusion modelling of fluorescence in tissue,” in Handbook of Biomedical Fluorescence, M. Mycek and W. Pogue , eds., (Marcel Dekker, 2003), pp. 29–60.

42. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “Time-resolved diffuse optical tomography using a modified generalized pulse spectrum technique,” IEICE Trans. Inf. and Syst. E85-D, 133–142 (2002).

43. R. A. J. Groenhuis, H. A. Ferwerda, and J. J. T. Bosch, “Scattering and absorption of turbid material determined from reflection measurements. 1. theory,” Appl. Opt 22, 2456–2462 (1983). [CrossRef] [PubMed]

44. W. G. Egan and T. W. Hilgeman, Optical properties of inhomogeneous materials (Academic, 1979).

45. P. A. Jansson, Deconvolution with Applications in Spectroscopy (Academic Press Inc., 1984).

46. P. A. Jansson, R. H. Hunt, and E. K. Plyler, “Resolution enhancement of spectra,” J. Opt. Soc. Am. 60, 596–599 (1970).

47. V. Ntziachristos, X. Ma, A. G. Yodh, and B. Chance, “Multichannel photon counting instrument for spatially resolved near infrared spectroscopy,” Rev. Sci. Instrum . 70, 193–201 (1999). [CrossRef]

48. F. Gao, H. Zhao, and Y. Yamada, “Improvement of image quality in diffuse optical tomography by use of full time-resolved data,” Appl. Opt. 41, 778–791 (2002). [CrossRef] [PubMed]

49. F. Gao, P. Poulet, and Y. Yamada, “Simultaneous mapping of absorption and scattering coefficients from a three-dimensional model of time-resolved optical tomography,” Appl. Opt. 39, 5898–5910 (2000). [CrossRef]

50. M. Sevick-Muraca, A. Godavarty, J. Houston, A. Thompson, and R. Roy “Near-infrared imaging with fluorescent contrast agents,” in Handbook of Biomedical FluorescenceM. Mycek and W. Pogue , eds. (Marcel Dekker, 2003), pp. 445–528.

	µ _a (mm^-1)	ε (mm^-1µM^-1)	N(µM)	µ _at (mm^-1)	D(mm)
Background	0.0006	0.000	0.0	0.0006	0.39
Targets	0.0066	0.013	1.0	0.0196	0.39

Total light approach of time-domain fluorescence diffuse optical tomography

Abstract

1. Introduction

2. Forward model - formulation using total light

2.1. Photon transport model – formulation using total light

2.2. Generating simulated measurement datasets

3. Inverse process

3.1. Setup of the measurement datasets for reconstruction

3.2 Inversion process for reconstruction of fluorophore concentration

4. Simulation model and solving method

5. Results and discussions

5.1 Generation of measurement datasets

5.2 Results of reconstruction

5.3 Reconstruction from noisy data

6. Conclusions

Acknowledgments

References and links

Cited By

Figures (14)

Tables (1)

Equations (23)

Optics Express