1. Introduction

Nowadays, the precise determination of the lifetime of fluorescent emitters has become essential for a wide range of applications. Indeed, enhancing the spontaneous rate of single emitters is a topical challenge in nanophotonics [1, 2], leading to the measurement of strongly reduced fluorescence lifetimes [3–8]. In biology, the contrast induced by lifetime variations is used to map different parameters on biological samples [9–11] such as the viscosity, the potential of hydrogen (pH) or the interaction between two emitters due to Förster resonance energy transfer (FRET). Whatever application one is interested in, the best precision that can be achieved depends on both the experimental conditions and the efficiency of the estimators. While the performances of lifetime estimators can numerically be studied using Monte Carlo experiments [12–16], finding the right experimental conditions requires a reliable benchmark to compare the performances of different experimental setups. Such a benchmark is set by Fisher information theory and the Cramér-Rao inequality, which gives the lower bound on the variance of unbiased estimators [17]. In other words, the Cramér-Rao inequality allows one to calculate the best precision that can be achieved in the estimation of one or several parameters, taking into account the various constraints induced by an experiment. In the recent years, this gauge became a standard used to assess the limit of localisation precision in the context of single-molecule microscopy [18–22]. Other recent applications of this formalism include the investigation of the dynamics of single molecules on the millisecond time scale [23] and the comparison of different imaging modalities in fluorescence diffuse optical tomography [24,25]. Among the different techniques that can be implemented for lifetime measurements [11,26,27], TCSPC is commonly used to exploit low level light signals with picosecond resolution [28]. Using the Cramér-Rao inequality, a relation between the estimation precision and the number of collected photons was obtained in 1992 for a simplified TCSPC model [29]. Here, we perform an extensive Cramér-Rao analysis to unravel the dependence of the lifetime estimation precision on experimentally relevant parameters, in the general case of a bi-exponential distribution along with any non-uniform background signal and considering the finite IRF of the setup. In addition, we provide an open-source code ( Code 1 [30]) which computes the Cramér-Rao bound for any set of experimental parameters. Within this framework, we thus provide a versatile tool which can be used for different purposes such as determining the shortest lifetime that can be probed with a TCSPC setup or achieving optimal contrast for FLIM-based applications.

2. Results and discussion

2.1. Excited-state lifetime distribution

Let us consider that the measured decay histogram follows a bi-exponential distribution with an additional contribution due to background noise. Usually this noise originates from dark counts due to the detection process, unfiltered excitation laser or luminescence of the substrate. Assuming that the noise follows a known probability density function (PDF) noted q_b(t), the total signal can be modeled by using a set of 5 unknown parameters noted θ, namely, the decay rate of each component (Γ₁ and Γ₂), the average number of detections for each component (N₁ and N₂), and the average number of detections due to background noise (N_b). Furthermore, we consider that the excitation laser has a finite pulse duration, and that the jitter of the detection system induces a loss of precision over the photon detection time. These two effects can be accounted for by measuring the IRF of the system, which is described by a PDF noted q_irf(t). Then, if the detected photon rate does not exceed the maximum counting speed of the detector, the expected number of events f_i detected in the i-th bin of the decay histogram (and associated with delays in-between t_i and t_i+1) reads

2.2. Cramér-Rao lower bound

In estimation theory, a well-known result is that the variance of any estimator θ̂ must satisfy the Cramér-Rao inequality [17], which reads

The magnitude of the off-diagonal elements determines the extent to which the Cramér-Rao lower bound on a given parameter is affected by the estimation of the other parameters. While these off-diagonal elements generally vanish in the context of single-molecule localisation [20,21], the cross-terms of the information matrix must be considered here. Indeed, the precision of lifetime estimations can be strongly influenced by a lack of information about other parameters. This is notably the case for bi-exponential decay histograms characterized by two decay rates of the same order of magnitude (Γ₁ ∼ Γ₂).

Dimensionless quantities for the parameters involved in the calculations can be obtained by performing the change of variable u = Γ₁t, and by normalizing parameters by N₁ and Γ₁ (this choice was notably made in [29]). Hence, we define the normalized repetition period r = Γ₁T, the number of data points per period k = n/r and the normalized expected number of counts due to background noise β = N_b/(rN₁). As an example, let us consider a common organic dye, Alexa Fluor 488 (τ = 4.1 ns) from which 2, 000 emitted photons have been collected by the detection system. Assuming a repetition rate of 80 MHz, a board resolution of 16 ps and 1, 000 detections due to background noise, the value taken by the dimensionless parameters r, k and β are respectively of the order of r = 3, k = 256 and β = 0.16. Bi-exponential decays can be parametrized using the ratio of the decay rates γ = Γ₂/Γ₁ and the ratio of the expected number of detections η = N₂/N₁. Table 1 summarizes the parameters used in the model. With these parameters, Eq. (3) reads

