Time Domain Fluorescent Diffuse Optical Tomography: analytical expressions

S. Lam; F. Lesage; X. Intes

doi:10.1364/OPEX.13.002263

1. Introduction

Optical techniques based on the Near-Infrared (NIR) spectral window have made significant progress in biomedical research in recent years. The relative low absorption and low scattering in the 600–1000 nm spectral range allow detection of photons that have traveled through several centimeters of biological tissue [1]. Coupled with accurate models of light propagation, NIR techniques enable imaging of deep tissue with boundary measurements using non-ionizing, low dose radiation.

The interest in NIR techniques is fueled by the ability of the techniques to monitor functional tissue parameters such as oxy- and deoxy-hemoglobin [2,3] and the development of appropriate low cost instrumentation [4]. Based on these qualities, NIR optical imaging is expected to play a key role in breast cancer detection, characterization [6–11] and monitoring through therapy [12]; brain functional imaging [13,] and stroke monitoring [15,16]; muscle physiological and peripheral vascular disease imaging [17]. For all these applications, NIR techniques rely on endogenous contrast such as tissue hemodynamics. Another potential application of NIR techniques is to monitor exogenous contrast [18]. Especially, we see the emergence of an optical molecular imaging field that bears great promises in clinical applications [19].

NIR fluorescence optical imaging is rapidly evolving as a new modality to monitor functional and/or molecular events in either human or animal tissue. The developments of new contrast agents that target specific molecular events [20,21] are particularly promising. By specifically binding [23,24] or being activated in tumors [25], molecular probe detection can be achieved in the early stages of molecular changes prior to structural modification [26]. Moreover, the endogenous fluorescence in the NIR spectral window is weak leading to good fluorescence sensitivity [27].

As of today, NIR molecular imaging is confined to small animal models [28] and the translation to human imaging is foreseen as imminent. However, the technical problems encountered in imaging larger tissue volumes are challenging. Besides sensitive instrumentation, robust and accurate models for fluorescent light propagation are needed. Tomographic algorithms in the continuous mode [29] and in the frequency domain [30,31] have been proposed. Both numerical and analytical models exist and have been applied successfully to experimental data.

In this work, we propose an analytical fluorescent diffuse optical algorithm extended to the time domain. For this purpose, we derive analytical solutions of the heterogeneous fluorescent diffusion equation for the 0^th, the 1^st and the 2^nd moments of the fluorescent time point spread function (FTPSF). The formalism of the fluorescent normalized Born approximation [29] is used to obtain the mathematical expressions of the forward model employed for the diffuse optical tomography (DOT) procedure. The algorithm is tested with synthetic data mimicking relevant breast geometry. These simulations highlight the advantage of incorporating higher moments in the fluorescent DOT problem.

2. Theory

In this section we describe the theoretical framework of light propagation in biological tissue and the specific theoretical frame we used to obtain the analytical solutions of the fluorescent time domain diffusion equation.

2.1. Light propagation in tissue

Light propagation in tissue is well modeled by the diffusion equation. In the time domain the mathematical expression modeling light propagation in a homogenous medium is:

\frac{1}{v} \frac{\partial}{\partial t} Φ (r, t) - D Δ^{2} Φ (r, t) + μ_{a} Φ (r, t) = S (r, t)

where Φ(r, t) is the photon fluence rate, D is the diffusion coefficient expressed as D=1/3µ′_s with µ′_s being the reduced scattering coefficient, µ_a is the linear absorption coefficient, v is the speed of light in the medium and S(r, t) is the source term (assumed to be a δ function in our case).

From Eq. (1), we can estimate the value of the field in each position in the investigated medium. In turn, the knowledge of the value of the field locally allows modeling accurately the reemission of a fluorescent field by endogenous or exogenous markers. Indeed, the fluorescent field is due to excited molecules that reemit photons at a constant wavelength. This phenomenon of reemission can be modeled as a source term embedded in the medium. The propagation from these sources to the detector is then modeled in the same frame as in Eq. (1). The temporal behavior of the excited fluorophore population at a given point is expressed by [32]:

\frac{\partial}{\partial t} N_{ex} (r, t) = - \frac{1}{τ} N_{ex} (r, t) + σ \cdot Φ^{λ 1} (r, t) [N_{tot} (r, t) - 2 N_{ex} (r, t)]

where N_ex(r, t) is the concentration of excited molecules at position r and time t, N_tot(r, t) is the concentration of total molecules of fluorophores (excited or not), τ is the radiative lifetime of the fluorescent compound (sec. or nanoseconds.), σ is the absorption cross section of the fluorophore (cm²) and Φ^λ1(r, t) is the photon fluence rate (Nb photons.s^-1.cm^-2) at the excitation wavelength λ₁. Considering that the number of excited molecule is low compared to the total molecules and taking the Fourier transform yield the expression for the concentration of excited molecules:

\frac{N_{ex} (r, ω)}{τ} = \frac{σ \cdot N_{tot} (r)}{1 - i ω τ} \cdot Φ^{λ 1} (r, ω)

where ω₁ is the angular frequency associated with t.

