Cramer-Rao analysis of steady-state and time-domain fluorescence diffuse optical imaging.

Using a Cramer-Rao analysis, we study the theoretical performances of a time and spatially resolved fDOT imaging system for jointly estimating the position and the concentration of a point-wide fluorescent volume in a diffusive sample. We show that the fluorescence lifetime is a critical parameter for the precision of the technique. A time resolved fDOT system that does not use spatial information is also considered. In certain cases, a simple steady-state configuration may be as efficient as this time resolved fDOT system.

1. Introduction

Fluorescence diffuse optical tomography (fDOT) has become a valuable tool for in vivo molecular imaging, and various fDOT systems have be proposed during the last two decades [1-3]. In this paper, we focus more specifically on two common modalities that can be identified in the literature, namely the continuous wave (or steady-state) and the time domain fDOT systems. In a continuous wave fDOT setup (CW-fDOT), one illuminates the sample with one or several continuous excitation sources and the emitted fluorescent light is detected on an array of detectors that is usually a CCD sensor [1, 4–6]. In a time domain fDOT system (TD-fDOT), the sample is illuminated with short light pulses and a time-resolved detection is performed using photon-counting methods with a single detector [7, 8] or with optimized arrangements of detectors [9, 11, 12, 14, 43] to get simultaneously temporal and spatial informations.

A key issue is the determination of the performances of a given setup, in terms of spatial resolution (location of fluorescent sources) and signal quantification (fluorophore concentration). The precision in the determination of these two parameters from tomographic measurements, using inverse reconstruction methods [10], depends on the sensitivity of the signal to these parameters and on the noise level [13]. TD-fDOT has received increasing interest, due to its (often assumed) ability to overtake CW-fDOT. Indeed, improvements have been reported in specific situations [9, 14, 15]. Yet, a systematic study of the performances of these imaging systems, and, in particular, of the influence of the lifetime of the fluorescence on the precision of the technique is still lacking.

In this paper, we present a rigorous methodology to calculate the precision limit of different imaging modalities. We first develop the analysis of a theoretical fDOT system based on ultrafast CCD sensor, allowing spatially and temporally resolved imaging; This analysis sheds light on the important role of the fluorophore lifetime on the depth localisation precision of the TD system. Then, two simplified systems are considered: an integrated time domain fDOT (ITD-fDOT) setup with no spatial information, and a CW-fDOT with no temporal information. These two “degraded” systems use the fluorescent signal in very different ways and comparing their ability to localize the depth of the fluorescent inclusion is of interest.

2. Light propagation models and fluctuations in the counting signal

Let us consider a fluorescent point-like volume located at r and embedded in a 9mm thick homogeneous diffusive slab (see Fig. 1). The determination of the location and intensity of a luminescent point source in such a diffusive medium is a standard problem in molecular imaging [16, 17].

Fig. 1

Schematic description of the reflection fDOT setup. L is the depth of the point-like fluorescent volume. The dots below the upper surface represent the 7×7 excitation laser sources. The detection is performed on a 32×32 pixel detector (CW-fDOT) or on a single detector with same area (TD-fDOT). The field of view is 25mm×25mm.

2.1. Tomographic configurations in reflection geometry

We consider a simple and theoretical contact optical tomography set-up in epi-illumination configuration and a reflection geometry as recently used in [18]. However, the methodology presented in this paper is general and can be adapted to deal with other instrumental configurations —e.g., setups in transmission geometry or with non contact detection devices (provided that one accounts for the free-space propagation of the light [19]). In this paper, a set of 7 × 7 excitation sources are placed on the upper side (z = 0 mm) of the slab over a 7 mm×7 mm surface. They deliver either a continuous light beam (CW-fDOT) or a Dirac delta-function light pulse (TD-fDOT). These sources are modelled by point-like isotropic sources placed inside the medium at a depth corresponding to one transport mean free path ℓ* = [μ(1 − g)]−1 [14, 20] with μ the scattering coefficient and g the anisotropy factor. The fluorescence signal is collected using 32 × 32 joined detectors (pixels) with area ΔΣ placed on the same side as the illumination and covering a 25 mm × 25 mm surface. In our TD-fDOT experiment, the spatially resolved fluorescence signal is also time resolved, the emitted photons being collected within time bins with duration Δ. We consider also a simplified spatially integrated time-domain experiment (ITD-fTOD) in which all the emitted photons impinging on the surface of the CCD camera are detected within the time bin Δ.

