Fluorescence lifetime imaging from time resolved measurements using a shape-based approach

Diego Álvarez; Paúl Medina; Miguel Moscoso

doi:10.1364/OE.17.008843

1. Introduction

Because of the strong interaction of light with biological tissue, in vivo optical imaging is showing great promise for studying cell and tissue function. Among the optical techniques for noninvasive macroscopic imaging, fluorescence techniques are becoming increasly important [1]. This is due to the use of fluorescent markers that provide higher specificity than absorption or scattering-based techniques [2, 3, 4, 5]. The biological target is labeled with a fluorescent molecule that is excited by a external light source enhancing the image and carrying information about the environment (local pH, calcium or sodium ion concentrations, oxygen saturation, ...). One approach is to use molecules with fluorescence spectra that depends on the local environment, and to analyze the spectrally resolved response of the medium. Another way of obtaining this information is using fluorescence lifetime measurements. The fluorescence lifetime of the excited molecules changes due to different quenching mechanisms, such as by collisions or energy transfer that deactivate the excited state. We will consider in this paper, this later approach: fluorescence lifetime-based imaging.

Fluorescence lifetime imaging (FLIM) involves the reconstruction of the spatial distribution of the fluorescence decay rates within the tissue samples. It uses the lifetime of the fluorophore signal, rather than its intensity, to create these images. Due to their relative simplicity, the most common techniques are based on continuous wave (CW) excitation [6, 7, 8, 9] and frequency modulated (FD) excitation [10, 11, 12, 13]. However, current research also focuses on pulsed excitation in the time domain (TD) [14, 15, 16, 17, 18, 19]. The TD approach requires a more sophisticated technology (short laser pulses for excitation and time resolved detectors) but offers several advantages over the CW and FD approaches [20]. Time-domain detection adds a new dimension to fluorescence data and the possibility to analyze more complex decay profiles in biological samples ranging from single cells to bulk tissue.

The problem of fluorescence lifetime-based imaging in several centimiter thick tissue closely resembles that of fluorescence microscopy and that of phosphorescence lifetime imaging (PLI). However, it poses several challenges compared to these later problems. In fluorescence microscopy, the tissue samples are optically very thin. Hence, the scattering is neglegible and the modeling is much simpler since there is no need to account for it. On the other hand, in phosphorescence lifetime imaging the photon migrations times in tissue (of the order of nanoseconds) are neglegible compared to the lifetimes of phosphorescence probes (of the order of tens of microseconds) leading to a simpler problem as well. However, the typical lifetimes of the fluorophores commonly used to obtain functional information of biological samples range from picoseconds to nanoseconds.

Several challenges arise in fluorescence-lifetime-based imaging in optically thick media due to the above mentioned issues: (i) one must account for multiple scattering of light in tissue giving rise to a more complex description of the light propagation phenomena, and (ii) the photon migrations times are comparable to the fluorescence decay rates of the fluorescent probes. Therefore, the temporal distribution of the resulting fluorescence signal is the convolution of the pulse dispersion in the tissue and the fluorescence decay rates. Since multiple scattering causes severe loss of spatial resolution, one cannot make use of direct images and must develope new image-reconstruction algorithms to solve the associated inverse problem that improve our visualization capability. In fact, the appearence of sophisticated fluorescent labels for specific biological targets, such that genetically expressed fluorophore labels, is motivating the development of new optical technologies for fluorescence imaging that get more information out of the fluorescent molecules. We note that in fluorescence imaging, the fluorescent probes act as internal sources, and the visualization process amounts to the solution of an inverse source problem from boundary measurements.

In this paper we propose a shape-based approach [21, 22, 23] to reconstruct the support of the fluorescent markers, together with a pixel-based gradient method to retrieve their fluorescence lifetimes. The joint search of the fluorophore shapes and the decay rates is carried out solving an optimization problem. In other words, we simultaneously estimate the fluorophore distribution using a shape based approach, and the fluorescence lifetimes using a gradient based one. Our numerical results validate this combined strategy. Instead of classical a pixel-based iterative scheme that usually lead to ’over-smoothed’ images of the optical tissue properties, a shape-based strategy yields to well defined boundaries and to more reliable estimates of the optical parameters of interest. This is so, because a shape-based strategy involves and implicit regularization that reduces the dimensionality of the inverse problem and thereby helps to stabilize the reconstruction process.

Our algorithm uses the so called ’level set technique’ for representing the shapes of the tissue regions with different fluorescence decay rates. The level set technique was originally introduced by Osher and Sethian for modeling flame propagation [24], and was applied first by Santosa to the solution of an inverse problem [25]. Since then, level set techniques have been widely used for problems of imaging in different applications with a lot of success (see [22] for a recent overview on this issue). The main advantage of a level set technique is that it frees us from topological restrictions when reconstructing an a-priori arbitrary number of fluorophores. To our knowledge, this technique has not been carried out previously for imaging the fluorescence decay rate map within biological tissue.