Then, the total fluorescent field is the sum of the contributions of all the secondary fluorescent sources over the entire volume. In the case of a point source located at r_s, the fluorescent field detected at a position r _d is modeled by:

Φ^{λ 2} (r_{s}, r_{d}, ω) = η ∭_{volume} N_{ex} (r, ω) \cdot Φ^{λ 2} (r, r_{d}, ω) \cdot d^{3} r

where Φ^λ2(r, r _d, ω) represent a propagation term of the fluorescent field from the element of volume at r to the detector position r _d et the reemission wavelength λ₂. Then, by using Eq. (4) we obtain the complete expression:

Φ^{λ 2} (r_{s}, r_{d}, ω) = ∭_{volume} Φ^{λ 1} (r_{s}, r, ω) \cdot \frac{Q_{eff} \cdot N_{tot} (r)}{1 - i ω τ} \cdot Φ^{λ 2} (r, r_{d}, ω) \cdot d^{3} r

where Q_eff=q·η.σ is the quantum efficiency, product of q the quenching factor, η the quantum yield and σ the absorption cross section of the fluorophore, Φ^λ1(r, t) is the photon fluence rate at the excitation wavelength λ₁ and at time t, and τ is the radiative lifetime of the fluorescent compound. Note that the product σ·N_tot(r) corresponds to the absorption coefficient of the fluorochrome and can be expressed also as ε·C_tot(r) where ε is the extinction coefficient of the fluorophore (cm^-1. Mol^-1) and C_tot(r) the concentration of the fluorochrome (Mol) at a position r.

2.2. Fluorescent moment analytical expression.

Following the derivation of Eq. (5) performed by O’Leary [33], Ntziachristos and Weissleder [29] proposed an approach to fluorescent diffuse optical tomography less susceptible to system source-detector couplings and tissue inhomogeneities. They cast the forward model in the frame of the normalized first order Born approximation that is mathematically expressed as:

\frac{Φ^{λ 2} (r_{s}, r_{d}, ω)}{Φ_{0}^{λ 1} (r_{s}, r_{d}, ω)} = \frac{1}{Φ_{0}^{λ 1} (r_{s}, r_{d}, ω)} ∭_{volume} Φ_{0}^{λ 1} (r_{s}, r, ω) \cdot \frac{Q_{eff} \cdot C_{tot} (r)}{1 - i ω τ} \cdot Φ_{0}^{λ 2} (r, r_{d}, ω) \cdot d^{3} r

The difference between Eq. (5) and (6) resides in the normalization achieved with the homogeneous excitation field reaching the detector. The gain attained here is that the left hand side can be determined purely from measurements thereby canceling any source-detector couplings. Following the expression of M. O’Leary [33], this expression is used to construct the forward model for DOT and then the analytical expression of the weight function is:

\frac{Φ^{λ 2} (r_{s}, r_{d}, ω)}{Φ_{0}^{λ 1} (r_{s}, r_{d}, ω)} = \frac{D^{λ 1}}{G^{λ 1} (r_{s}, r_{d}, ω)} \sum_{voxels} \frac{1}{D^{λ 1}} G^{λ 1} (r_{s}, r_{v}, ω) \cdot \frac{Q_{eff} \cdot C_{tot} (r_{v})}{1 - i ω τ} \cdot \frac{1}{D^{λ 2}} G^{λ 2} (r_{v}, r_{d}, ω) \cdot h^{3}

where $G^{λ j} (r_{1}, r_{2}, ω) = \frac{e^{({ik}_{j} ∣ r_{1} - r_{2} ∣)}}{∣ r_{1} - r_{2} ∣}$ is the system’s Green’s function with $k_{j}^{2}$ =(-v $μ_{a}^{λj}$ +iω)/D^λj, at the considered wavelength λ_j∈[λ₁, λ₂].

Fig. 1. Typical TPSF and respective moments. The IFS curve corresponds to a typical instrument response function.

Download Full Size | PDF

The Eq. (7) is defined in the frequency domain. We propose to find similar analytical expressions in the time domain. Such analytical solutions for the absorption case have been proposed in the past for the 0^th, 1^st and 2^nd moment of the TPSF [34]. The correspondence of these moments to the TPSF is illustrated in Fig. 1. The 0^th moment corresponds to the integration of the counts (equivalent to the continuous mode), the 1^st moment corresponds to the mean time of arrival of the photon and the 2^nd moment to the variance of arrival of the photon.

The normalized moments of order k>1 of a distribution function p(t) are defined by [35]:

m_{k} = 〈 t^{k} 〉 = \int_{- \infty}^{+ \infty} t^{k} \cdot p (t) dt ⁄ \int_{- \infty}^{+ \infty} p (t) dt

In our notations, for k=0, we use the simple integration without division. We employed this formalism in the case of the normalized first order born approximation. Hence the normalized 0^th moment is expressed as :

m_{0}^{λ 2} (r_{s}, r_{d}) = Φ_{N}^{λ 2} (r_{s}, r_{d}, ω = 0) = \sum_{voxels} \frac{G^{λ 1} (r_{s}, r_{v}, ω = 0) \cdot G^{λ 2} (r_{v}, r_{d}, ω = 0)}{G^{λ 1} (r_{s}, r_{d}, ω = 0)} \times \frac{Q_{eff} h^{3}}{D^{λ 2}} \times C_{tot} (r_{v})