2.2. The expected counting signal

The Cramer-Rao analysis presented in this paper requires a light propagation model in the slab. This subsection describes this model that rests upon the diffusion approximation [27] and the method of images [28]. We emphasize that such a model may be inaccurate, especially if the energy density is computed near a boundary or close to a source, see [30, 31]. In such cases, the methodology presented in this paper may still be implanted provided that an appropriate light transport model is adopted, see for instance [13].

Let us consider first a single pointwise excitation source located at the position s and emitting a temporal pulse of energy ℰ0 at t = 0. Given an arbitrary density distribution of fluorescent sources n(r), the energy density on the detection surface, at the position r, reads [21, 24, 25, 29] where ★ is the temporal convolution operator, and η and σ is the quantum efficiency and the absorption cross-section, respectively. The time response function R(t; τ) describes dynamics of photon emission by a fluorophore. Assuming single-exponential decay (which is the case for organic fluorophores widely used in biomedical optics), one has R(t; τ) = exp(t/τ)/τ, where τ is the lifetime of the excited state. The expression of the time-dependent fluorescence intensity given by Eq. (1) has been widely used to model fDOT experiments, either in the time-domain [23, 24, 29] and in the frequency domain [21, 25]. The photon transport, at both the excitation and emission wavelengths, is described in the diffusion approximation [27] by the diffusive Green function for the slab geometry G(r, t) computed using the method of images [28]. We assume that the characteristics of the model transport are the same at the emission and excitation wavelengths.

For a monochromatic excitation located in s, the average (or expected) photon count detected in a temporal bin centered in t and in the pixel centered in r is where 〈...〉 denotes the ensemble average and with Π(t) and Π(r) the support of the temporal bin and of the pixel, respectively. h and ν are the Plank constant and the frequency of the collected photons, D = cℓ*/3 is the diffusion coefficient with c the energy velocity in the medium, and z* is the extrapolation length in the diffusion approximation [27]. Provided that the duration of the bin Δ and the area of the pixel ΔΣ are small enough, the counting signal induced by the pointwise fluorescent volume buried in the slab is on average where 𝒜 = cησD/(hνz), and with r the position of the fluorophores and ρ the fluorophore concentration defined by n(r) = ρ δ (r – r) with δ(r) the Dirac delta function in space. Finally, let us introduce and the number of detected photons on the k-th pixel of a CCD camera in the CW-fDOT situation, or the number of detected photons in the p-th temporal bin in the ITD-fDOT case. The expected number of detected photon induced by the excitation in s is given by and with T and 𝒫0 the exposition time and the power of the laser illumination for the steady-state experiment.

In practice, the choice of the incident energy ℰ0 (in the TD or ITD cases), or of both the acquisition time T and the incident power 𝒫0 (in the CW case) are set by the experimental conditions. The integration time is limited by the drift of the sample and the photobleaching of the fluorophores while the power of the excitation pulse is limited by the available sources and the possible damage caused to the sample. However, for a fair comparison of the considered techniques, we assumed that for a given position of the fluorescent volume, the total number of detected photons is the same whatever the fDOT configurations. This condition ensures that all the experiments exhibit an equivalent signal to noise ratio. Note that, without this constraint, one could increase the integration time of the detector in the CW set-up and obtain infinite signal to noise ratio (and infinite precision). We shall also consider configurations that have been designed so that the truncation and the discretization of the space or time varying signals yield negligible loss of information. As a result, for a given position of the fluorescent volume and an excitation source, the three setups yield (in average) the same total amount of detected photons, i.e. we get with K and P the number of pixels in the field-of-view of the camera and the total number of time bins, respectively.

