## Abstract

In this study, time-domain fluorescence diffuse optical tomography in biological tissue is numerically investigated using a total light approach. Total light is a summation of excitation light and zero-lifetime emission light divided by quantum yield. The zero-lifetime emission light is an emitted fluorescence light calculated by assuming that the fluorescence lifetime is zero. The zero-lifetime emission light is calculated by deconvolving the actually measured emission light with a lifetime function, an exponential function for fluorescence decay. The object for numerical simulation is a 2-D 10 mm-radius circle with the optical properties simulating biological tissues for near infrared light, and contains regions with fluorophore. The inverse problem of fluorescence diffuse optical tomography is solved using time-resolved simulated measurement data of the excitation and total lights for reconstructing the absorption coefficient and fluorophore concentration simultaneously. The mean time of flight is used as the featured data-type extracted from the time-resolved data. The reconstructed images of fluorophore concentration show good quantitativeness and spatial reproducibility. By use of the total light approach, computation is performed much faster than the conventional ones.

© 2008 Optical Society of America

## 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]

where *r* [mm] is the position vector, *t* [ps] time, *c* [0.225 mm/ps] the speed of light, Φ [W/mm^{2}] the diffuse photon fluence rate in a tissue volume Ω, *q* [W/mm^{3}] 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

*µ*

*[mm*

_{a}^{-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\left(t\right)=\frac{1}{\tau}{e}^{-\frac{t}{\tau}}.$ . 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]

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,

where

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),

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]

The excitation light source of an impulse light is incident on the position at the surface, **r**
_{0}, at time zero, and expressed by a delta function, *q*
* _{x}*(

**r**,

*t*)=δ(

**r**-

**r**

_{0}, 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]

where **r**
* _{b}* is the position of the surface, ∇

*denotes the gradient outward normal to the surface,*

_{n}*A*=(1+

*R*

*)/(1 -*

_{f}*R*

*) is a factor for internal reflection, and the subscript*

_{f}*v*stands for

*x*or

*t*.

*R*

*is the internal reflection coefficient at the surface and is given by Eq. (10) [44]*

_{f}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,

Φ* _{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

^{2}], are given by the fluxes of their fluence rates as Eq. (12)

## 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*

*) and*

_{s}*ξ*

_{d}(

*d*= 1, 2, …,

*N*

*) (*

_{d}*N*

*and*

_{s}*N*

*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 Ψ*

_{d}*(*

_{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.

## 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]

where *k* represents the iteration number and * denotes convolution. The dependences of Ψ* * _{m}*, Ψ

*and*

_{m}*T*on

*ξ*

*,*

_{d}*ζ*

*and*

_{s}*t*are not explicitly written in Eq. (13). The first guess of Ψ*

*is given as Ψ*

_{m}*. Figure 1 shows the lifetime function,*

_{m}*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}*ζ*

*) and*

_{s}**M**

*(*

_{t}*ξ*

*,*

_{d}*ζ*

*) are given by Eq. (14),*

_{s}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}*ζ*

*) and*

_{s}**M**

*(*

_{t}*ξ*

*,*

_{d}*ζ*

*), we start the inversion process for reconstructing*

_{s}*µ*

*(*

_{at}*r*) and

*µ*

*(*

_{a}**r**) independently.

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}_{1},

*µ*

_{a}_{2}, …,

*µ*

*]*

_{aM}^{T}, µ

*=[*

_{at}*µ*

_{at}_{1}, µ

_{at}_{2}, …,

*µ*

*]*

_{atM}*and*

^{T}**u**(

**r**)=[

*u*

_{1}(

**r**),

*u*

_{2}(

**r**), …,

*u*

*(*

_{M}**r**)]

*are the vectors of the absorption coefficients and the basis function, respectively. The lifetime,*

^{T}*τ*, 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]

where

The absorption coefficients are updated by the iteration as *χ*
_{k+1}=*χ*
* _{k}*+

*δ*

*. Here,*

_{χk}*χ*

*denotes the vectors of the absorption coefficients at the*

_{k}*k*

^{th}iteration,

*µ*

*or*

_{a}*µ*

*, at the nodes of the FEM meshes, and their perturbations,*

_{at}*δ*

*, are calculated by Algebraic Reconstruction Technique (ART),*

_{χk}$${\delta}_{{\mathbf{\chi}}_{k}^{j+1}}=0,{\mathbf{\chi}}_{k+1}={\mathbf{\chi}}_{k}+{\delta}_{{\mathbf{\chi}}_{k}}^{{N}_{d}\times {N}_{s}}$$

$$j=0,1,2,\dots {N}_{d}\times {N}_{s}$$

where the superscript *j* for *δ*
* _{χk}* is the updating number and the updating continues until

*j*=

*N*

*×*

_{d}*N*

*.*

_{s}*δ*

*satisfies the equation*

_{χk}**M**

*-*

_{ν}*F*

*(*

_{ν}*χ*

*)=*

_{k}**J**

_{v}

*δ*

_{χ}*where*

_{k}*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**

*(*

_{ν}*χ*

*)], and*

_{k}**J**

_{v}

^{(j)}(

*χ*

*) is the*

_{k}*j*

^{th}row of the Jacobian matrix

**J**

_{v}(

*χ*

*). We do not employ any relaxation parameter in this study. As shown in Fig. 2, after we reconstruct the distributions of*

_{k}*µ*

*(*

_{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.

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*

*= 12) fiber bundles are assumed to be attached on the boundary with an equal spacing.*

_{d}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

*µ*

*= 0.0006 mm*

_{a}^{-1}for the background. Although not shown in this paper we also tried background

*µ*

*up to 0.1 mm*

_{a}^{-1}, and we found that the results with higher background

*µ*

*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.*

_{a}The fundamental partial differential equations, Eq. (3) for Φ* _{x}* and Eq. (5) for Φ

*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}*t*= 10 ps.

## 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.

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, Ψ

*, and the zero-lifetime emission light divided by the quantum yield, Ψ**

_{x}*/*

_{m}*. As shown in Figs. 6(a) and 6(c), Ψ*

_{γ}*is three orders of magnitude larger than Ψ**

_{x}*. Even after Ψ**

_{m}*is divided by*

_{m}*γ*= 0.016, Ψ*

*/*

_{m}*is still one order of magnitude smaller than Ψ*

_{γ}*. Therefore, Ψ*

_{x}*shown in Fig. 6(d) is slightly larger than Ψ*

_{t}*shown in Fig. 6(a).*

_{x}## 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

*µ*

*, the reconstructed*

_{a}*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**

*and*

_{x}**M**

*during 40 iterations. 40 iterations took less than 100 seconds in computation. The convergence of*

_{t}**M**

*was slightly better than that of*

_{x}**M**

*. Both*

_{t}**M**

*and*

_{x}**M**

*well converged after 40 iterations, and*

_{t}*µ*

*and*

_{at}*µ*

*were simultaneously reconstructed. The image of*

_{a}*N*(

**r**) shown in Fig. 7(c) was obtained by subtracting

*µ*

*image from*

_{a}*µ*

*image and dividing it by*

_{at}*ε*. 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

*µ*

*, the reconstructed*

_{a}*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.

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**

*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.*

_{t}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.

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.

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

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).

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}*/

^{20}‖

**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 Fluorescence*M. Mycek and W. Pogue
, eds. (Marcel Dekker, 2003), pp. 445–528.