This expression corresponds to Eq. (7) for the continuous mode. From now on, we assume that the excitation and emission wavelengths are similar (λ₁=λ₂=λ), i.e. typically the shift incurred under the fluorescent process will not be significant and the same propagator describes the propagation of light for both excitation and emission. Then normalizing the 1^st and the 2^nd moment to this first moment yields the analytical solutions: Normalized 1^st moment

m_{0}^{λ 2} (r_{s}, r_{d}) \times m_{1}^{λ 2} (r_{s}, r_{d}) = \sum_{voxels} {\begin{matrix} (τ + \frac{∣ r_{s} - r_{v} ∣ + ∣ r_{s} - r_{d} ∣}{2 . v \sqrt{μ_{a}^{λ} D^{λ}}} - \frac{∣ r_{v} - r_{d} ∣}{2 . v \sqrt{μ_{a} D^{λ}}}) \times \\ \frac{G^{λ} (r_{s}, r_{v}, ω = 0) \cdot G^{λ} (r_{v}, r_{d}, ω = 0)}{G^{λ} (r_{s}, r_{d}, ω = 0)} \times \frac{Q_{eff} h^{3}}{D^{λ}} \times C_{tot} (r_{v}) \end{matrix}}

Normalized 2^nd moment

m_{0}^{λ 2} (r_{s}, r_{d}) \cdot m_{2}^{λ 2} (r_{2}, r_{d}) = \sum_{voxels} {\begin{matrix} (\begin{matrix} τ^{2} + \frac{∣ r_{s} - r_{v} ∣ + ∣ r_{s} - r_{d} ∣}{4 . v^{2} μ_{a}^{λ} \sqrt{μ_{a}^{λ} D^{λ}}} + \frac{∣ r_{v} - r_{d} ∣}{4 . v^{2} μ_{a}^{λ} \sqrt{μ_{a}^{λ} D^{λ}}} \\ + {τ + \frac{∣ r_{s} - r_{v} ∣}{2 . v \sqrt{μ_{a}^{λ} D^{λ}}} + \frac{∣ r_{v} - r_{d} ∣}{2 . v \sqrt{μ_{a}^{λ} D^{λ}}}}^{2} - t^{- λ 2} (r_{s}, r_{d}) \cdot {τ + \frac{∣ r_{s} - r_{v} ∣ + ∣ r_{v} - r_{d} ∣}{2 . v \sqrt{μ_{a}^{λ} D^{λ}}} - \frac{∣ r_{v} - r_{d} ∣}{2 . v \sqrt{μ_{a}^{λ} D^{λ}}}} \end{matrix}) \\ \times \frac{G^{λ} (r_{s}, r_{v}, ω = 0) \cdot G^{λ} (r_{v}, r_{d}, ω = 0)}{G^{λ} (r_{s}, r_{d}, ω = 0)} \times \frac{Q_{eff} h^{2}}{D^{λ}} \times C_{tot} (r_{v}) \end{matrix}}

Where t^-λ2(r_s,r_d) corresponds to the fluorescent mean time for the particular source-detector pair considered. For simplicity, we do not present in this paper the full derivation of these analytical solutions.

The fluorescent DOT problem in the time domain is based on the analytical expression derived above and summarized in the set of linear equations:

∣ \begin{matrix} m_{0}^{λ 2} (r_{s 1}, r_{d 1}) \\ ⋮ \\ m_{0}^{λ 2} (r_{sm}, r_{dm}) \\ m_{0}^{λ 2} (r_{s 1}, r_{d 1}) \cdot m_{1}^{λ 2} (r_{s 1}, r_{d 1}) \\ ⋮ \\ m_{0}^{λ 2} (r_{sm}, r_{dm}) \cdot m_{1}^{λ 2} (r_{sm}, r_{dm}) \\ m_{0}^{λ 2} (r_{s 1}, r_{d 1}) \cdot m_{2}^{λ 2} (r_{s 1}, r_{d 1}) \\ ⋮ \\ m_{0}^{λ 2} (r_{sm}, r_{dm}) \cdot m_{2}^{λ 2} (r_{sm}, r_{dm}) \end{matrix} ∣ = ⋮ ∣ \begin{matrix} W_{11}^{m_{0}^{λ 2}} & \dots & W_{\ln}^{m_{0}^{λ 2}} \\ ⋮ & ⋱ & ⋮ \\ W_{ml}^{m_{0}^{λ 2}} & \dots & W_{mn}^{m_{0}^{λ 2}} \\ W_{11}^{m_{0}^{λ 2} \cdot m_{1}^{λ 2}} & \dots & W_{\ln}^{m_{0}^{λ 2} \cdot m_{1}^{λ 2}} \\ ⋱ & ⋮ \\ W_{m 1}^{m_{0}^{λ 2} \cdot m_{1}^{λ 2}} & \dots & W_{mn}^{m_{0}^{λ 2} \cdot m_{1}^{λ 2}} \\ W_{11}^{m_{0}^{λ 2} \cdot m_{2}^{λ 2}} & \dots & W_{\ln}^{m_{0}^{λ 2} \cdot m_{2}^{λ 2}} \\ ⋮ & ⋱ & ⋮ \\ W_{m 1}^{m_{0}^{λ 2} \cdot m_{2}^{λ 2}} & \dots & W_{mn}^{m_{0}^{λ 2} \cdot m_{2}^{λ 2}} \end{matrix} ∣ \cdot ∣ \begin{matrix} C_{tot} (r_{vl}) \\ ⋮ \\ ⋮ \\ ⋮ \\ ⋮ \\ ⋮ \\ ⋮ \\ ⋮ \\ ⋮ \\ C_{tot} (r_{vn}) \end{matrix} ∣