For the sake of simplicity, we restrict our localization problem to a one-dimensional problem such that the depth of the source (denoted L) is the only position parameter that is unknown. Hence, in the sequel of the paper, both the position L and the concentration ρ are unknown parameters of interest, whereas the time response R(t;τ) is assumed to be known —lifetime τ included.

As an illustrative example, Fig. 2 shows the averaged fluorescence signals (normalized by the excitation energy ℰ0) computed in the CW and ITD configurations for a single excitation source located in (0, 0, ℓ) and for two different locations of the fluorescent source L = 3 mm and L = 4.5 mm. In both cases, the scattering and absorption coefficients are μ = 36.7 mm−1 and μ = 0.029 mm−1, with an anisotropy factor g = 0.9, which is a coarse description of a mouse brain where the optical properties have been homogenized [13]. The chosen characteristics of fluorescent probes are standard for fluorophores such as Cyanines 5 [26]: The quantum yield is η = 0.3, the absorption cross-section is σ = 10−16 cm−2 and the lifetime is τ = 1 ns. The solid lines correspond to a concentration ρ = 1 whereas the dashed lines were obtained with adjusted values of ρ so that the maxima of the expected signals coincide for both depths. For the CW setup, the plots are radial profiles of the averaged fluorescence signal detected on the CCD camera.

Fig. 2

Expected counting signal in CW-fDOT (a) or in ITD-fDOT (c) obtained for a single excitation source (i.e. N = 1) located in s = (0, 0, ℓ) and a fluorescent volume located L = 3 mm in the slab —these signals are normalized by the incident excitation energy ℰ0. The shaded area corresponds to the standard deviation due to photo noise. (b,d): Same as (a) and (c) for L = 4.5 mm (blue solid curve). The red dashed line is the expected signal for L = 3 mm, but with a concentration adjusted so that the maxima of the signals coincide for both depths. For the ITD configuration (c,d), the considered fluorophore lifetime is τ = 1 ns.

As shown in Fig. 2, the parameters ρ and L affect the signal in a very different way: the source depth L modifies the shape of the curve while the source concentration ρ acts as a multiplicative factor. For the time resolved experiments, the expected signals Eq. (3) and Eq. (4) are expressed as a multiple convolution in time involving the time response of a fluorophore. For this reason, the ability to recover both the position and the concentration of the fluorescent source is certainly dependant on the lifetime τ. The next subsection illustrates the potential difficulty in estimating accurately both these parameters from data degraded by photon noise.

2.3. The noise fluctuations in the data-set

In the following, let be the data-set drawn from any of the instrumental setup considered this paper, with N the total number of measurements [32]. Similarly, is a generic average number of detected photon that corresponds to Eq. (3) in the TD case, or to Eq. (5) and Eq. (4) in the CW and ITD case, respectively.

To investigate the sensitivity of the considered imaging configurations, we assume the ideal case where the experiments are only limited by photon noise. In this case, the number of detected photons m obeys Poisson statistics. Assuming further statistical independence of the data-set, the probability law for the observed data-set for any instrumental setup that is considered in this paper reads

Assuming that the geometry and optical properties of the diffusive sample are perfectly known, the estimation could be done by fitting the time and/or space varying signals with a two-parameters model of the experiment. However, the noise spoiling the data may render this a priori simple task most difficult. In Fig. 2(b) and 2(c), we plotted the expected fluorescence signals in the CW and ITD configurations for L = 3 mm and L = 4.5 mm. One standard-deviation error-bars are also depicted by a shaded area [33]. Two different concentrations of fluorophores ρ are considered so that the maxima of both signals coincide —the concentration for L = 3 mm being much smaller than that taken for L = 4.5 mm.

We observe that, although the exact spatial or temporal dependence of the signals are different for both depths, the general shape remain similar and the signals at L = 3 mm and L = 4.5 mm are included inside the error bars. Hence, whatever the chosen imaging configuration, (CW or ITD), the signal of a weak fluorescent source at a given depth could be erroneously associated to that of a source with stronger intensity located deeper in the slab. This simple observation shows that the precision of the joint estimation of L and ρ can be very poor. Furthermore, for time resolved experiments (TD or ITD configurations), the time response of the fluorescent source R(t;τ) will “blur” the temporal counting signal, that in turn increases the ambiguity in depth of the fluorescent source. Figure 3 shows an illustration of this phenomenon: the counting signals from the depths L = 5 mm and L = 2.5 mm are hardly distinguishable if τ = 2 ns whereas τ = 0 ns allows a clear separation of the signals.