Table 1. Parameters Involved in the TCSPC Data Model.

View Table

Equation (6) explicitly gives the fundamental limit on the precision of decay rate (and lifetime) estimations. In these expressions, F is calculated by inverting the information matrix and describes the influence of the different parameters involved in the model on the value of the Cramér-Rao bound. A low value of F indicates an experimental setup with a high sensitivity: the F-value is always greater than unity and equals unity when the shot noise limit is reached. Optimisation of a TCSPC setup is therefore achieved when the F-value reaches unity, indicating that the precision of lifetime estimation is limited by the number of detected fluorescence photons. For these reasons, the F-value is used as a figure of merit to quantify the performance of a lifetime imaging technique [31, 32]. In the following sections, we will perform a parametric study of the F-value. To do so, we will consider as a reference situation the ideal case for which k = 500, r = 100, β = 0 and q̃_irf is a Dirac delta function. With these parameters, the F-value is approximately of unity. Each parameter will then be individually varied, in order to highlight the influence of each parameter upon the F-value.

2.3. Mono-exponential case

Let us consider a signal following a mono-exponential distribution (η = 0) with a uniform background noise. The set of unknown parameters is θ = (N, Γ). Figure 1 shows the dependence of the F-value upon the number of counts due to background noise (Fig. 1(a)), the number of data points per lifetime (Fig. 1(b)), the number of fluorescence lifetimes per repetition period (Fig. 1(c)) and the standard deviation of the IRF (Fig. 1(d)) which is assumed to follow a inverse Gaussian distribution. All these parameters fully characterize the experimental conditions, and must be optimized in order to perform shot-noise limited estimations.

 

Fig. 1 Precision of lifetime estimations as a function of (a) the number of counts due to background noise, (b) the number of data point, (c) the repetition period and (d) the standard deviation of the IRF. Each parameter is individually varied with respect to the ideal case in which k = 500, r = 100, β = 0 and q̃_irf is a Dirac delta function. Dashed lines represent the asymptotic values.

Download Full Size | PDF

An analytical expression of F was obtained by Köllner and Wolfrum [29], for the special case in which the IRF is modelled by a Dirac delta function and the number of counts due to background noise is β = 0:

The results shown in Fig. 1 can be straightforwardly applied to identify the parameters limiting the precision of an experimental setup. For instance, if N = 400 fluorescence photons are detected from an emitter, it follows from Eq. (6) that a relative error on Γ of 6% can be reached for F ≤ 1.2. Considering each of the four situations depicted in Fig. 1 individually, such precision can be achieved for a number of counts due to the background noise β ≤ 10⁻² (Fig. 1(a)), for a number of data points per lifetime k ≥ 0.5 (Fig. 1(b)), for a number of fluorescence lifetimes per repetition period r ≥ 4 (Fig. 1(c)), and for a standard deviation of the IRF Γσ_irf ≤ 0.7 (Fig. 1(d)). When these four conditions on the parameters are simultaneously verified, one is ensured to obtain a lower bound on σ_Γ1/Γ₁ smaller than 9.7%.

2.4. Bi-exponential case

Let us now consider a signal following a bi-exponential distribution. If N₂ and Γ₂ can be precisely estimated from independent measurements, the set of unknown parameters is θ = (N₁, Γ₁). Such situation can occur for instance if the second component originates from the luminescence of the substrate, because N₂ and Γ₂ can be estimated by performing a reference measurement without the emitter. In contrast, if the second component originates from the emitter itself, N₂ and Γ₂ are generally estimated from the data, and the set of unknown parameters to be considered becomes θ = (N₁, Γ₁, N₂, Γ₂).

In both cases, the precision of lifetime estimations depends on η = N₂/N₁ and γ = Γ₂/Γ₁. Figure 2 shows the dependence of the F-value upon these parameters for both situations. As we aim at estimating Γ₁, the second decay should be interpreted as a noise source which can only degrade the estimation precision. Hence, the estimation precision is best for N₂ < N₁. Interestingly, when the parameters of the second component can be independently estimated (Fig. 2(a)), the estimation precision is best when Γ₂ > Γ₁, i.e. the lifetime of the second component is shorter than the lifetime to be estimated. Indeed, Γ₁ can be correctly estimated from the last points of the decay histogram, for which the contribution of the second component has vanished (see Section 3 of Appendix: Numerical methods). This arguments also hold when the parameters of the second component must be estimated from the data (Fig. 2(b)) but, in that case, Γ₁ and Γ₂ must be significantly different in order to enable a correct estimation.

 