where $W_{ij}^{m_{0}^{λ 2}}$ , $W_{ij}^{m_{0}^{λ 2} \cdot m_{1}^{λ 2}}$ and $W_{ij}^{m_{0}^{λ 2} \cdot m_{2}^{λ 2}}$ , the weight function for the i^th source-detector pair and the j^th voxel are directly derived respectively from Eqs. (9), (10) and (11). In this inverse problem, the object function is defined as the fluorophore concentration. For the cases presented herein, we implemented boundary conditions using the extrapolated boundary conditions [36] and image sources.

2.3. Inverse problem

Many different approaches exist to tackle the inverse problem. In this work we choose to employ the algebraic reconstruction technique (ART) due to its modest memory requirements for large inversion problems and the calculation speed it attains.

Algebraic techniques are well known and broadly used in the biomedical community [37]. These techniques operate on a system of linear equations such as the ones seen in Eq. (12). We can rewrite Eq. (12) as:

b = A \cdot x

where b is a vector holding the measurements for each source-detector pair, A is the matrix of the forward model (weight matrix), and x is the vector of unknowns (object function). ART solves this linear system by sequentially projecting a solution estimate onto the hyperplanes defined by each row of the linear system. The technique is used in an iterative scheme and the projection at the end of the k^th iteration becomes the estimate for the (k+1)^th iteration. This projection process can be expressed mathematically as [38]:

x_{j}^{(k + 1)} = x_{j}^{(k)} + ξ \frac{b_{i} - \sum_{i} a_{ij} x_{j}^{(k)}}{\sum_{i} a_{ij} a_{ij}} \sum_{i} a_{ij}

where $x_{j}^{(k)}$ is the k^th estimate of j^th element of the object function, b_i the i^th measurement, aij the i-j^th element of the weight matrix A and ξ the relaxation parameter.

The relaxation parameter adjusts the projection step for each iteration. A small ξ value makes the inversion more robust but also slows down convergence. The selection of ξ is most of the time, done empirically [39]. We have chosen ξ=0.1 based on previous studies [40]. Also, a positive constraint was imposed on the object function. This hard constraint is adequate with fluorescent measurements as long as negative concentrations are unphysical. The stopping criterion of the iterative inversion algorithm was set a posteriori using the reconstructions of Fig. 5 for the Cy 5.5 case. The criterion retained was the number of iterations, which was fixed to 500 and selected using the merit functions described in [40].

2.3. Simulations

We tested the formulation derived in Section 2.2 with simulations. First we constructed a synthetic phantom with parameters relevant to the softly compressed human breast in dimension (6cm thickness) and for the optical endogenous properties. Second we simulated a homogeneous fluorochrome distribution over the volume with 1 cm³ heterogeneities exhibiting a contrast of 10 in concentration. The different parameters of the simulations are provided in Table 1.

The fluorescent signal is dependent on the intrinsic characteristics of the fluorochrome employed. Simulations were carried out with three representative compounds: Cy 7, Cy 5.5 and Cy 3B. These fluorochromes were selected due to the span of lifetimes they do exhibit, which is characteristic of cyanine dyes [22]. The different properties of these fluorochromes are provided in Table 2.

Table 1. Parameters used in the simulations

View Table | View all tables in this article

Fig. 2. Configuration used for the simulations herein. The source (detectors) locations are depicted by red (blue) dots

Download Full Size | PDF

The synthetic phantom was probed with a 25×25 constellation of source detectors. This constellation was distributed evenly 1.5 cm apart in both dimensions. The phantom configuration is provided in Fig. 2.

Table 2. Fluorochrome investigated herein.

View Table | View all tables in this article

2.3. Noise model

Higher order moments are sensitive to noise. Thus, it is important to investigate the performance of the algorithm in the presence of noise. Analytical noise models exist for the intrinsic NIR higher moments for homogeneous cases [35]. However, the derivation of the same analytical model for tomographic purposes is overly complex. We decided thus to employ a heuristically derived noise model.

We generated synthetic homogeneous TPSF and considered a Poisson noise of the temporal distribution of photon time of flights. The TPSF was normalized at 500 counts at the maximum bin mimicking real acquisition scenarios. From the noised TPSF, we estimated one set of energy, meantime and variance. The same estimation was performed over 1,000 trials. The statistics of these estimates were used as our noise model. An example of noisy moments value distribution is given in Fig. 3.