The paper is organized as follows. In Section 2 we describe the diffusion equations modelling the transport of excitation and emission of light. In Section 3 we present the proposed algorithm and we briefly explain the level set technique. In Section 4 we give some details of the reconstruction algorithm. In Section 5 we show our numerical expriments. Section 6 contains our conclusions. Finally, in Appendix I and II we derive the gradient directions of the cost functional with respect to the fluorophore absorption coefficient and with respect to the fluorescence lifetime value, respectively. The cost functional describes the mistmatch between the ’true’ data and the data predicted by the model.

2. Light propagation model

Modeling fluorescence in tissues must account for the following stages: (i) propagation of excitation light from the tissue’s surface into its interior, (ii) absorption by fluorophores, (iii) conversion to fluorescence, (iv) emission of the fluorescent light from the fluorophores, and (v) propagation of that light back up to the tissue surface (we neglect here re-emission of the absorbed light). Besides, we assume that the absorption and emission spectra of the fluorescent molecules do not overlap. Hence, the excitation and emission processes take place at different wavelengths denoted by λ_x and λ_m > λ_x, respectively.

The diffusion equation has been widely used in diffuse optical tomography for imaging several-centimiter-thick turbid tissue by several authors (see [3], and references therein) and, in particular, for solving the inverse fluorescent source problem in optical molecular imaging [10, 26, 27, 28, 29, 30, 19]. It is a good approximation to the radiative transport equation, that describes how light propagates through an absorbing scattering medium, when the medium is optically thick and the fluorescent sources are not close to the detectors on the tissue surface. An alternative approximation in tissue optics is provided by the Fokker-Planck approximation that takes the forward-peaked scattering into account analitically by replacing the integral operator in the transport equation by a simple differential operator [31, 32, 33].

The diffusion equations modelling the transport of excitation and emission light are

\frac{1}{v} \frac{\partial U_{x} (r, t)}{\partial t} - D_{x} \nabla^{2} U_{x} (r, t) + (μ_{a}^{x} + μ_{a}^{x \to m}) U_{x} (r, t) = δ (r - r_{s}) δ (t - t_{0})

\frac{1}{v} \frac{\partial U_{m} (r, t)}{\partial t} - D_{m} \nabla^{2} U_{m} (r, t) + μ_{a}^{m} U_{m} (r, t) = S_{x \to m} (r, t) .

In these equations, U_x and U_m are the average diffuse intensities corresponding to the excitation and emission light propagating at wavelengths λ_x and λ_m, respectively, with velocity v. They depend on position r ∈ Ω ⊂ ℝⁿ (n = 2,3) and time t ∈ [0, T]. In this paper, we consider a two-dimensional geometry. In a two-dimensional domain, the diffusion coefficients at the excitation and emission wavelengths are

D_{x} = \frac{1}{2 [(1 - g) μ_{s}^{x} + (μ_{a}^{x} + μ_{a}^{x \to m})]}, D_{m} = \frac{1}{2 [(1 - g) μ_{s}^{m} + μ_{a}^{m}]} .

In Eq. (3), μ_a^x and μ_s^x, and μ_a^m and μ_s^m, denote the absorption and scattering coefficients at the excitation and emission wavelengths, respectively. The absorption by fluorophores in Eq. (1) is given by the fluorophore absorption coefficient μ_a^x→m. In Eq. (3), g denotes the anisotropy factor representing the average of forward (g > 0) or backward (g < 0) scattering.

In this paper, we consider that the tissue is illuminated at time t ₀ by a short laser pulse incident upon a point r _s on its boundary. This is modelled by the right hand side of Eq. (1). The generation of emission light at position r and time t defines the source term