Fig. 2 Precision of lifetime estimations as a function of η = N₂/N₁ and γ = Γ₂/Γ₁ for a bi-exponential distribution when (a) the parameters N₂ and Γ₂ are known and when (b) these parameters must be estimated from the data. (c) Probability density function associated with the first component (solid lines) and with the second component (dashed lines) for Γ₂ = 2 Γ₁ and Γ₂ = 0.5 Γ₁.

Download Full Size | PDF

As an example, we can investigate the number of fluorescence photons required to obtain a relative error of 6% if the number of photons detected from each component is the same (N₂/N₁ = 1) and for the two situations represented in Fig. 2(c). To begin with, we consider that the parameters of the second component are independently estimated. While N₁ = 381 photons are required if Γ₂ = 2 Γ₁, N₁ = 852 photons are needed if Γ₂ = 0.5 Γ₁. This example illustrates that the situation Γ₂ > Γ₁ is favorable in order to precisely estimate the decay rate Γ₁. As expected, the number of photons required to properly estimate Γ₁ drastically increases when the parameters of the second component must be estimated from the data. Indeed, N₁ = 5, 237 photons are required if Γ₂ = 2 Γ₁ and N₁ = 11, 450 photons are required if Γ₂ = 0.5 Γ₁.

2.5. Numerical tests

In this section, we demonstrate the versatility of the approach by analyzing two examples, one for the estimation of lifetimes on the order of several picoseconds, the other for the estimation of lifetimes on the order of the nanosecond. In both cases, we consider an exponential decay with known background noise, so that the parameters to be estimated are θ = (N, Γ). Each example takes into account an IRF that was experimentally measured on a TCSPC setup. In these two examples, we compare the Cramér-Rao bound on lifetime estimators to numerical results obtained from Monte Carlo experiments. To this end, we numerically generate a set of 10,000 decay histograms for each experimental condition that will be investigated. This is performed by using Eq. (1) to calculate the cumulative distribution function, from which decay histograms can be randomly generated based on the inversion principle [33]. The lifetime τ = 1/Γ is then estimated from each histogram by using a maximum-likelihood (ML) method. The estimation bias is not significant for these numerical experiments (see Section 4 of Appendix: Numerical methods). Moreover, the variance of ML estimators asymptotically approaches the Cramér-Rao bound for large sample statistics, which allows the root-mean square (RMS) deviation of the estimated lifetimes to be close to the fundamental limit for the experimental conditions that are investigated. Note that the implementation of the second example is provided in the open source Code 1 [30].

Precision of picosecond lifetime estimations

We consider an experimental setup designed to estimate picosecond lifetimes, with a board resolution of 1 ps, with a repetition rate of 80 MHz (T = 12.5 ns) and with an IRF characterised by a full with at half maximum (FWHM) of 38 ps. The expected number of fluorescence photons is set to N = 1, 000. Moreover, the background signal is assumed to be uniform, and the expected number of detections due to this background is set to N_b = 500. We calculated the lower bound on the relative standard error σ_τ/τ as a function of the lifetime for these experimental conditions (Fig. 3(a), blue curve) as well as in the case the IRF is a Dirac delta function (Fig. 3(a), green curve). For τ ∼ 1 ns, the lower bound on σ_τ/τ is slightly larger than the shot-noise limited bound due to the background noise. For an ideal IRF, the lower bound on σ_τ/τ slightly decreases when τ decreases since the number of counts per lifetime due to background noise (noted β) decreases. In contrast, when the actual IRF is considered, the lower bound on σ_τ/τ strongly increases for lifetime shorter than half the FWHM, which clearly highlights the relevance of taking into account the IRF for picosecond lifetime estimations. Note that, in this regime, the ML estimator becomes less efficient and slightly deviates from the Cramér-Rao bound.

 

Fig. 3 Relative standard error on lifetime estimates for an experimental setup designed for (a) picosecond and (b) nanosecond lifetime estimations. Solid lines represent the lower bounds on σ_τ/τ and data points are the RMS deviations of lifetime estimates obtained from randomly generated data.

Download Full Size | PDF

Optimal contrast in fluorescence lifetime microscopy