Fig. 3. Results of repartition of energy, meantimes and variance of 1,000 randomly generated noised TPSF.

Download Full Size | PDF

Gaussian distribution approximated the noise model. The different values of the noise model employed for the three moments evaluated herein is given in Table 3.

Table 3. Noise model used. The standard deviations are expressed in percent of the mean.

View Table | View all tables in this article

3. Results

In this section we present preliminary results concerning fluorescent DOT based on the formulation by normalized moments.

3.1. Sensitivy matrix

The pattern described by the photon migration is often referred as a banana shape. Especially, in the case of continuous mode, the measurements are highly sensitive to the surface. Such dependence of the data type can be visualized through the mapping of the sensitivity matrix. Indeed, each line of the linear system described in Eq. (12) represents the sensitivity to a local perturbation for the correspondent source-detector pair. Thus by mapping this local dependence, we render the spatial sensitivity of this particular source-detector pair for this specific configuration and specific data type.

We propose in Fig. 4 some examples of sensitivity matrices for the transmittance case. We limited ourselves to depict slices across the discrete volume, but by construction, the banana shapes are in 3D. The optical and fluorochrome properties characterizing this medium are provided in Table 1 and Table 2.

Fig. 4. Example of sensitivity matrices. a) and b) correspond respectively to $m_{0}^{λ 2}$ (r _s,r _d) and $m_{0}^{λ 2}$ (r _s,r _d)· $m_{2}^{λ 2}$ (r _s,r _d) for a 6cm thick slab with source-detector facing each other and a 0.1 µM background of Cy 7. c) and d) correspond to the same parameters for a 0.1 µM background of Cy 5.5. Last, e) and f) correspond to the same parameters for a 0.1 µM background of Cy 3B.

Download Full Size | PDF

The examples in Fig. 4 underline interesting features of the time domain moment fluorescent DOT. First, as seen in Figs. 4(a),–(c) and (e), the normalized 0^th order Born approximation in continuous mode is highly sensitive to surface voxels. This is a well-known behavior that is both present in absorption and fluorescent mode. This also demonstrates the poor sensitivity of planar fluorescent techniques to deep fluorescent inclusions due to overwhelming dependence on surface interactions.

Secondly, we see that the spatial dependence profile of the 2^nd normalized fluorescent moment possesses distinctive features. The 2^nd normalized fluorescent moment still exhibits some strong dependence from the surface voxels, but also from deeper voxels. The profile presents a distinguishing depression in the line connecting the source-detector pair. This fact is striking in the case of Fig. 4(b) where we used the properties of Cy 7 for the simulated chromophore. In this specific case, the 2^nd normalized fluorescent moment is characterized by a sharp and well-demarcated hollow dependence. Such typical features are related to the fact that the fluorescent mean time t-^λ2(r _s,r _d) is subtracted in Eq. (11). Indeed, the measured mean-time is always greater than the mean time of propagation for the shorter path, i.e. for the voxels located on the line connecting the source-detector pair. Then if the contribution of the lifetime is small enough, the 2^nd normalized fluorescent moment will exhibit reduced (eventually negative) contribution for these voxels. This property is dependent on the lifetime of the fluorochrome investigated. This hollow distribution is still present for the Cy 5.5 case but disappears for the Cy 3B simulations. In this last case, the contribution of the lifetime is predominant for these shorter path voxels and the spatial distribution of 2^nd normalized fluorescent moment is not markedly different than the 0^th normalized fluorescent moment.

From this set of examples, we note that the 2^nd normalized fluorescent moment provides a different kind of information compared to the 0^Th normalized moment of the fluorescent TPSF (we overlooked here the 1^st moment for simplicity). The incorporation of this additional information in fluorescent DOT is expected to produce more accurate reconstructions. In the next section we propose an example of 3D fluorescent reconstructions validating this hypothesis.

3.2. Reconstructions

One should note that the background fluorochrome concentration is non-negligible. This background simulates non-perfect compound uptake/trapping and represents a challenging case for all FDOT approaches.

We use then the formulation of (12) to generate synthetic measurements from the phantom. We simulated a 25 source and 25 detectors array as described in Fig. 2(a). The value of the fluorescent mean time and the fluorescent variance were evaluated to be around ~3 ns and 1 ns respectively. These values are in agreement with expected values for real cases. In this simulation no noise was added. The reconstructions obtained by constrained ART are provided in Fig. 5. We propose in this figure the reconstructions based on the 0^th normalized moment of the fluorescent TPSF and with the combined three normalized moments.

In all three cases presented, the inclusions were successfully reconstructed. Their locations were well retrieved and the objects clearly discriminated from the background. However, we can notice some differences in terms of reconstruction quality between the three different inverse problems solved. Especially, in the case using only the 0^th normalized moment of the fluorescent TPSF, the reconstructions exhibit strong artifacts on the boundary, artifacts that scale with the reconstructed heterogeneity. On the other hand the reconstructions based on the three moments combined (as reconstructions based on the 2^nd normalized fluorescent moment solely; results not shown here) do not exhibit such strong surface artifacts. In this last case, the homogenous background fluorophore is more accurately reconstructed.