Fig. 3

(Blue solid curve) Averaged temporal signal in the ITD configuration for a single excitation source (i.e. N = 1) located in s = (0, 0, ℓ) and a fluorescent volume located L = 5 mm with τ = 0 ns (a) or τ = 2 ns (b) —these signals are normalized by the incident excitation energy ℰ0. The shaded area corresponds to the standard deviation due to photo noise. (Red dashed line) Averaged signal for a source in L = 2.5 mm with τ = 0 ns (a) or τ = 2 ns (b), but with a concentration adjusted so that the maxima of the averaged signals coincide for both depths.

Now these illustrative examples are not sufficient to draw quantitative precision limits concerning the imaging modes. One needs a formal quantitative tool that is able to measure the precision of a given imaging configuration. The Cramer-Rao bound (CRB) analysis provides a perfect rationale for this task.

3. Cramer-Rao analysis

Because of noise, different realizations of the same experiment yield different estimations for L and ρ. The depth and concentration estimates can be considered as random variables and one can define their standard deviations, δ and δ, which measure the precision of the experiment. Using the Cramer-Rao formalism [34], it is possible to calculate lower bounds for δ and δ that depend solely [35] on the model describing the behavior of the detected signals and the noise Eqs. ((3)–(6)).

3.1. The vectorial Cramer-Rao bound

Let us introduce the Fisher information matrix where . This matrix is readily deduced from Eqs. (6) and reads with ∂ = ∂/∂L. The (Cramer-Rao) bounds for δ and δ provide lower bounds for the variance in the unbiased estimation of L and ρ. Thus, we introduce the following precision limits [34] where CRB and CRB are defined via the diagonal elements of F−1: These expressions generalize the single-parameter result previously described in [13]. The coefficient [resp. ] represents the precision limit for the depth [resp. concentration] that is obtained when the intensity [resp location] of the source is known [13]. The actual precision limit 𝒫 is larger than , except when L and ρ affect the signal in a uncoupled way, which is a very particular case. We stress that Eq. (9) describes quantitatively the cross-coupling between the precision limits for the depth L and for the concentration of fluorophores ρ.

The value of 𝒫 and 𝒫 depend on the number of collected photons. Indeed, for photon noise limited imaging, theses precision limits decreases with the energy of the excitation source ℰ0. Furthermore, the depth precision limit 𝒫 decreases with the concentration ρ, the inverse being true for the concentration precision limit 𝒫 —i.e. increasing the concentration ρ deteriorates the bound for this parameter [36]. Let us note that the CRB depends also on the size of the elementary elements chosen for discretizing the signal (pixel size Δ or time bin Δ). As already explained in subsection 2.2, we wish to avoid considering these additional parameters in the discussion. Hence, in order to take advantage of all the information contained in the signals, we have chosen elements that are small enough for the sums in Eq. (7) to reach their asymptotic values. More specifically, we chose ΔΣ = 0.64 mm2 and Δ = 5 ps. Finally, let us underline that the CRB provides the best precision that the imaging experiment can achieves, but it does not provide any tool to retrieve it. Since the pioneering work of R.A. Fisher [37], the estimation from noisy data-set is usually performed by maximizing the likelihood function with respect to (ρ, L) ; in our context, the likelihood function is simply the expression (6) as a function of the unknown parameters (ρ,L). In the asymptotic limit of high counting rates, this estimation technique is unbiased and reaches the lower bounds defined by Eq. (9) [38].

3.2. The intrinsic depth precision limit