S_{x \to m} (r, t) = \frac{η}{τ (r)} μ_{a}^{x \to m} (r) \int_{0}^{t} e^{- \frac{(t - t')}{τ (r)}} U_{x} (r, t') d t'

in Eq. (2), with τ(r) denoting the fluorescence lifetime distribution. It results from the convolution over time of an exponential term describing the decay of emission, and the average excited intensity at time t. The stimulated emission (4) is directly proportional to the quantum efficiency η of the fluorophore, which quantifies the conversion to fluorescence, and the fluorophore absorption coefficient μ_a^x→m(r). The product of the quantum efficiency and the fluorophore absorption is called the fluorescence yield.

Equations (1)-(2) are solved with appropriate boundary and initial conditions. To model zero total inward directed flux we use the Robin-type boundary condition

U_{x, m} (r, t) + \frac{π}{2} D_{x, m} \frac{\partial U_{x, m} (r, t)}{\partial \hat{n}} = 0, r \in \partial Ω,

where n̂ is the outward unit normal to the tissue surface ∂Ω. In addition, we assume that no energy is present in the medium at time t = 0, so the initial conditions are

U_{x, m} (r, 0) = 0 r \in Ω .

Equations (1)-(2) subject to conditions (5)-(6) comprise a complete description of the direct problem of fluorescent diffuse optical imaging. We solve this set of equations using a second order finite difference Crank-Nicolson scheme for the spatial operators in a mesh of 80 × 80 pixels. To evaluate the integral involved in the source term S_x→m(r,t) we use a classical trapezoidal rule.

3. Reconstructing the fluorescent source

Our objective here is to reconstruct the fluorophore distribution μ_a^x→m(r) and the fluorescence lifetime distribution τ(r) of the fluorescent source (4) using the physically measured data

{\tilde{G}}_{jl} = - D \frac{\partial {\tilde{U}}_{m}}{\partial \hat{n}},

taken from the tissue surface at positions r _l, l = 1,…, l, for the exciting source in Eq. (1) with r _s = r _j. In Eq. (7), Ũ_m represents the average intensity correponding to the ’physically correct’ parameters $\tilde{μ}$ _a^x→m and τ̃. We will denote the set of real measurements {G̃_jl}_l=1,…,l, for a fixed source j, by the vector G̃_j.

The goal of reconstructing the fluorescence source distribution requires the detection of well-defined boundaries between the fluorphores and the homogeneous backgound. Because traditional reconstructions techniques that seek to recover the fluorescence intensity at each pixel/voxel of the image often lead to blurring of these boundaries, we propose a shape-based approach that use level-set techniques. Shape-based reconstruction techniques preserve the sharp edges and enhance the contrast in the image.

In the level-set approach of shape reconstruction, the unknown shape S of the fluorophore distribution is implicitly represented by a sufficiently smooth ’level-set function’ ϕ as follows:

μ_{a}^{x \to m} (r) = {\begin{array}{c} μ_{a, in}^{x \to m} inside S, where ϕ (r) \leq 0, \\ 0 outside S, where ϕ (r) > 0 . \end{array}

Here, μ_a,in^x→m denotes the fluorophore absorption coefficient which we will assume to be constant inside the fluorophore (in the discussion below we have set the quantum efficiency η to 1). The boundary of the fluorophore, δS, consists of all points where ϕ(r) = 0, and the fluorophore distribution is different than zero at all points where ϕ(r) ≤ 0. We will indicate the dependence of the parameter ^μ_ax→m on the level-set function ϕ, by μ_a^x→m[ϕ]. The main advantage of this implicit representation of the unknown shape by a level-set function is its capability of automatically splitting and merging shapes during the reconstruction.

With the above definition we can formulate the shape reconstruction problem using level set as follows. Assuming that the optical properties of the tissue and the fluorophore are known, and given the measured data (7) on the tissue surface, find a level set function ϕ̂ such that μ_a^x→m[ϕ̂] matches the data.

To solve the shape reconstruction problem we will follow a time evolution approach [23, 34]. We introduce an artificial time ξ (related to the iterative step of the reconstruction process), and use the evolution law

\frac{dϕ}{dξ} = f (r; ξ)

for the unknown level set function ϕ. The purpose of the evolution law (9) is to reduce, and eventually minimize, the least-squares cost functional

퓙 (ϕ (ξ)) = \frac{1}{2} {∥ ℛ (μ_{a}^{x \to m} [ϕ (ξ)]) ∥}_{2}^{2},

where

ℛ (ϕ (ξ)) = G_{j} (ξ) - {\tilde{G}}_{j},

denotes the mismatch between the ’physically’ measured data G̃_j, and the predicted data

G_{j} = - D \frac{\partial U_{m}}{\partial \hat{n}}

using our model with parameter distribution given by μ_a^x→m [ϕ(ξ)].

The unknown forcing term f(r, ξ) in Eq. (9) will be chosen such that the cost functional (10) decreases for a sufficiently small time step of the algorithm. Note that in this formulation the least-squares cost functional (10) depends on the artificial time ξ. Formally differentiating the least squares cost functional 𝓙(μ_a^x→m[ϕ(ξ)]) with respect to the time ξ and applying the chain rule yields

\frac{d 퓙}{d ξ} = \frac{\partial 퓙}{\partial μ_{a}^{x \to m}} \frac{\partial μ_{a}^{x \to m}}{\partial ϕ} \frac{d ϕ}{d ξ} = < {grad}_{μ} 퓙, \frac{\partial μ_{a}^{x \to m}}{\partial ϕ} \frac{d ϕ}{d ξ} >_{P} =

where <>_P represents the canonical inner product in the parameter space P. Using this equation, we can select a descent direction for the cost functional by choosing

f (r; ξ) = - {grad}_{μ} 퓙 (r; ξ) for all r \in Ω .

Note that in order to solve Eq. (9), we have extended f to the whole domain (rigorously, f is only derived from Eq. (13) at those points where ϕ(r) = 0, i.e., at the boundary shape). The chosen ’extended velocity’ f(r;ξ) has the useful property that it allows for the creation of objects at any point in the domain, by lowering a positive-valued level set function until its values become negative. We refer to [22] for more details about other possible extensions of f(r;ξ).

An efficient way to compute the (pixel-based) gradient direction of 𝓙 in μ_a^x→m is to use the adjoint formulation. We give here the main result (see Appendix I for details).

The gradient direction of 𝓙 with respect to μ_a^x→m at each articial time ξ is given by

{grad}_{μ} 퓙 (r; ξ) = \frac{η}{τ (r)} \int_{0}^{T} dt W (r, t) \int_{0}^{t} d t' e^{- \frac{(t - t')}{τ (r)}} U_{x} (r, t'),

where U_x solves Eq. (1) with Eq. (8), and W solves the following adjoint equation

- \frac{1}{v} \frac{\partial W (r, t)}{\partial t} - D_{m} \nabla^{2} W (r, t) + μ_{a}^{m} W (r, t) = 0

W (r, t = T) = 0

W (r, t) = ℛ (μ_{a}^{x \to m} [ϕ]) where r \in \partial Ω .

Numerically discretizing Eq. (9) by a straightforward finite difference time-discretization with time-step Δξ ⁽ⁿ⁾ > 0 in step n, and interpreting ϕ ⁽ⁿ⁺¹⁾ = ϕ(ξ ⁽ⁿ⁾ + Δξ ⁽ⁿ⁾) and ϕ ⁽ⁿ⁾ = ϕ(ξ ⁽ⁿ⁾) yields the iteration rule

ϕ^{(n + 1)} = ϕ^{(n)} + Δ ξ^{(n)} f^{(n)} (r), ϕ^{(0)} = ϕ_{0},

where f ⁽ⁿ⁾(r) at time step ξ ⁽ⁿ⁾ is given by Eq. (14).

Note that, so far, we have assumed that we knew τ on those points where μ_a^x→m(r) ≠ 0. To simultaneously search for τ(r) and the fluorophore distribution μ_a^x→m(r) we apply the following update at each artificial time step ξ ⁽ⁿ⁾ (see Appendix II for details):

τ^{(n + 1)} (r) = τ^{(n)} (r) + Δ τ (r),

where

Δ τ (r) = - \frac{η μ_{a}^{x \to m} (r)}{{(τ^{(n)} (r))}^{2}} \int_{0}^{T} dt W (r, t) \int_{0}^{t} d t' (- 1 + \frac{t - t'}{τ^{(n)} (r)}) e^{- \frac{(t - t')}{τ^{(n)} (r)}} U_{x} (r, t'),

where U_x solves Eq. (1) with Eq. (8), and W solves the adjoint problem (16)-(18). Note that (i) since the supports of μ_a^x→m and τ coincides, the updates (21) are only different than zero within the fluorophores, and (ii) these updates are at almost no computational cost because the only require the knowledge of U_x and W which anyway are part of the calculation of the updates of the level set function ϕ.

In the case that the fluorophore distribution, defined by μ_a^x→m(r), is composed of an arbitrary collection of (connected or disconnected) disjoint compact subdomains S_p, p = 1,...,p, each having a different, but constant, τ_p in its interior, one can further regularize the inversion applying

Δ τ_{p} = \int_{S_{p}} d r Δ τ (r) .

4. Reconstruction algorithm

In this section we outline our algortihm to reconstruct the fluorophore and fluorescence lifetime distributions. The algorithm consists of two stages. During the first stage we apply an iterative pixel-based strategy to achieve successive improvements of an initial guess by minimizing the data misfit cost functional. During the second stage we assume the ’shape character’ of the fluorophores, and we apply the level-set technique.

I. Pixel-based stage.

Start with an initial fluorophore profile μ_a^x→m(r) = 0 at all points of the domain. Activate the first exciting source j = 1.
Solve the direct problem (1)-(2) subject to Eqs. (5)-(6) to compute the predicted data G_j. In the first step of the algorithm G_j = 0 since μ_a^x→m(r) = 0 in the whole domain.
Compute the residuals ℛ = G_j - G̃_j, and solve the adjoint problem (16)-(18).
Use expresions (15) and (21) to update the fluorophore and fluorescence lifetime distributions, respectively, at each pixel of the domain.
Check the stopping criterium. If it is not satisfied, activate the next source j = j + 1 and go to step (2). Our stopping criterium is that the cost functional becomes quasi-stationary.
If the stopping criterium is satisfied save the reconstructed profiles and go to the shape-based stage. A result of such a reconstruction at the end of this stage can be seen, for example, in the top right image of Fig. 1.

II. Shape-based stage.

Replace the pixel profiles by piecewise distributions. To this end, compute a threshold value of the fluorophore absorption coefficient μ_a^x→m(r) equal to 1/e times its maximum value over the whole tissue, and introduce a level set function ϕ ⁽ⁿ⁾(r), with n = 0, which is positive at the locations of lower fluorophore absorption value and negative at the locations of higher fluorophore absorption value. Assign to each fluorophore region the assumed known value μ_a^x→m of the fluorescence markers and the average values of τ(r) over each region. The result is a piecewise distribution defined by the introduced level set function ϕ ⁽⁰⁾(r), as shown in the center left image of Fig. 1.
Activate the first source j = 1 and solve Eqs. (1)-(2) subject to conditions (5)-(6) for the direct problem with the latest best guess μ_a^x→m[ϕ ⁽ⁿ⁾]. This gives rise to the predicted data G_j ⁽ⁿ⁾.
Compute the residuals ℛ = G_j ⁽ⁿ⁾ - G̃_j, and solve the adjoint problem (16)-(18). The gradient direction of 𝓙(μ_a^x→m(ϕ ⁽ⁿ⁾)) is given by Eq. (15).
Compute the velocity direction f ⁽ⁿ⁾(r) given by Eq. (14) and apply the update (19) to the level set function ϕ ⁽ⁿ⁾. Smooth f ⁽ⁿ⁾(r) by solving a heat equation [34]. The step size Δξ ⁽ⁿ⁾ in Eq. (19) depends on the activated source, so the variations provided by each one are similar.
Apply Eq. (20) with Eqs. (21)-(22) to update the fluorescence lifetime value in each fluorophore region.
Check the stopping criterium. If it is not satisfied, activate the next source j = j + 1 and go to step (2) with n = n + 1.

5. Numerical experiments

In the numerical experiments presented in this section we consider a 2D tomographic configuration. The experimental set-up consists of a 5 cm × 5 cm. square with mean background optical parameters of scattering μ¯_s^x = μ¯_a^s = 10cm ^-1, absorption μ¯_a^x = μ¯_a^m = 0.1cm ^-1, and anisotropy factor g = 0.8. We consider four sources only, located at the centers of each edge of the domain.

Fig. 1. First numerical experiment: a multi-object lifetime reconstruction. On the first row: reference fluorescence lifetime profile (left) and pixel by pixel reconstruction at the end of the first stage of the algorithm (right). On the second row: initial fluorescence lifetime profile for the second stage (left) and reconstructed profile at the 13^th iteration. On the third row: final fluorescence lifetime reconstruction (left) and evolution of the cost during the two stages of the algorithm.

Download Full Size | PDF

They illuminate the tissue one after the other. The experiments were conducted for a variation of fluorescent inclusions of different sizes and positioned at different places. The inclusions all have the same fluorescence yield ημ_a^x→m equal to 0.1 cm ^-1, but distinct fluorescence lifetimes between 0.5 ns and 1 ns. To model tissue autofluorescence we add random endogenous sources at each pixel of the background using a uniform random function over the interval [0,0.001]. Besides, 1% white Gaussian noise was added directly to the data. All together, the autofluorescence and the Gaussian noise, yield to an estimate noise in the data between 3% and 4%. For each source, we record the temporal distribution of the outgoing flux at the edges using the diffusion forward model explained in section 2. We note that in these tissue samples, multiple scattering is very important resulting in severe image blurring.

The first numerical experiment deals with the reconstruction of a multi-object consisting of three disconnected fluorophores with distinct lifetimes, 0.5,0.75 and 1 ns, as shown in the top left image of Fig. 1 where we represent the reference fluorescence lifetime profile.

We start the inversion process by a classical pixel-based iterative approach. The reconstruction at the end of this first stage of the algorithm is shown in the upper right image of Fig. 1. It is apparent that at the end of this stage it is very difficult (or even impossible) to accurately estimate the presence of the three distinct fluorophores. We note, though, that other pixel-based schemes might be optimized (for instance, using more than the 4 sources in our experiments) yielding possibly slightly better results than those of our straightforward pixel-based implementation. Nevertheless, the ’oversmoothing effect’ of regularized pixel-based schemes will still apply making difficult to extract the important key characteristics of these images.

Fig. 2. Second numerical experiment: a T-shape fluorescence lifetime profile. Form left to right: reference fluorescence lifetime profile, final reconstruction at the end of the first stage (pixel by pixel reconstruction), and final fluorescence lifetime reconstruction at the end of our algorithm.

Download Full Size | PDF

To continue with a shape-based approach we first replace the reconstructed pixel by pixel profile achieved during the previous stage by the piecewise distribution shown in the center left image, as it was explained in Section 4. Once we have obtained this initial guess of the shapes of the fluorophores, we proceed with our shape-based reconstruction. The reconstructed fluorescence lifetime profile at the 13th iteration is shown in the center right image. Observe that one object is breaking into two pieces. Each piece evolves during the next iterations changing its position, shape and fluorescence lifetime value. At the end of the process we obtain a very good estimate of the key characteristics of the three different fluorophores, as it can be seen in the bottom left image of Fig. 1.

It is interesting to compare the quality of our final reconstruction with the ’classical’ pixel-based reconstruction (using our own implementation) presented in the top right image. The shape-based approach provides a much clearer identification of the different fluorophores and a more accurate estimation of their shapes and fluorescence lifetime values. This is so because the shape-based approach incorporates important prior information regarding the finite support of the fluorophores, which cannot easily be incorporated in a classical pixel-based scheme.

In the bottom right image of Fig. 1 we show the evolution of the cost during the two stages of the algorithm. At each stage, the cost is mainly monotonically decreasing until after a certain iteration becomes quasi-stationary.

In our second numerical experiment we want to recover a complex-shape fluorophore with τ = 0.75 ns, as shown in the left image of Fig. 2. As we can see, the reconstruction at the end of the first stage using the pixel-based approach, shown in the center image, gives rise to a rough estimate of the position of the fluorophore. However, because the diffusion of light results in the loss of imaging resolution we can not identify its shape. The shape-based algorithm applied to the next iterations provides a much clearer identification of the fluorophore key characteristics such as its shape and fluorescence lifetime value, as shown in the right image.

Finally, we investigate a situation in which two small fluorophores are very close to each other, as shown in the left image of Fig. 3. Their fluorescence lifetime is τ = 0.75 ns. The center image shows the pixel by pixel reconstruction at the end of the first stage. The image is too blurred and the two targets are not resolved. However, the shape-based reconstruction achieved at the end of our algorithm (see the right image) recovers the positions and shapes of the two fluorophores acurrately. The retrieved fluorescence lifetime value is also very good.

Fig. 3. Third numerical experiment: testing the spatial resolution. Form left to right: reference fluorescence lifetime profile, final reconstruction at the end of the first stage (pixel by pixel reconstruction) and final fluorescence lifetime reconstruction at the end of the algorithm.

Download Full Size | PDF

6. Conclusions

We have presented a novel shape-based algorithm for fluorescence lifetime tomography. To this end, we have also derived the pixel-based gradients using an adjoint formulation. For the shape reconstruction problem we have employed level-set techniques that allow for an implicit representation of the fluorophore shapes and frees us from topological restrictions during the reconstruction process. Our shape-based approach preserves the sharp edges of the fluorophores thereby enhancing the contrast of the images. We believe that the results reported in this paper are important for the design of lifetime-based optical molecular imaging systems.

Appendix I

We derive here equations (15)-(18).

To compute a descent direction of the cost (10) when μ_a^x→m varies, subject to conditions (1)-(2) and (5)-(6), we build the Lagrange function defined by

ℒ (μ_{a}^{x \to m}; Z, W) = 퓙 (μ_{a}^{x \to m}) -

and we apply the adjoint method. Obviously, when U_x and U_m satisfy Eqs. (1) and (2), A = B = 0 and ℒ = 𝓙 for all Z and W, so we have incorporated conditions (1)-(2) into one function. We now take the derivative of Eq. (23)

δ ℒ = \frac{\partial ℒ}{\partial μ_{a}^{x \to m}} δ μ_{a}^{x \to m} - \frac{\partial ℒ}{\partial U_{x}} δ U_{x} - \frac{\partial ℒ}{\partial U_{m}} δ U_{m} .

Note that if the two last terms are zero then δℒ = δ𝓙. Since, so far, Z and W are arbitray functions, we can choose them so that these two last terms vanish, i.e., we impose that

\frac{\partial ℒ}{\partial U_{x}} δ U_{x} = 0,

\frac{\partial ℒ}{\partial U_{m}} δ U_{m} = 0 .

From these conditions we derive the adjoint equations.

Let us now find the first ajoint equation. To first order we can write that

\frac{\partial ℒ}{\partial U_{x}} δ U_{x} = ℒ (μ_{a}^{x \to m}, U_{x} + δ U_{x}, U_{m}) - ℒ (μ_{a}^{x \to m}, U_{x}, U_{m}) =

Integrating by parts, appying the divergence theorem, and using Eqs. (5) and (6) we obtain

\frac{\partial ℒ}{\partial U_{x}} δ U_{x} = \int_{0}^{T} dt \int_{Ω} d r δ U_{x} {- \frac{1}{v} \frac{\partial Z}{\partial t} - D_{x} \nabla^{2} Z + (μ_{a}^{x} + μ_{a}^{x \to m}) Z},

where we have imposed that

Z (r, t = T) = 0 in Ω and, Z (r, t) = 0 on \partial Ω .

Since $\frac{\partial ℒ}{\partial U_{x}} δ U_{x}$ have to be zero for all δU_x, we arrive to an equation for Z. It reads

- \frac{1}{v} \frac{\partial Z}{\partial t} - D_{x} \nabla^{2} Z + (μ_{a}^{x} + μ_{a}^{x \to m}) Z = 0 .

The minus sign in front of the time derivative in Eq. (30) means backwards time integration. The solution to this equation with ’initial’ and boundary conditions given by Eq. (29), vanishes in all the domain, i.e.,

Z (r, t) = 0 for all r \in Ω, and t \in [0, T] .

Similarly, to determine W we impose Eq. (26). To first order we obtain

\frac{\partial ℒ}{\partial U_{m}} δ U_{m} = ℒ (μ_{a}^{x \to m}, U_{x}, U_{m} + δ U_{m}) - ℒ (μ_{a}^{x \to m}, U_{x}, U_{m}) =

and

\frac{\partial ℒ}{\partial U_{m}} δ U_{m} = \int_{0}^{T} dt \int_{Ω} d r δ U_{m} {- \frac{1}{v} \frac{\partial W}{\partial t} - D_{m} \nabla^{2} W + μ_{a}^{m} W},

where Eqs. (5) and (6) have been used, and where we have imposed that

W (r, t = T) = 0 in Ω and, W (r, t) = ℛ (μ_{a}^{x \to m}) on \partial Ω .

From Eq. (33) we get

- \frac{1}{v} \frac{\partial W}{\partial t} - D_{m} \nabla^{2} W + μ_{a}^{m} W = 0,

as we wanted.

Therefore, if we select Z and W so that they satisty Eqs. (30) and (35), supplemented with conditions (29) and (34), we have that

δ 퓙 = δ ℒ = \frac{\partial ℒ}{\partial μ_{a}^{x \to m}} δ μ_{a}^{x \to m} = < {grad}_{μ} ℒ, δ μ_{a}^{x \to m} >_{P},

from where

δ 퓙 = < {{grad}_{μ} 퓙, δ μ_{a}^{x \to m} >}_{P} = \int_{Ω} d r \int_{0}^{T} dt W \frac{\partial S}{\partial μ_{a}^{x \to m}} δ μ_{a}^{x \to m},

where < >_P represents the inner product in the parameter space. To obtain Eq. (37), we have taken the derivative of Eq. (23) with respect to μ_a^x→m and used the fact that Z vanishes in all the domain. Note that D_x and S(r,t) depend on μ_a^x→m. Finally, taking the partial derivative of the source S(r,t) with respect to μ_a^x→m in Eq. (37), we get that the gradient direction of 𝓙 with respect to μ_a^x→m is given by

{grad}_{μ} 퓙 (r) = \frac{η}{τ (r)} \int_{0}^{T} dt W (r, t) \int_{0}^{t} dt' e^{- \frac{(t - t')}{τ (r)}} U_{x} (r, t') .

With this we have proven equations (15)-(18).

Appendix II

In order to find a correction Δτ for the fluorescence lifetime distribution τ, we linearize the (nonlinear) residual operator (11), neglect terms of order O(∥Δτ∥²), and solve ℛ[τ+Δτ] = 0. This amounts to solving

ℛ' [τ] Δ τ = - ℛ [τ],

where ℛ́ represents the derivative of ℛ at τ. The minimum-norm solution of this equation can be written as

Δ τ = - R' {[τ]}^{*} {(R' [τ] R' {[τ]}^{*})}^{- 1} R [τ],

where Ŕ[τ]^* is the adjoint operator of Ŕ[τ]. Since in our application the operator C = (Ŕ[τ]Ŕ[τ]^*)^-1 is very expensive to calculate we approximate it by the identity operator I (note that C acts as a ’filter’ from the parameter space to the parameter space). Using this approximation, we just have to calculate

Δ τ = - R' {[τ]}^{*} R [τ] .

Note that Ŕ[τ]^* R[τ] coincides with the gradient direction of 𝓙 with respect to τ. Indeed,

𝓙 [τ + Δ τ] = 𝓙 [τ] + \underset{δ 𝓙}{\underset{︸}{< {grad}_{τ} 𝓙, Δ τ >_{P}}} + O ({∥ Δ τ ∥}_{P}^{2}) =

From Eq. (42), we can extract the gradient and write

{grad}_{τ} 퓙 = R' {[τ]}^{*} R [τ] .

Following the same steps as in Appendix I we readily find that

{grad}_{τ} 𝓙 (r) = \int_{0}^{T} dt W (r, t) \frac{\partial S}{\partial τ} (r, t),

With this, and using Eqs. (43) and (41), we have proven Eq. (21).

References and links

1. V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng. 8, 1–33 (2006). [CrossRef] [PubMed]

2. O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. 14, 1107–1130 (1998). [CrossRef]

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

4. M. Schweiger, S. R. Arridge, O. Dorn, A. Zacharopoulos, and V. Kolehmainen, “Reconstructing absorption and diffusion shape profiles in optical tomography using a level set technique,” Opt. Lett. 31, 471–473 (2006). [CrossRef] [PubMed]

5. P. Gonzalez-Rodriguez, A. D. Kim, and M. Moscoso, “Reconstructing a thin absorbing obstacle in a half space of tissue,” J. Opt. Soc. Am. A 24, 3456–3466 (2007). [CrossRef]

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

7. R. B. Schulz, J. Ripoll, and V. Ntziachristos, 1Experimental fluorescence tomography of tissues with noncontact measurements,” IEEE Trans. Med. Imaging 23, 492–500 (2004). [CrossRef] [PubMed]

8. A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. 202, 323–345 (2005). [CrossRef]

9. A. D. Kim and M. Moscoso, “Radiative transport theory for optical molecular imaging,” Inverse Probl. 22, 23–42 (2006). [CrossRef]

10. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21, 158–160 (1996). [CrossRef] [PubMed]

11. E. Shives, Y. Xu, and H. Jiang, “Fluorescence lifetime tomography of turbid media based on an oxygen-sensitive dye,” Opt. Express 10, 1557–1562 (2002). [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. A. Godavarty, E. M. Sevick-Muraca, and M. J. Eppstein, “Three-dimensional fluorescence lifetime tomography,” Med. Phys. 32, 992–1000 (2005). [CrossRef] [PubMed]

14. B. B. Das, F. Liu, and R. R. Alfano, “Time-resolved fluorescence and photon migration studies in biomedical and model random media,” Rep. Prog. Phys. 60, 227 (1997). [CrossRef]

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

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

17. S. Bloch, F. Lesage, L. McIntosh, A. Gandjbakhche, 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]

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

19. A. T. N. Kumar, S. B. Raymond, G. Boverman, D. A. Boas, and B. J. Bacskai, “Time resolved fluorescence tomography of turbid media based on lifetime contrast,” Opt. Express 14, 12255–12270 (2006). [CrossRef] [PubMed]

20. A. T. N. Kumar, S. B. Raymond, B. J. Bacskai, and D. A. Boas, “Comparison of frequency-domain and time-domain fluorescence lifetime tomography,” Opt. Lett. 33, 470–472 (2008). [CrossRef] [PubMed]

21. F. Santosa, “A Level set approach for inverse problems involving obstacles,” ESAIM Control, Optimization and Calculus of Variations 1, 17–33 (1996). [CrossRef]

22. O. Dorn and D. Lesselier, “Level set methods for inverse scattering,” Inverse Probl. 22, R67–R131 (2006). [CrossRef]

23. N. Irishina, M. Moscoso, and O. Dorn, “Microwave imaging for early breast cancer detection using a shape-based strategy,” IEEE Trans. Biomed. Eng. 56, 1143–1153 (2009). [CrossRef] [PubMed]

24. S. Osher and J. A. Sethian, “Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations,” J. Comput. Phys. 79, 12–49 (1988). [CrossRef]

25. F. Santosa, “A level set approach for inverse problems involving obstacles,” ESAIM Control, Optim. Calculus Variations 1, 17–33 (1996). [CrossRef]

26. D. Y. Paithankar, A. U. Chen, B. W. Pogue, M. S. Patterson, and E. M. Sevick-Muraca, “Imaging of fluorescent yield and lifetime from multiply scattered light reemitted from random media,” Appl. Opt. 36, 2260–2272 (1997). [CrossRef] [PubMed]

27. J. Chang, H. L. Graber, and R.L. Barbour, “Imaging of fluorescence in highly scattering media,” IEEE. Trans. Biomed. Eng. 44, 810 (1997). [CrossRef] [PubMed]

28. H. Jiang, “Frequency-Domain Fluorescent Diffusion Tomography: A Finite-Element-Based Algorithm and Simulations,” Appl. Opt. 37, 5337–5343 (1998). [CrossRef]

29. V. Ntziachristos and R. Weissleder, “Charge-coupled-device based scanner for tomography of fluorescent near-infrared probes in turbid media,” Med. Phys. 29, 803 (2002). [CrossRef] [PubMed]

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

31. A. D. Kim and J. B. Keller, “Light propagation in biological tissue,” J. Opt. Soc. Am. A 20, 92–98 (2003). [CrossRef]

32. A. D. Kim and M. Moscoso, “Beam propagation in sharply peaked forward scattering media,” J. Opt. Soc. Am. A 21, 797–803 (2004). [CrossRef]

33. A. D. Kim and P. Tranquilli, “Numerical solution of a boundary value problem for the Fokker-Planck equation with variable coefficients,” J. Quant. Spectrosc. Radiat. Transfer 109, 727–740 (2007).

34. D. Álvarez, O. Dorn, N. Irishina, and M. Moscoso, “Crack reconstructions using a level-set strategy,” J. Comput. Phys. (to be published).

Fluorescence lifetime imaging from time resolved measurements using a shape-based approach

Abstract

1. Introduction

2. Light propagation model

3. Reconstructing the fluorescent source

4. Reconstruction algorithm

I. Pixel-based stage.

II. Shape-based stage.

5. Numerical experiments

6. Conclusions

Appendix I

Appendix II

References and links

Cited By

Figures (3)

Equations (50)

Optics Express