Theses findings are related to the results of Section 3.1. For the 0^Th normalized moment reconstructions, a high sensitivity to surface voxels leads to artifacts placed in front of the individual sources and detectors. This is emphasized in our case where a non-negligible fluorophore homogeneous background concentration was simulated. The contribution of this homogenous background to the measurements is reconstructed as strong concentrations localized in front of the optodes. Reconstructions based on the 2^nd normalized fluorescent moment suffer less from this ambiguity. In the latter case, the reconstruction does not exhibit artifact scaling with the reconstructed heterogeneity and the homogeneous background is reconstructed with more fidelity. The gain is even more substantial when the three moments are used simultaneously in the inverse problem. In this case, the object is accurately reconstructed in location and size with a more homogeneous background fluorophore concentration.

Last, the reconstructions based on the three different compounds are very similar when using only the 0^th normalized moment. However, and as expected from Section 3.1, the reconstructions employing the 2^nd normalized moment exhibit different performances. In the case of relatively short lifetimes, i.e., Cy 7and Cy 5.5, the reconstructions are similar and provide accurate recovery of the three heterogeneities. However, in the case of longer lifetime, i.e., Cy 3B, even though, the reconstructions are far superior when using the 3 moments simultaneously in the inverse problem, the objects are less well defined. This fact is linked to the close similarity between the spatial distributions of the 0^th normalized and the 2^nd normalized fluorescent moments. One should note also that the constellation of source-detector selected herein is quite sparse and such reconstructed structure is expected as seen in Ref [41].

Fig. 5. Reconstruction from synthetic data for Cy 7: a) 0th moment only, b) 0^th, 1^st and 2^nd moments; Cy 5.5 : c) 0^th moment only, d) 0^th, 1^st and 2^nd moments; and Cy 3B: e) 0^th moment only, f) 0^th, 1^st and 2^nd moments. The quantitative values are expressed in µM.

Download Full Size | PDF

3.2. Noisy reconstructions

The noise model described in section 2.3 was applied to the measurements for the Cy 5.5 case. The reconstructions based on this noisy simulation are provided in Fig. 6. We restricted the reconstruction to the Cy 5.5 case only for conciseness.

Fig. 6. Reconstruction from synthetic data for Cy 5.5 using all three moments noisy data.

Download Full Size | PDF

As one can see, the algorithm is still performing satisfactorily in the case of noise. Even though the 2^nd normalized moments are sensitive to noise, the incorporation of this information benefits the inverse problem. The objects are reconstructed with fidelity and the surface artifacts are still minimized due to the inherent spatial information of the 2^nd normalized moment.

4. Conclusion

In this contribution, we presented the analytical expressions of the 0^th, 1^st and 2^nd normalized moments of the fluorescent TPSF. To our knowledge, this is the first time that such higher moment expressions are derived for fluorescent diffuse optical tomography explicitly. The analytical expressions were tested with synthetic data. Simulated phantoms with a constant background of 0.1 µM of Cy 5.5, Cy 7 and Cy 3B bearing an inclusion of 1 µM of the same compound were created. The reconstructions based on a forward model using the 3 normalized moments produced superior results compared to the reconstructions based on individual moments. The higher moments of the fluorescent TPSF provide information that is less overwhelmed by the surface interactions. The gain is important when a background fluorophore concentration exists, as is generally the case in molecular imaging. Using this reconstruction technique, the background fluorophore distribution is not reconstructed as strong surface concentrations that are generally considered as plaguing artifacts in continuous wave fluorescent imaging. Especially, the incorporation of the 2^nd normalized moment resulted in superior reconstructions due to higher sensitivity to deep voxels even though the 2^nd moment is known to be sensitive to noise. Reconstructions based on relevant noise levels were still satisfactory. The next step will be to assess the usefulness of this forward model in simple phantoms and real case scenario.

References and links

1. A Yodh and B. Chance, “Spectroscopy and imaging with diffusing light,” Phys. Today 48, 34–40 (1995). [CrossRef]

2. X. Intes and B. Chance, “Non-PET Functional Imaging Techniques Optical,” Clin. No. Am. 43, 221–234 (2005).

3. F. Jobsis, “Noninvasive infrared monitoring of cerebral and myocardial sufficiency and circulatory parameters,” Science 198, 1264–1267 (1977). [CrossRef] [PubMed]

4. Y. Lin, G. Lech, S. Nioka, X. Intes, and B. Chance, “Noninvasive, low-noise, fast imaging of blood volume and deoxygenation changes in muscles using light-emitting diode continuous-wave imager,” Rev. Sci. Instrum. 73, 3065–3074 (2002). [CrossRef]

5. Y. Chen, C. Mu, X. Intes, D. Blessington, and B. Chance, “Frequency domain phase cancellation instrument for fast and accurate localization of fluorescent heterogeneity,” Rev. Sci. Instrum. 74, 3466–3473 (2003). [CrossRef]