The depth precision limit for the various fDOT configurations presented in section 2 are now considered. Note that the analysis of the precision limit on the concentration yields similar results. Indeed, we introduce the following normalized depth precision limits that does not depend neither on the excitation energy ℰ0 nor on the concentration ρ of the fluorescent source where “•” stands for TD, ITD or CW. More specifically, the precision limits and depend on both the depth L and the lifetime τ, whereas is only a function of the depth L. In addition, we note that these three precision limits depend on the quantum yield η and on the absorption cross-section σ [39]. Let us remind that in our study the CW situation corresponds to the use of many detectors to resolve spatially a steady-state signal, while the ITD case corresponds to a time-resolved detection on a single extended detector.

For TD-fDOT in the reflection geometry depicted in Fig. 1, the precision limit is depicted in Fig. 4. These curves are computed using Eqs. (3), (7) and (9) for typical values of the fluorescence lifetime τ. Firstly, we note that the precision of the depth estimation decreases with the depth of the fluorescent volume (the precision limit increases). This is due to a decrease of the signal-to-noise ratio and of the sensitivity of the signal to a variation of the source position (both the spatial and temporal signal flatten out for deeper sources). Secondly, the precision on the depth estimation substantially depends on the fluorescence lifetime τ. The longer the lifetime, the poorer the precision limit, the worst precision been achieved for τ → ∞ which corresponds to the precision limit in the CW configuration, i.e. . Qualitatively, the fluorescence lifetime affects the TD signal in a similar way as the source depth, reducing the signal sensitivity to the source depth. Finally, we note that the precision limits and are almost identical for τ = 0: In this limit case, the CCD sensor in the considered TD-fDOT system does not help in the depth estimation of a single source. However, for typical lifetime values, spatially resolved detectors does bring useful information for the localization of the fluorescent source. In particular, for a given number of incident photons on the system, the CW and the ITD setups in reflection geometry achieve comparable precision limits, see Fig. 5. In this case, and for the reflection geometry, we observe that the precision in CW is comparable to that in the ITD technique. The precision limit can even be smaller in CW for a depth L that is smaller than the photon diffusion length within the fluorescent lifetime . Of course, the precision of the the CW situation could be improved by increasing the integration time of the CCD. Finally, we note that in the limit τ → ∞, the ITD experiment reduces to a CW configuration, but with a single extended detector that collects all the emitted photons (no spatial sampling). For the estimation of L and ρ, this is a worst case scenario that still achieves a finite precision limit since this fDOT system provides N = 7 × 7 independent measurements to retrieve the unknown pair of parameter (L, ρ).

Fig. 4

(Black solid lines) Normalized precision limit for the estimation of L with different values of the fluorescence lifetime τ in the TD-fDOT situation, see Sec. 2 for details. (Red dashed line with circles) Normalized precision limit for τ = 0 ns in the ITD-fDOT setup. (Blue dashed line with triangles) Normalized precision limit for the CW-fDOT setup.

Fig. 5

(Black solid lines) Normalized precision limit for different values of the fluorescence lifetime τ in the ITD-fDOT situation, see Sec. 2 for details. (Blue dashed line with triangles) Normalized precision limit for the CW-fDOT setup.

4. Conclusion

In summary, we have introduced a rigorous methodology to compare the theoretical performances of a single source localisation in an homogeneous medium for the time-resolved and/or spatially-resolved fDOT setups in the reflection geometry depicted in Fig 1. However, for the considered time-resolved (i.e. TD and ITD) configurations, this analysis reveals the critical role of the fluorophore lifetime for the localization capability of the setup. In particular, our results suggest that the best depth precision in fDOT imaging are expected with the smallest lifetime. Furthermore, the depth precision capability of a time-resolved (single-detector) ITD modality is not necessarily better than that of a spatially-resolved CW technique.

Finally, let us note that the fDOT systems considered in this paper are simple examples chosen to illustrate the potentiality of a Cramer-Rao analysis as a performance index for a given setup and a given task. We do not claim to compare CW and time resolved fDOT systems in general since we think that a universal answer is probably out of reach. However, the proposed methodology does not restrict to the fDOT systems evaluated here and it can be adapted to other fDOT systems like frequency domain systems [42, 45–47] or time-gating techniques [9, 14].