We consider a typical TCSPC setup designed for nanosecond lifetime estimation, with a board resolution of 16 ps, a repetition rate of 80 MHz and an IRF characterized by a FWHM of 240 ps. The expected number of fluorescence photons is set to N = 1, 000. Moreover, we assume that a luminescence background signal is also detected, with the same intensity as the fluorescence signal of interest and a lifetime of 1 ns. In order to choose fluorescent emitters that will lead to the most contrasted FLIM images, we calculated the lower bound on σ_τ/τ as a function of the fluorescence lifetime of the emitter under these experimental conditions (Fig. 3(b), blue curve) as well as in the case of uniform background noise (Fig. 3(b), green curve). While the precision of lifetime estimations is optimal for τ ∼ 400 ps for the uniform background, the precision of lifetime estimations is optimal for τ ∼ 2 ns for the exponential background. Indeed, the precision of lifetime estimations depends on γ (see Fig. 2(a)) and the time dependence of the background noise allows estimations that are more precise for γ > 1, that is, τ > 1 ns. In both situations, for lifetimes larger than 5 ns, the lower bound on σ_τ/τ strongly increases as the number of fluorescence lifetimes per repetition period (noted r) decreases. In this regime, using a laser with a lower repetition rate must be considered in order to improve the precision of lifetime estimations. Moreover, for lifetimes smaller than 100 ps, the lower bound on σ_τ/τ also strongly increases as the fluorescence lifetime becomes comparable to the IRF. In this regime, the precision of the estimations could be improved by using a detection system with a smaller jitter and an excitation laser with a shorter pulse width.

3. Conclusion

In summary, we calculated the lower bound on the standard error on lifetime estimates depending on key experimental parameters. These results can be used as a benchmark for the evaluation of the precision of lifetime estimations. Moreover, they reveal the influence of different parameters upon the estimation precision, providing us with a powerful tool for the optimisation of a TCSPC setup as illustrated by two examples. We notably showed that a significant enhancement of the precision of lifetime estimations can be achieved by choosing the proper fluorescent emitter depending on the expected background noise, the IRF of the setup and the repetition rate of the excitation laser. Various dependent parameters such as the integration time, the power of the excitation laser or the measured spectral range can easily be studied by using the proposed formalism, which can also be extended to the study of time-gated photon counting techniques. We expect these results to be of great interest for current experimental challenges such as the reduction of the acquisition time in FLIM-based techniques, the characterization of photonic antennas with high Purcell factors as well as fluorescence lifetime measurements at the single molecule level.

4. Appendix: Numerical methods

In this section, we give explicit expressions of the Fisher information matrix and of the Cramér-Rao bound, and we provide additional numerical results.

4.1. Calculation of the information matrix

4.1.1. General expression

Using the dimensionless parameters defined in the manuscript, the expectation of each data item reads

(8)$$\begin{array}{r}\hfill f_i=N_1(\sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*e^{-u}\right]\text{d}u+\eta \sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*\gamma e^{-\gamma u}\right]\text{d}u\\ \hfill +\beta r\underset{u_i}{\overset{u_{i+1}}{\int }}{\tilde{q}}_b(u)\text{d}u).\end{array}$$

In order to calculate the information matrix, let us define J^I, J^II, K^I, K^II and J^B as follows:

(9)$$\begin{array}{l}{J}_i^{\text{I}}=\sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*e^{-u}\right]\text{d}u,\\ K_i^{\text{I}}=\sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*(1-u)e^{-u}\right]\text{d}u,\\ J_i^{\text{II}}=\sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*\gamma e^{-\gamma u}\right]\text{d}u,\\ K_i^{\text{II}}=\sum _{l=0}^{+\infty }\underset{u_i+lr}{\overset{u_{i+1}+lr}{\int }}\left[{\tilde{q}}_{\mathit{irf}}(u)*\gamma (1-\gamma u)e^{-\gamma u}\right]\text{d}u,\\ J_i^{\text{B}}=r\underset{u_i}{\overset{u_{i+1}}{\int }}{\tilde{q}}_b(u)\text{d}u.\end{array}$$

With these notations, Eq. (8) reads $f_i=N_1\left(J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}\right)$. Differentiating this expression by each parameters yields

(10)$$\begin{array}{l}\frac{\partial f_i}{\partial N_1}=J_i^{\text{I}},\\ \frac{\partial f_i}{\partial {\mathrm{\Gamma }}_1}=\frac{N_1}{{\mathrm{\Gamma }}_1}{K}_i^{\text{I}},\\ \frac{\partial f_i}{\partial N_b}=J_i^{\text{B}},\\ \frac{\partial f_i}{\partial N_2}=J_i^{\text{II}},\\ \frac{\partial f_i}{\partial {\mathrm{\Gamma }}_2}=\frac{\eta N_1}{\gamma {\mathrm{\Gamma }}_1}{K}_i^{\text{II}}.\end{array}$$

Let us recall the expression of the information matrix:

(11)$${\left[\mathit{ℐ}(\mathit{\theta })\right]}_{jk}=\sum _{i=1}^n\frac{1}{f_i}\left(\frac{\partial f_i}{\partial {\theta }_j}\right)\left(\frac{\partial f_i}{\partial {\theta }_k}\right).$$