6. B. Tromberg, N. Shah, R. Lanning, A. Cerussi, J. Espinoza, and T. Pham, et al., “Non-invasive in vivo characterization of breast tumors using photon migration spectroscopy,” Neoplasia 2, 26–40 (2000). [CrossRef] [PubMed]

7. D. Grosenick, H. Wabnitz, K. Moesta, J. Mucke, M. Moller, C. Stroszczunski, J. Stobel, B. Wassermann, P. Schlag, and H. Rinnerberg, “Concentration and oxygen saturation of haemoglobin of 50 breast tumours determined by time-domain optical mammography,” Phys Med. Biol. 49, 1165–1181 (2004). [CrossRef] [PubMed]

8. H. Jiang, N. Iftimia, J. Eggert, L. Fajardo, and K. Klove, “Near-infrared optical imaging of the breast with model-based reconstruction,” Acad. Radiol. 9, 186–194 (2002). [CrossRef] [PubMed]

9. M. Franceschini, K. Moesta, S. Fantini, G. Gaida, E. Gratton, and H. Jess, et al., “Frequency-domain techniques enhance optical mammography: Initial clinical results,” Proc. Nat. Acad. Sci. Am. 94, 6468–6473 (1997). [CrossRef]

10. S. Colak, M. van der Mark, G. Hooft, J. Hoogenraad, E. van der Linden, and F. Kuijpers, “Clinical optical tomography and NIR spectroscopy for breast cancer detection,” IEEE J. Sel. Top. Quatum Electron. 5, 1143–1158 (1999). [CrossRef]

11. X. Intes, S. Djeziri, Z. Ichalalene, N. Mincu, Y. Wang, P. St. -Jean, F. Lesage, D. Hall, D. A. Boas, and M. Polyzos, “Time-Domain Optical Mammography Softscan^®: Initial Results on Detection and Characterization of Breast Tumors”, Proc. SPIE 5578, 188–197 (2004). [CrossRef]

12. D. B. Jakubowski, A. E. Cerussi, F. Bevilacqua, N. Shah, D. Hsiang, J. Butler, and B. J. Tromberg, “Monitoring neoadjuvant chemotherapy in breast cancer using quantitative diffuse optical spectroscopy: a case study,” J Biomed Opt. 9, 230–238 (2004). [CrossRef] [PubMed]

13. G. Strangman, D. A. Boas, and J. Sutton, “Non-invasive neuroimaging using Near-Infrared light,” Biol. Psychiatry 52, 679–693 (2002). [CrossRef] [PubMed]

14. Y. Chen, D. Tailor, X. Intes, and B. Chance, “Quantitative correlation between Near-Infrared spectroscopy (NIRS) and magnetic resonance imaging (MRI) on rat brain oxygenation modulation,” Phys. Med. Biol. 48, 417–427 (2003). [CrossRef] [PubMed]

15. M. Stankovic, D. Maulik, W. Rosenfeld, P. Stubblefield, A. Kofinas, and E. Gratton, et al., “Role of frequency domain optical spectroscopy in the detection of neonatal brain hemorrhage- a newborn piglet study,” J. Matern. Fetal Med. 9, 142–149 (2000). [CrossRef] [PubMed]

16. J. C. Hebden, A. Gibson, T. Austin, R. M. Yusof, N. Everdell, D. T. Delpy, S. R. Arridge, J. H. Meek, and J. S. Wyatt, “Imaging changes in blood volume and oxygenation in the newborn infant brain using three-dimensional optical tomography,” Phys Med Biol. 49, 1117–1130 (2004). [CrossRef] [PubMed]

17. V. Quaresima, R. Lepanto, and M. Ferrari, “The use of near infrared spectroscopy in sports medicine,” J. Sports Med. Phys. Fitness 43, 1–13 (2003). [PubMed]

18. X. Intes, J. Ripoll, Y. Chen, S. Nioka, A. G. Yodh, and B. Chance, “In vivo continuous-wave optical breast imaging enhanced with Indocyanine Green,” Med. Phys. 30, 1039–1047 (2003). [CrossRef] [PubMed]

19. R. Weissleder and U. Mahmood, “Molecular imaging,” Radiology 219, 316–333 (2001). [PubMed]

20. J. V. Frangioni, “In vivo near-infrared fluorescence imaging,” Curr. Opin. Chem. Biol. 7, 626–634 (2003). [CrossRef] [PubMed]

21. K. Licha, “Contrast agents for optical imaging,” Topics in Current Chemistry 222, 1–29 (2002). [CrossRef]

22. G. Zheng, Y. Chen, X. Intes, B. Chance, and J. Glickson, “Contrast-Enhanced NIR Optical Imaging for subsurface cancer detection,” J. Porphyrin and Phthalocyanines 8, 1106–1118 (2004). [CrossRef]

23. S. Achilefu, R. Dorshow, J. Bugaj, and R. Rajagopalan, “Novel receptor-targeted fluorescent contrast agents for in-vivo tumor imaging,” Invest. Radiol. 35, 479–485 (2000). [CrossRef] [PubMed]