where θ = (N₁, Γ₁, N_b, N₂, Γ₂) are the parameters to be estimated. The elements of the information matrix are therefore expressed by

(12)$$\begin{array}{l}{ℐ}_{N_1{N}_1}=\frac{1}{N_1}\sum _{i=1}^n\frac{{\left(J_i^{\text{I}}\right)}^2}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_1}=\frac{N_1}{{\mathrm{\Gamma }}_1^2}\sum _{i=1}^n\frac{{\left(K_i^{\text{I}}\right)}^2}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_1{\mathrm{\Gamma }}_1}=\frac{1}{{\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{I}}{K}_i^{\text{I}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_b{N}_b}=\frac{1}{N_1}\sum _{i=1}^n\frac{{\left(J_i^{\text{B}}\right)}^2}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_1{N}_b}=\frac{1}{N_1}\sum _{i=1}^n\frac{J_i^{\text{I}}{J}_i^{\text{B}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{{\mathrm{\Gamma }}_1{N}_b}=\frac{1}{{\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{B}}{K}_i^{\text{I}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_2{N}_2}=\frac{1}{N_1}\sum _{i=1}^n\frac{{\left(J_i^{\text{II}}\right)}^2}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_1{N}_2}=\frac{1}{N_1}\sum _{i=1}^n\frac{J_i^{\text{I}}{J}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{{\mathrm{\Gamma }}_1{N}_2}=\frac{1}{{\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{II}}{K}_i^{\text{I}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_b{N}_2}=\frac{1}{N_1}\sum _{i=1}^n\frac{J_i^{\text{B}}{J}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{{\mathrm{\Gamma }}_2{\mathrm{\Gamma }}_2}=\frac{{\eta }^2{N}_1}{{\gamma }^2{\mathrm{\Gamma }}_1^2}\sum _{i=1}^n\frac{{\left(K_i^{\text{II}}\right)}^2}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_1{\mathrm{\Gamma }}_2}=\frac{\eta }{\gamma {\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{I}}{K}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{{\mathrm{\Gamma }}_1{\mathrm{\Gamma }}_2}=\frac{\eta N_1}{\gamma {\mathrm{\Gamma }}_1^2}\sum _{i=1}^n\frac{K_i^{\text{I}}{K}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_b{\mathrm{\Gamma }}_2}=\frac{\eta }{\gamma {\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{B}}{K}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}},\\ {ℐ}_{N_2{\mathrm{\Gamma }}_2}=\frac{\eta }{\gamma {\mathrm{\Gamma }}_1}\sum _{i=1}^n\frac{J_i^{\text{II}}{K}_i^{\text{II}}}{J_i^{\text{I}}+\eta J_i^{\text{II}}+\beta J_i^{\text{B}}}.\end{array}$$

Since Γ₁ is our reference, the Cramér-Rao inequality will be most conveniently expressed in terms of this parameter. By inverting the information matrix, we indeed obtain the following expression:

(13)$$\frac{{\sigma }_{{\mathrm{\Gamma }}_1}}{{\mathrm{\Gamma }}_1}\ge \frac{1}{\sqrt{N_1}}\times F(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b),$$

where σ_Γ1 is the standard error on the decay rate estimates and F can be calculated by numerical inversion of the information matrix.

4.1.2. Limiting cases

Ideal IRF

Whenever the IRF can be considered as a Dirac delta function, the coefficients defined by Eq. (9) become

(14)$$\begin{array}{l}{J}_i^{\text{I}}=\left[e^{-u_i}-e^{-u_{i+1}}\right]\sum _{l=0}^{+\infty }{e}^{-lr},\\ K_i^{\text{I}}=\left[u_{i+1}{e}^{-u_{i+1}}-u_i{e}^{-u_i}\right]\sum _{l=0}^{+\infty }{e}^{-lr}+\left[e^{-u_{i+1}}-e^{-u_i}\right]\sum _{l=0}^{+\infty }lr{e}^{-lr},\\ J_i^{\text{II}}=\left[e^{-\gamma u_i}-e^{-\gamma u_{i+1}}\right]\sum _{l=0}^{+\infty }{e}^{-\gamma lr},\\ K_i^{\text{II}}=\left[\gamma u_{i+1}{e}^{-\gamma u_{i+1}}-\gamma u_i{e}^{-\gamma u_i}\right]\sum _{l=0}^{+\infty }{e}^{-\gamma lr}+\left[e^{-\gamma u_{i+1}}-e^{-\gamma u_i}\right]\sum _{l=0}^{+\infty }\gamma lr{e}^{-\gamma lr}.\end{array}$$

$J_i^{\text{B}}$ remains unchanged, as it does not depend on the IRF. These expressions can be further simplified, by using the following properties of geometric series:

(15)$$\begin{array}{l}\sum _{l=0}^{+\infty }{e}^{-lr}=\frac{1}{1-e^{-r}},\\ \sum _{l=0}^{+\infty }lr{e}^{-lr}=\frac{r{e}^{-r}}{{\left(1-e^{-r}\right)}^2}.\end{array}$$

We obtain

(16)$$\begin{array}{l}{J}_i^{\text{I}}=\frac{e^{-u_i}-e^{-u_{i+1}}}{1-e^{-r}},\\ K_i^{\text{I}}=\frac{u_{i+1}{e}^{-u_{i+1}}-u_i{e}^{-u_i}}{1-e^{-r}}+\frac{r{e}^{-r}\left(e^{-u_{i+1}}-e^{-u_i}\right)}{{\left(1-e^{-r}\right)}^2},\\ J_i^{\text{II}}=\frac{e^{-\gamma u_i}-e^{-\gamma u_{i+1}}}{1-e^{-\gamma r}},\\ K_i^{\text{II}}=\frac{\gamma u_{i+1}{e}^{-\gamma u_{i+1}}-\gamma u_i{e}^{-\gamma u_i}}{1-e^{-\gamma r}}+\frac{\gamma r{e}^{-\gamma r}\left(e^{-\gamma u_{i+1}}-e^{-\gamma u_i}\right)}{{\left(1-e^{-\gamma r}\right)}^2}.\end{array}$$

Uniform background noise

Whenever the background noise is uniform over the repetition period, the coefficient $J_i^{\text{B}}$ simply becomes

(17)$$J^{\text{B}}=r/n.$$

4.2. Cramér-Rao bound on the lifetime estimator

From Eq. (13), it is straightforward to calculate the Cramér-Rao bound on the lifetime estimator. To do so, we define the excited-state lifetime τ₁ = 1/Γ₁ and we perform a transformation of parameter, as detailed in Appendix 3B of [17]. This reads

(18)$$\text{Var}({\widehat{\tau }}_1)\ge {\left[\frac{\partial (1/{\mathrm{\Gamma }}_1)}{\partial {\mathrm{\Gamma }}_1}\right]}^2{\left[\frac{{\mathrm{\Gamma }}_1}{\sqrt{N_1}}\times F(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b)\right]}^2,$$

where τ̂₁ is the lifetime estimator. This expression simplifies to

(19)$$\frac{{\sigma }_{{\tau }_1}}{{\tau }_1}\ge \frac{1}{\sqrt{N_1}}\times F(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b),$$

where σ_τ1 is the standard error on the lifetime estimates. This demonstrates that the Cramér-Rao lower bounds are the same for the relative standard error on the decay rate and lifetime estimators.

4.3. Partial and cumulated information

In general, elements of the information matrix are calculated as a sum over the n points of the decay histogram. From Eq. (11), it can therefore be expressed in a formal way as follows:

(20)$$\mathit{ℐ}=\sum _{i=1}^n{\mathit{ℛ}}_i(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b).$$

We define the partial information matrix ℐ_p as the information matrix calculated from the sum over p points of the decay histogram, where p is the number of parameters to be estimated. This reads

(21)$${\mathit{ℐ}}_p=\sum _{i=j}^{j+(p-1)}{\mathit{ℛ}}_i(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b).$$

We also define the cumulated information matrix ℐ_c as the information matrix calculated from the sum over the m first points of the decay histogram, with m ≤ n. This reads

(22)$${\mathit{ℐ}}_c=\sum _{i=j}^m{\mathit{ℛ}}_i(\eta ,\gamma ,r,k,{\tilde{q}}_{\mathit{irf}},\beta ,{\tilde{q}}_b).$$

From these definitions, the partial information and the cumulated information can be used to compare the information carried by different sets of points of the data histogram about the decay rate Γ₁ to be estimated. Figure 4(a) shows the F-value calculated from the partial information matrix as a function of j/k, considering a bi-exponential decay histogram (η = 1) with a large number of data point per lifetime (k = 500), a large number of lifetime per repetition period (r = 100), no background noise (β = 0) and an ideal IRF. The first points of the decay histogram carry more information for γ = 0.2 than for γ = 5. Indeed, in this latter case, the signal due to the second decay is concentrated on the first points of the decay histogram, making the estimation of Γ₁ more difficult. However, this signal vanishes faster and, indeed, we can see that the information contained in the last points of the histogram is larger for γ = 5 than for γ = 0.2. Figure 4(b) shows the F-value calculated from the cumulated information matrix as a function of m/k. For m/k ∼ 2.3, the F-value becomes smaller for γ = 5 than for γ = 0.2, resulting in a better precision for γ = 5 than for γ = 0.2 whenever the whole histogram is considered.

 

Fig. 4 Precision of lifetime estimations calculated for a bi-exponential distribution assuming that the parameters of the second decay component are known. (a) F-value calculated from the partial information matrix, as a function of j/k. (b) F-value calculated from the cumulated information matrix, as a function of m/k.

Download Full Size | PDF

4.4. Bias of maximum-likelihood estimations

The bias B_τ of an estimator τ̂ describes whether the estimator can, on average, recover the true value of the parameter τ. It is defined by

(23)$$B_{\tau }=\text{E}(\widehat{\tau })-\tau .$$

Figure 5 shows the relative bias B_τ/τ of the maximum-likelihood estimator used in the two examples analyzed in the manuscript. In both cases, we can see that the relative bias increases when τ decreases, which can be attributed to the influence of the IRF on the estimation process. However, the relative bias is always one order of magnitude smaller than the relative standard error (Fig. 3 of the manuscript). This justifies the relevance of using the Cramér-Rao bound as a benchmark for the estimation precision (the Cramér-Rao bound applies only to unbiased estimators).

 

Fig. 5 Relative bias on lifetime estimates for an experimental setup designed for (a) picosecond and (b) nanosecond lifetime estimations.

Download Full Size | PDF

Funding

Ville de Paris (Emergences 2015); Agence Nationale de la Recherche (ANR) (ANR-10_LABX-24, ANR-10-IDEX-0001-02 PSL*, ANR-17-CE09-0006).

Acknowledgments

The authors thank Hans Gerritsen for helpful comments about the manuscript and Sebastien Bidault for providing the IRF used in the first example.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. L. Novotny and N. v. Hulst, “Antennas for light,” Nat. Photonics 5, 83 (2011). [CrossRef]  

2. A. F. Koenderink, “Single-photon nanoantennas,” ACS Photonics 4, 710–722 (2017). [CrossRef]  

3. A. Kinkhabwala, Z. Yu, S. Fan, Y. Avlasevich, K. Müllen, and W. E. Moerner, “Large single-molecule fluorescence enhancements produced by a bowtie nanoantenna,” Nat. Photonics 3, 654–657 (2009). [CrossRef]  

4. C. Belacel, B. Habert, F. Bigourdan, F. Marquier, J.-P. Hugonin, S. Michaelis de Vasconcellos, X. Lafosse, L. Coolen, C. Schwob, C. Javaux, B. Dubertret, J.-J. Greffet, P. Senellart, and A. Maitre, “Controlling spontaneous emission with plasmonic optical patch antennas,” Nano Lett. 13, 1516–1521 (2013). [CrossRef]   [PubMed]  

5. G. M. Akselrod, C. Argyropoulos, T. B. Hoang, C. Ciracì, C. Fang, J. Huang, D. R. Smith, and M. H. Mikkelsen, “Probing the mechanisms of large Purcell enhancement in plasmonic nanoantennas,” Nat. Photonics 8, 835–840 (2014). [CrossRef]  

6. T. B. Hoang, G. M. Akselrod, and M. H. Mikkelsen, “Ultrafast room-temperature single photon emission from quantum dots coupled to plasmonic nanocavities,” Nano Lett. 16, 270–275 (2016). [CrossRef]  

7. S. Bidault, A. Devilez, V. Maillard, L. Lermusiaux, J.-M. Guigner, N. Bonod, and J. Wenger, “Picosecond lifetimes with high quantum yields from single-photon-emitting colloidal nanostructures at room temperature,” ACS Nano 10, 4806–4815 (2016). [CrossRef]   [PubMed]  

8. J. Cambiasso, G. Grinblat, Y. Li, A. Rakovich, E. Cortés, and S. A. Maier, “Bridging the gap between dielectric nanophotonics and the visible regime with effectively lossless gallium phosphide antennas,” Nano Lett. 17, 1219–1225 (2017). [CrossRef]   [PubMed]  

9. T. W. J. Gadella, T. M. Jovin, and R. M. Clegg, “Fluorescence lifetime imaging microscopy (FLIM): Spatial resolution of microstructures on the nanosecond time scale,” Biophys. Chem. 48, 221–239 (1993). [CrossRef]  

10. M. Y. Berezin and S. Achilefu, “Fluorescence lifetime measurements and biological imaging,” Chem. Rev. 110, 2641–2684 (2010). [CrossRef]   [PubMed]  

11. W. Becker, “Fluorescence lifetime imaging–techniques and applications,” J. Microsc. 247, 119–136 (2012). [CrossRef]   [PubMed]  

12. Ž. Bajzer, T. M. Therneau, J. C. Sharp, and F. G. Prendergast, “Maximum likelihood method for the analysis of time-resolved fluorescence decay curves,” Eur. Biophys. J. 20, 247–262 (1991). [CrossRef]  

13. J. Tellinghuisen and C. W. Wilkerson, “Bias and precision in the estimation of exponential decay parameters from sparse data,” Anal. Chem. 65, 1240–1246 (1993). [CrossRef]  

14. K. K. Sharman, A. Periasamy, H. Ashworth, and J. N. Demas, “Error analysis of the rapid lifetime determination method for double-exponential decays and new windowing schemes,” Anal. Chem. 71, 947–952 (1999). [CrossRef]   [PubMed]  

15. M. I. Rowley, P. R. Barber, A. C. C. Coolen, and B. Vojnovic, “Bayesian analysis of fluorescence lifetime imaging data,” Proc. SPIE 7903, 790325 (2011). [CrossRef]  

16. B. Kaye, P. J. Foster, T. Y. Yoo, and D. J. Needleman, “Developing and testing a Bayesian analysis of fluorescence lifetime Measurements,” PLOS ONE 12, e0169337 (2017). [CrossRef]   [PubMed]  

17. S. Kay, Cramer-Rao Lower Bound, in Fundamentals of Statistical Processing, Volume I: Estimation Ttheory, (Prentice Hall, 1993), pp. 27–82.

18. R. J. Ober, S. Ram, and E. S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. 86, 1185–1200 (2004). [CrossRef]   [PubMed]  

19. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7, 377–381 (2010). [CrossRef]   [PubMed]  

20. H. Deschout, F. C. Zanacchi, M. Mlodzianoski, A. Diaspro, J. Bewersdorf, S. T. Hess, and K. Braeckmans, “Precisely and accurately localizing single emitters in fluorescence microscopy,” Nat. Methods 11, 253–266 (2014). [CrossRef]   [PubMed]  

21. J. Chao, E. S. Ward, and R. J. Ober, “Fisher information theory for parameter estimation in single molecule microscopy: tutorial,” JOSA A 33, B36–B57 (2016). [CrossRef]   [PubMed]  

22. D. Bouchet, J. Scholler, G. Blanquer, Y. D. Wilde, I. Izeddin, and V. Krachmalnicoff, “Probing near-field light-matter interactions with single-molecule lifetime imaging,” Optica 6, 135–138 (2019). [CrossRef]  

23. L. P. Watkins and H. Yang, “Information bounds and optimal analysis of dynamic single molecule measurements,” Biophys. J. 86, 4015–4029 (2004). [CrossRef]   [PubMed]  

24. M. Boffety, M. Allain, A. Sentenac, M. Massonneau, and R. Carminati, “Analysis of the depth resolution limit of luminescence diffuse optical imaging,” Opt. Lett. 33, 2290–2292 (2008). [CrossRef]   [PubMed]  

25. M. Boffety, M. Allain, A. Sentenac, M. Massonneau, and R. Carminati, “Cramer-Rao analysis of steady-state and time-domain fluorescence diffuse optical imaging,” Biomed. Opt. Express 2, 1626–1636 (2011). [CrossRef]   [PubMed]  

26. E. Gratton, S. Breusegem, J. D. B. Sutin, Q. Ruan, and N. P. Barry, “Fluorescence lifetime imaging for the two-photon microscope: time-domain and frequency-domain methods,” J. Biomed. Opt. 8, 381–391 (2003). [CrossRef]   [PubMed]  

27. K. Suhling, P. M. W. French, and D. Phillips, “Time-resolved fluorescence microscopy,” Photochem. & Photobiol. Sci. 4, 13–22 (2005). [CrossRef]  

28. W. Becker, Advanced Time-Correlated Single Photon Counting Applications (Springer, 2015). [CrossRef]  

29. M. Köllner and J. Wolfrum, “How many photons are necessary for fluorescence-lifetime measurements?” Chem. Phys. Lett. 200, 199–204 (1992). [CrossRef]  

30. D. Bouchet, “Computation of Cramér-Rao lower bound (CRLB) on decay rate estimates,” figshare 2019, https://doi.org/10.6084/m9.figshare.8332733

31. H. Gerritsen, A. Draaijer, D. van den Heuvel, and A. Agronskaia, Fluorescence Llifetime Iimaging in Scanning Microscopy, in Handbook of Biological Confocal Microscopy, J. Pawley, ed. (Springer, 2006), pp. 516–534.

32. J. Philip and K. Carlsson, “Theoretical investigation of the signal-to-noise ratio in fluorescence lifetime imaging,” JOSA A 20, 368–379 (2003). [CrossRef]   [PubMed]  

33. L. Devroye, The Inversion Method, in Non-Uniform Random Variate Generation, (Springer-Verlag, 1986), pp. 27–39.