24. Y. Chen, G. Zheng, Z. Zhang, D. Blessington, M. Zhang, and H. Li, et al., “Metabolism Enhanced Tumor Localization by Fluorescence Imaging: In Vivo Animal Studies,” Opt. Lett. 28, 2070–2072 (2003). [CrossRef] [PubMed]

25. R. Weissleder, C. H. Tung, U. Mahmood, and A. Bogdanov, “In vivo imaging with protease-activated near-infrared fluorescent probes,” Nat. Biotech. 17, 375–378 (1999). [CrossRef]

26. R. Weinberg, “How Does Cancer Arise,” Sci. Am. 275, 62–71 (1996). [CrossRef] [PubMed]

27. X. Intes, Y. Chen, X. Li, and B. Chance, “Detection limit enhancement of fluorescent heterogeneities in turbid media by dual-interfering excitation,” Appl. Opt. 41, 3999–4007 (2002). [CrossRef] [PubMed]

28. 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.o Cancer 38, 2173–88 (2002). [CrossRef]

29. 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]

30. M. J. Eppstein, D. J. Hawrysz, A. Godavarty, and E. M. Sevick-Muraca, “Three-dimensional, Bayesian image reconstruction from sparse and noisy data sets: Near-infrared fluorescence tomography,” Proc. Nat. Acad. Sci. Am. 99, 9619–9624 (2002). [CrossRef]

31. 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]

32. X. Li, “Fluorescence and diffusive wave diffraction tomographic probes in turbid media,” PhD University of Pennsylvania (1996).

33. M. O’Leary, “Imaging with diffuse photon density waves,” PhD University of Pennsylvania (1996).

34. E. Hillman, “Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications,” PhD University College London (2002).

35. A. Liebert, H. Wabnitz, D. Grosenick, M. Moller, R. Macdonald, and H. Rinnerberg, “Evaluation of optical properties of highly scattering media by moments of distributions of times of flight of photons,” Appl. Opt. 42, 5785–5792 (2003). [CrossRef] [PubMed]

36. R. C. Haskell, L. O. Svaasand, T. Tsay, T. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am A 11, 2727–41 (1994). [CrossRef]

37. R. Gordon, R. Bender, and G. Herman, “Algebraic reconstruction techniques (ART) for the three dimensional electron microscopy and X-Ray photography,” J. Theoret. Biol. 69, 471–482 (1970). [CrossRef]

38. A. Kak and M. Slaney, “Computerized tomographic Imaging”, IEEE Press, N-Y (1987).

39. D. Ros, C. Falcon, I. Juvells, and J. Pavia, “The influence of a relaxation parameter on SPECT iterative reconstruction algorithms,” Phys. Med. Biol. 41, 925–937 (1996). [CrossRef] [PubMed]

40. X. Intes, V. Ntziachristos, J. Culver, A. G. Yodh, and B. Chance, “Projection access order in Algebraic Reconstruction Techniques for Diffuse Optical Tomography,” Phys. Med. Biol. 47, N1–N10 (2002). [CrossRef] [PubMed]

41. E. Graves, J. Culver, J. Ripoll, R. Weissleder, and V. Ntziachristos, “Singular-value analysis and optimization of experimental parameters in fluorescence molecular tomography,” J. Opt. Soc. Am. A 21, 231–241 (2004). [CrossRef]

$µ_{a}^{λ 1}$ (cm^-1)	0.06	Dimensions (cm)	9×6×9
$µ_{a}^{λ 2}$ (cm^-1)	0.06	C_background (µM)	0.1
µ $'_{a}^{λ 2}$ (cm^-1)	10.00	C_inclusion (µM)	1.0
µ $'_{a}^{λ 2}$ (cm^-1)	10.00	Voxel size (cm)	0.36×0.3×0.36

Compound	τ(ns)	ε (cm^-1.M^-1)	η(%)
Cy 7	<0.3	200 000	0.28
Cy 5.5	1.0	190 000	0.23
Cy 3-B	2.8	130 000	0.67

Measure	Energy	Meantime	Variance
σ(%)	2	0.2	2

$µ_{a}^{λ 1}$ (cm^-1)	0.06	Dimensions (cm)	9×6×9
$µ_{a}^{λ 2}$ (cm^-1)	0.06	C_background (µM)	0.1
µ $'_{a}^{λ 2}$ (cm^-1)	10.00	C_inclusion (µM)	1.0
µ $'_{a}^{λ 2}$ (cm^-1)	10.00	Voxel size (cm)	0.36×0.3×0.36

Compound	τ(ns)	ε (cm^-1.M^-1)	η(%)
Cy 7	<0.3	200 000	0.28
Cy 5.5	1.0	190 000	0.23
Cy 3-B	2.8	130 000	0.67

Time Domain Fluorescent Diffuse Optical Tomography: analytical expressions

Abstract

1. Introduction

2. Theory

2.1. Light propagation in tissue

2.2. Fluorescent moment analytical expression.

2.3. Inverse problem

2.3. Simulations

2.3. Noise model

3. Results

3.1. Sensitivy matrix

3.2. Reconstructions

3.2. Noisy reconstructions

4. Conclusion

References and links

Cited By

Figures (6)

Tables (3)

Equations (14)

Optics Express