1.

Introduction

Fluorescence spectroscopy has been extensively explored as a technique for detecting biochemical changes in tissue resulting from pathological transformations.^{1, 2} Tissue autofluorescence is produced by the relative concentration and distribution of endogenous fluorophores. Each fluorophore possesses specific spectral and lifetime (time-resolved) characteristics. Thus, changes in tissue composition will be reflected in the tissue autofluorescence pattern. Time-resolved fluorescence spectroscopy (TRFS) measures fluorescence lifetime as a function of emission wavelength. Lifetime is sensitive to the biological microenvironment and allows temporal discrimination of fluorophores with overlapping emission spectra.³ Since lifetime measurements are independent of the absolute emission intensity, TRFS is less sensitive to intensity artifacts than steady-state measurements; thus, it is more robust for clinical applications. Over the past years, TRFS has been evaluated as a nondestructively and minimally invasive optical technology for clinical diagnosis.^{4, 5, 6, 7, 8, 9, 10, 11} A key element for the translation of TRFS into the clinical arena is the development of robust and automated computational methods for the real-time analysis and interpretation of TRFS data.

In the context of time-resolved time-domain fluorescence measurements, the fluorescence impulse response function (IRF) contains all the temporal information of a single fluorescence decay measurement.^{12, 13} Mathematically, the measured fluorescence intensity decay is the convolution of the fluorescence IRF with the instrument response. Thus, to estimate the fluorescence IRF of a sample, the instrument response must be deconvolved from the measured fluorescence intensity pulse.^{12, 13} The most commonly applied deconvolution method is the so-called least-square iterative reconvolution (LSIR).¹³ In LSIR, the fluorescence IRF is modeled as a sum of exponential terms, and a nonlinear least-square optimization method (i.e., Gauss-Newton or Levenberg-Marquardt) is applied to estimate the parameters of the multiexponential function that would fit best its convolution with the instrument response to the fluorescence decay data. Since the optimization process involves iterative convolutions, LSIR is computationally expensive. Moreover, due to the well known correlation of the fitting parameters in the multiexponential model, different multiexponential models can be fitted equally well to the same fluorescence decay data.³ In the particular case of tissue TRFS data, fluorescence emission typically originates from several endogenous fluorophores, each one potentially presenting complex decay dynamics.³ For such a complex medium, it is not entirely adequate to analyze time-resolved fluorescence decays in terms of multiexponential components, since the resulting exponential terms cannot be directly associated with a particular tissue fluorophore.³ Thus, for tissue TRFS data analysis, there is an advantage in avoiding to predefine any a priori physical model of the fluorescence IRF decay.

An alternative model-free deconvolution method, based on the Laguerre expansion of the kernel (LET) technique, has been recently introduced.^{14, 15} LET was adapted and popularized in the early 1990s by Marmarelis for the linear and nonlinear modeling of physiological systems.¹⁶ LET is based on the expansion of kernels (IRF, for linear systems) on orthonormal sets of discrete Laguerre functions (DLFs), and allows fast converging kernel estimation from short input-output data time series. This method has been extensively applied since then to the modeling of different physiological systems, including cerebral autoregulation and cardiac autonomic control.^{17, 18} Taking advantage of the asymptotic exponential decline characteristics of the DLFs, the LET was most recently adapted for the deconvolution of TRFS decay data and the estimation of fluorescence IRFs.^{14, 15} The resulting Laguerre deconvolution method has been proven to be a fast, robust, and model-free alternative for the analysis of TRFS, properties highly desirable for tissue diagnostic applications. In addition, the Laguerre expansion coefficients derived from this method can be further analyzed to estimate the relative concentration of fluorophores in complex fluorescence systems,¹⁵ and they can be also directly used as features for designing classification algorithms aimed to perform TRFS-based tissue diagnosis.^{8, 9, 10, 11}

One limitation of the original implementation of the Laguerre deconvolution method, however, is the need to choose a priori optimal values of the Laguerre parameter $\alpha$ to guaranty accurate estimations of the fluorescence IRF.^{14, 15} The Laguerre parameter $(0<\alpha <1)$ determines the rate of exponential decline of the DLFs.¹⁶ Thus, small $\alpha$ values are usually needed for expanding short-lived fluorescence decays, while large $\alpha$ values are required for estimating long-lived decays. In most of the previous studies using the Laguerre deconvolution method, the selection of both the Laguerre parameter $\alpha$ and the expansion order was performed heuristically by trial and error, limiting the applicability of this method to offline analysis.^{8, 9, 10, 11} Maarek proposed a method for estimating the Laguerre parameter $\alpha$ , in which simultaneous nonlinear optimization of both $\alpha$ and the Laguerre expansion coefficients is performed.¹⁴ One disadvantage of this method, however, is the unnecessary application of nonlinear optimization for estimating the Laguerre expansion coefficients. Since the expansion coefficients are linear terms in the fluorescence IRF expansion, they can be estimated extremely fast using linear least-square approaches.¹⁵ Thus, we propose a more computationally efficient approach for optimizing the Laguerre parameter, in which nonlinear optimization is only applied to estimate $\alpha$ , while the expansion coefficients are computed using least-square estimation.

The expansion order indicates the number of Laguerre functions needed for accurate modeling of the fluorescence IRF. In general, more DLFs (i.e., higher order) are needed to expand complex decay dynamics; on the other hand, too many DLFs in the expansion can incur model overfitting. Thus, optimal selection the expansion order would also affect the accuracy of the fluorescence IRF estimation. In addition, accurate deconvolution and fluorescence IRF estimation strongly depends on the ability to properly record the instrument response. Ideally, the instrument response should be measured at the same wavelength as the emission fluorescence to minimize any wavelength dependence distortion introduced at the TRFS instrument transmission medium (e.g., optical fiber). In practice, however, the instrument response is usually measured from the scattering of the excitation light pulse traveling the collection optical pathway and reaching the detector, resulting on a zero-time shift or delay of the instrument response with respect to the fluorescence decay.¹⁹ Thus, for exact estimation of the fluorescence IRF and lifetimes, the delay between the scattering-derived instrument response and the fluorescence emission decay has to be properly quantified.

In this work, a new implementation of the Laguerre deconvolution method for fluorescence IRF estimation is introduced. This new automated Laguerre deconvolution approach performs a nonlinear least-square optimization of the Laguerre parameter $\alpha$ , finds the optimal expansion order based on a minimum description length (MDL) criterion, and corrects for the zero-time delay between the recorded instrument response and fluorescence decay based on a normalized mean square error (NMSE) criterion. Results of the method validation on fluorescence lifetimes standards, tissue endogenous fluorphores, and biological tissues are presented. Unlike its previous implementations, the proposed automated Laguerre deconvolution method converges extremely fast to an accurate fluorescence IRF, without the need of selecting a-priori the expansion parameters. Online analysis is thus possible to perform, making our proposed method more suitable for real-time applications, such as tissue diagnosis based on TRFS and/or fluorescence lifetime imaging.

2.

Methods

In this section, the automated Laguerre deconvolution method is first presented. Then, a description of the experiments performed for validating our proposed deconvolution method and comparing it with the standard least-square iterative reconvolution (LSIR) is provided.

2.1.

Automated Laguerre Deconvolution Technique

2.1.1.

Laguerre expansion technique for deconvolution of time-resolved fluorescence spectroscopy data

In the context of TRFS, the measured fluorescence intensity decay data $y\left(n\right)$ is given by the convolution of the fluorescence IRF $h\left(n\right)$ with the instrument response $x\left(n\right)$ :

Eq. 1

$y\left(n\right)=T\cdot \sum _{m=0}^{K-1}h\left(m\right)x(n-m),n=0,\dots ,N-1.$

The parameter $K$ in Eq. 1 determines the time length of the fluorescence IRF, $T$ is the sampling interval, and $N$ is the number of samples measured in $y\left(n\right)$ and $x\left(n\right)$ . The Laguerre deconvolution technique expands the fluorescence IRF on an orthonormal set of discrete time Laguerre functions (DLF) $b_j^{\alpha }\left(n\right)$ :

Eq. 2

$h\left(n\right)=T.\sum _{j=0}^{L-1}{c}_j^{\alpha }{b}_j^{\alpha }\left(n\right).$

In Eq. 3, $c_j^{\alpha }$ are the unknown Laguerre expansion coefficients (LEC), which are to be estimated from the input-output data; $b_j^{\alpha }\left(n\right)$ denotes the $j$ ’th order orthonormal DLF;¹⁶ $L$ is the number of DLFs used to model the IRF, thus defining the order of the expansion. The Laguerre parameter $\alpha$ $(0<\alpha <1)$ determines the rate of asymptotic decline of the DLFs. Thus, larger $\alpha$ values imply longer convergence time to zero. By inserting Eq. 2 into Eq. 1, the convolution Eq. 1 becomes:

3.

$y\left(n\right)=\sum _{j=0}^{L-1}{c}_j^{\alpha }{v}_j^{\alpha }\left(n\right)$

$v_j^{\alpha }\left(n\right)=\sum _{i=0}^{K-1}{b}_j^{\alpha }\left(m\right)x(n-m).$

The functions $v_j^{\alpha }\left(n\right)$ , representing the digital convolution of the input $x\left(n\right)$ with each of the Laguerre functions, are denoted as the “key variables.” The system of linear Eq. 3 can be expressed in a matrix notation as follows:

Eq. 4

Finally, the expansion coefficients ${\widehat{c}}_{\alpha }$ can be estimated by least-square estimation as:

Eq. 5

${\widehat{c}}_{\alpha }={\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}.$

Once the expansion coefficient ${\widehat{c}}_{\alpha }$ has been estimated, the fluorescence IRF $h\left(n\right)$ can be computed from Eq. 2 and the average lifetimes can be calculated as:³

Eq. 6

$\tau =\frac{T.{\sum }_{n=0}^{N}n.h\left(n\right)}{{\sum }_{n=0}^{N}h\left(n\right)}.$

2.1.2.

Nonlinear least square optimization of the Laguerre parameter

The choice of the Laguerre parameter $\alpha$ is critical in achieving accurate fluorescence IRF expansions. As a first approach, the parameter $\alpha$ can be selected based on the kernel memory length $K$ and the number of DLFs used for the expansion, so that all the functions converge to zero by the end of the impulse response.^{15, 16} However, this empirical approach does not warranty an optimal expansion, and in practice a good choice for the value of $\alpha$ is usually found by trial-and-error procedures. Here, we propose a computationally efficient method, where the Laguerre parameter is treated as a free parameter within a nonlinear least-square optimization scheme. This automates the procedure for the determination of suitable Laguerre parameters, guided by the actual TRFS experimental data.

In the context of nonlinear least-square optimization, the model free parameter ( $\alpha$ in this case) is chosen so that an objective function is minimized. If we write Eq. 4 as:

Eq. 7

$\overline{y}=V_{\alpha }{\overline{c}}_{\alpha },$

then we know that the least-Squares solution of this equation is given by:

Eq. 8

$\widehat{y}=V_{\alpha }{\widehat{c}}_{\alpha }=V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}.$

We now define the error in our estimation as:

Eq. 9

${\epsilon }_N=\overline{y}-\widehat{y}=\overline{y}-V_{\alpha }{\widehat{c}}_{\alpha }=\overline{y}-V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y},$

hence our cost function can be defined as:

Eq. 10

$F\left(\alpha \right)={\epsilon }_N^T{\epsilon }_N.$

Our task is to find a value of $\alpha$ that minimizes the prior cost function:

Eq. 11

$\widehat{\alpha }=\arg {\min }_{\alpha }F\left(\alpha \right).$

Since the cost function is a nonlinear function of $\alpha$ , minimization requires an iterative scheme for determining $\widehat{\alpha }$ . Thus, an iterative scheme for $\alpha$ is defined as:

Eq. 12

${\widehat{\alpha }}^{(k+1)}={\widehat{\alpha }}^{\left(k\right)}-{\gamma }^{\left(k\right)}{\left[\frac{dF\left(\alpha \right)}{d\alpha }\right]}^{\left(k\right)},$

where ${\gamma }^{\left(k\right)}$ is a scaling factor chosen in a way such that the value of $\widehat{\alpha }$ is incremented by a fixed amount after each iteration. Note that the superscript $k$ denotes the iteration number. Estimation of the Jacobian $dF\left(\alpha \right)∕d\alpha$ can be achieved by numerical methods. Alternatively, closed form recursive formulation of $dF\left(\alpha \right)∕d\alpha$ can also be derived as follows.

Substituting ${\epsilon }_N$ from Eq. 9 into Eq. 10, we get:

Eq. 13

$F\left(\alpha \right)={[\overline{y}-V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}]}^T[\overline{y}-V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}].$

This can be further simplified to:

Eq. 14

$F\left(\alpha \right)={\overline{y}}^T[I_N-V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T]\overline{y}={\overline{y}}^T[\overline{y}-V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}],$

where $I_N$ is the identity matrix of order $N$ . By identifying the second term on the RHS as ${\epsilon }_N$ , it can be further simplified to:

Eq. 15

$F\left(\alpha \right)={\overline{y}}^T{\epsilon }_N.$

Thus, the Jacobian can be simply expressed as:

Eq. 16

$\frac{dF\left(\alpha \right)}{d\alpha }={\overline{y}}^T\frac{d\left({\epsilon }_N\right)}{d\alpha }.$

Using the definition of ${\epsilon }_N$ as in Eq. 9, its derivative with respect to $\alpha$ becomes:

Eq. 17

$\frac{d\left({\epsilon }_N\right)}{d\alpha }=\frac{d(\overline{y}-V_{\alpha }{\widehat{c}}_{\alpha })}{d\alpha }=-\frac{d\left(V_{\alpha }{\widehat{c}}_{\alpha }\right)}{d\alpha }=-[\frac{d\left(V_{\alpha }\right)}{d\alpha }{\widehat{c}}_{\alpha }+V_{\alpha }\frac{d\left({\widehat{c}}_{\alpha }\right)}{d\alpha }].$

The key variables $v_j^{\alpha }\left(n\right)$ forming the matrix $V_{\alpha }$ can be computed recursively using the following formulation:¹⁶

18.

$v_j^{\alpha }\left(n\right)=\sqrt{\alpha }{v}_j^{\alpha }(n-1)+\sqrt{\alpha }{v}_{j-1}^{\alpha }\left(n\right)-v_{j-1}^{\alpha }(n-1),$

$v_0^{\alpha }\left(n\right)=\sqrt{\alpha }{v}_0^{\alpha }(n-1)+\sqrt{1-\alpha }x\left(n\right).$

From here, we can derive a recursive formulation for the terms in $d\left(V_{\alpha }\right)∕d\alpha$ :

19.

$\frac{d{v}_j^{\alpha }\left(n\right)}{d\alpha }=\frac{1}{2\sqrt{\alpha }}[v_j^{\alpha }(n-1)+v_{j-1}^{\alpha }\left(n\right)]+\sqrt{\alpha }[\frac{d{v}_j^{\alpha }(n-1)}{d\alpha }+\frac{d{v}_{j-1}^{\alpha }\left(n\right)}{d\alpha }]-\frac{d{v}_{j-1}^{\alpha }(n-1)}{d\alpha },$

$\frac{d{v}_0^{\alpha }\left(n\right)}{d\alpha }=\frac{1}{2\sqrt{\alpha }}\left[v_0^{\alpha }(n-1)\right]+\sqrt{\alpha }\left[\frac{d{v}_0^{\alpha }(n-1)}{d\alpha }\right]+\frac{1}{2\sqrt{1-\alpha }}x\left(n\right).$

Similarly, by taking the derivative with respect to $\alpha$ to the least-square solution of the expansion coefficients as defined in Eq. 5, we obtain:

Eq. 20

$\frac{d{\widehat{c}}_{\alpha }}{d\alpha }=\frac{d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}{V}_{\alpha }^T\overline{y}\right]}{d\alpha }=\{\frac{d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\right]}{d\alpha }{V}_{\alpha }^T\overline{y}+{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\frac{d\left[V_{\alpha }^T\overline{y}\right]}{d\alpha }\}.$

To calculate $d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\right]∕d\alpha$ , let us first consider that:

$\frac{d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\left(V_{\alpha }^T{V}_{\alpha }\right)\right]}{d\alpha }=\frac{d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\right]}{d\alpha }\left(V_{\alpha }^T{V}_{\alpha }\right)+{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\frac{d\left[\left(V_{\alpha }^T{V}_{\alpha }\right)\right]}{d\alpha }.$

Since ${\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\left(V_{\alpha }^T{V}_{\alpha }\right)=I_N$ , the LHS becomes zero, thus:

Eq. 21

$\frac{d\left[{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\right]}{d\alpha }=-{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\frac{d\left[\left(V_{\alpha }^T{V}_{\alpha }\right)\right]}{d\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}.$

Substituting Eqs. 5, 21 into Eq. 20, we get:

Eq. 22

$\frac{d{\widehat{c}}_{\alpha }}{d\alpha }={\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\{\frac{d\left(V_{\alpha }^T\right)}{d\alpha }\overline{y}-[\frac{d\left(V_{\alpha }^T\right)}{d\alpha }{V}_{\alpha }+V_{\alpha }^T\frac{d\left(V_{\alpha }\right)}{d\alpha }]{\widehat{c}}_{\alpha }\}.$

By inserting Eq. 22 into Eq. 17, the Jacobian can be expressed as:

Eq. 23

$\frac{dF\left(\alpha \right)}{d\alpha }=-{\overline{y}}^T(\frac{d\left(V_{\alpha }\right)}{d\alpha }{\widehat{c}}_{\alpha }+V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\{\frac{d\left(V_{\alpha }^T\right)}{d\alpha }\overline{y}-[\frac{d\left(V_{\alpha }^T\right)}{d\alpha }{V}_{\alpha }+V_{\alpha }^T\frac{d\left(V_{\alpha }\right)}{d\alpha }]{\widehat{c}}_{\alpha }\}).$

Finally, our last iterative relation for calculating $\widehat{\alpha }$ becomes:

Eq. 24

${\widehat{\alpha }}^{(k+1)}={\widehat{\alpha }}^{\left(k\right)}-{\gamma }^{\left(k\right)}{\overline{y}}^T{(\frac{d\left(V_{\alpha }\right)}{d\alpha }{\widehat{c}}_{\alpha }+V_{\alpha }{\left(V_{\alpha }^T{V}_{\alpha }\right)}^{-1}\{\frac{d\left(V_{\alpha }^T\right)}{d\alpha }\overline{y}-[\frac{d\left(V_{\alpha }^T\right)}{d\alpha }{V}_{\alpha }+V_{\alpha }^T\frac{d\left(V_{\alpha }\right)}{d\alpha }]{\widehat{c}}_{\alpha }\})}^{\left(k\right)},$

where the terms of $d\left(V_{\alpha }\right)∕d\alpha$ and $d\left(V_{\alpha }^T\right)∕d\alpha$ can be calculated using the recursive relations in Eq. 19. Iterative updating of the Laguerre parameter $\widehat{\alpha }$ is performed until the cost function $F\left(\alpha \right)$ converges to a global minimum.

2.1.3.

Expansion order selection

The deconvolution approach based on the Laguerre expansion technique can be taken as a system identification problem, in which the model complexity (i.e., number of Laguerre functions used in the expansion) is to be determined. Increasing the model complexity will decrease the systematic errors. However, at a certain complexity, additional model parameters no longer reduce the systematic errors but start to follow the actual noise realization on the data, a phenomenon known as overfitting. To avoid this unwanted behavior, a model selection criterion usually includes a model complexity term to penalize for overfitting conditions. Among the various model selection criterion used, one of the most robust and popular ones is the minimum description length (MDL),²⁰ which is defined as:

Eq. 25

$\mathrm{MDL}\left(L\right)=\frac{F\left(\alpha \right)}{N}\exp \left[\frac{\ln \left(N\right)L}{N}\right].$

Here, $N$ is the number of measured samples, $F\left(\alpha \right)$ is the sum of squared errors as defined in Eq. 10, and $L$ is the number of model parameters, which in our case corresponds to the expansion order (number of DLFs used in the expansion). To estimate the optimal expansion order, $L$ takes values from 2 to 8, and $\mathit{MDL}\left(L\right)$ is computed. The value of $L$ that yields the minimum $\mathit{MDL}$ value is then chosen as the optimal expansion order.

2.1.4.

Zero-time delay estimation

As discussed before, when the instrument response is recorded experimentally from the scattering of the excitation light pulse, a zero-time shift or delay of the instrument response with respect to the fluorescence decay is often observed.¹⁹ Here, we describe a systematic method to accurately estimate the delay between the scattered instrument response and fluorescence emission decay. Basically, a number of potential delay values is assumed, and for each one the Laguerre deconvolution is applied using a fixed order (e.g., 4), and the normalized mean-square error (NMSE) values is computed. Then, the value yielding the lowest NMSE is taken as the optimal delay between the instrument response and the analyzed fluorescence decay at a given emission wavelength.

3.

Results

3.1.

Validation on Fluorescence Lifetime Standards and Tissue Endogenous Fluorphores

The proposed automated deconvolution method was first validated on the TRFS measurements from fluorescence lifetime standards. Results for this first validation are summarized in Fig. 1 and Table 1 . Figure 1 shows the fluorescence pulses of the standards measured at their peak emission wavelength, together with the instrument response, and the estimated fluorescence pulses by the automated Laguerre and the multiexponential deconvolution methods. The estimated fluorescence IRF, the normalized residuals, and the autocorrelation of the residuals are also shown. The accurate deconvolution of the fluorescence decay of 9-cyanoanthracene (9CA) in ethanol demonstrated the capability of the automated Laguerre technique to estimate long fluorescence lifetimes [Figs. 1a and 1b]. Both the Laguerre and exponential methods yielded similar fluorescence IRFs with lifetime values in agreement with previous reports (Table 1).^{3, 22, 23, 24} The normalized residuals $(<10\%)$ and their autocorrelation functions showed a significant random behavior, indicating excellent performance by the two methods. Optimal values of $\alpha =0.96$ and order $L=4$ were obtained, and a single exponential component was needed for proper multiexponential deconvolution. Similarly, measurements of rhodamin B in ethanol demonstrated the ability of the automated Laguerre method to accurately resolve nanosecond fluorescence lifetimes [Figs. 1c and 1d]. Again, both methods yielded similar fluorescence IRFs with lifetime values also in agreement with previous reports (Table 1), and excellent estimation accuracy was also evident by the low normalized residual $(<5\%)$ and flat autocorrelation function. For this case, optimal values of $\alpha =0.84$ and order $L=4$ were obtained, and again one exponential component was also needed for proper multiexponential deconvolution. Short-lived fluorescence IRF, like that from rose Bengal in methanol, with lifetimes ranging in the few hundreds of picoseconds could also be reliably retrieved by our proposed deconvolution technique [Figs. 1e and 1f]. Optimal values of $\alpha =0.57$ and order $L=3$ were obtained, and one exponential component was also used for this fluorophore. As in the previous cases, both methods yielded similar fluorescence IRFs with lifetime values in agreement with the literature (Table 1).

Fig. 1

Deconvolution results of the TRFS data from the fluorescence lifetime standards (a) and (b) 9CA (at 450 nm), (c) and (d) rhodamin B (at 580 nm), and (e) and (f) rose Bengal (at 570) using the automated Laguerre (left) and multiexponential (right) methods. Main panels show the measured fluorescence pulse (solid black), instrument response (solid gray), and estimated fluorescence pulse (dotted black). Smaller panels show the fluorescence IRF (top), normalized residuals (middle), and residual autocorrelation (bottom). The fluorescence IRFs were accurately estimated by both methods.

Table 1

Validation results from fluorescence standards and tissue endogenous fluorphore TRFS data. Computation time for LSIR is a single exponential decay function for the first five entries. The last three show a double exponential decay function. Values are shown s ( mean±SE ) from a total of five measurements per wavelength.

Sample	Solvent	Wavelength(nm)	Optimalalpha	Optimalorder	Lifetime (ns)			Computation time (ms)
					Opt. LAG	LSIR	Literature	LAG	Opt. LAG	LSIR
9-CA	Ethanol	450	$0.96\pm 0.02$	4	$12.61\pm 0.58$	$12.18\pm 0.73$	11.7 to 12.28	7	266	153
Rhodamin B	Ethanol	580	$0.84\pm 0.01$	4	$2.52\pm 0.10$	$2.44\pm 0.14$	2.60 to 3.01	5	258	67
Rhodamin B	${\mathrm{H}}_2\mathrm{O}$	580	$0.74\pm 0.01$	4	$1.55\pm 0.03$	$1.58\pm 0.05$	1.48 to 1.67	5	193	68
Rose Bengal	Ethanol	570	$0.62\pm 0.01$	3	$0.65\pm 0.02$	$0.61\pm 0.01$	0.84	3	69	45
Rose Bengal	Methanol	570	$0.57\pm 0.01$	3	$0.48\pm 0.01$	$0.48\pm 0.01$	0.54	3	32	73
Collagen	N/A	400	$0.91\pm 0.02$	3	$3.49\pm 0.14$	$3.61\pm 0.15$	2 to 4	3	175	366
NADH	PBS	450	$0.67\pm 0.01$	4	$0.45\pm 0.03$	$0.36\pm 0.01$	0.30 to 0.50	4	54	139
FAD	PBS	530	$0.93\pm 0.02$	3	$3.11\pm 0.13$	$3.17\pm 0.12$	3 to 4	3	157	779

Following the validation with fluorescence lifetime standards, we also measured and analyzed TRFS data from three important tissue constituents relevant for potential fluorescence-based tissue diagnosis: collagen, NADH, and FAD. Results for this second validation are summarized in Fig. 2 and Table 1. Both the automated Laguerre and multiexponential methods accurately estimated the collagen fluorescence IRF and lifetime values at its peak emission wavelength of $400\mathrm{nm}$ , as shown in Figs. 2a and 2b, and Table 1. From the residuals plots, however, it seemed that the Laguerre approach performed better ( $<5\%$ error) than the multiexponential deconvolution ( $<10\%$ error). Optimal values of $\alpha =0.91$ and order $L=3$ were obtained, and two exponential components were needed for proper multiexponential deconvolution ( $\mathrm{tau}1=4.82$ , $\mathrm{tau}2=1.21$ , $\mathrm{a}1∕\mathrm{a}2=0.34$ ). The NADH fluorescence IRF and lifetime value at the peak emission of $450\mathrm{nm}$ were also properly estimated by both methods, as reflected in Figs. 2c and 2d, and Table 1. For this fluorophore, optimal values of $\alpha =0.67$ and order $L=4$ were obtained, and one exponential component was used. Similarly, both methods accurately estimated the fluorescence IRF and lifetime of FAD at its peak emission wavelength of $530\mathrm{nm}$ . Optimal values for $\alpha =0.93$ and order $L=3$ were obtained, and two exponential terms were used ( $\mathrm{tau}1=3.32$ , $\mathrm{tau}2=0.81$ , $\mathrm{a}1∕\mathrm{a}2=3.87$ ).

Fig. 2

Deconvolution results of the TRFS data from the tissue constituents (a) and (b) collagen (at 400 nm), (c) and (d) NADH (at 450 nm), and (e) and (f) FAD (at 530 nm) using the automated Laguerre (left) and multiexponential (right) methods. The fluorescence IRFs were accurately estimated by both methods.

To illustrate the sensitivity of the figure of merits used for optimizing orders and delays, we show in Fig. 3 the MDL value as a function of order, and the NMSE value as a function of delay, corresponding to the FAD fluorescence decay deconvolution. The left panel shows how the MDL value decreases significantly from an order 2 to an order 3 estimation, but does not change significantly for orders higher than 3. Based on this, an optimal order $L=3$ was successfully found. The right panel shows how the NMSE reaches a global minimum at a delay of $0.8\mathrm{ns}$ , clearly showing that this is the zero-time difference between the measured fluorescence pulse and instrument response. These results indicate that the MDL and NMSE values are adequate criteria for searching optimal expansion order and zero-time delay, respectively.

Fig. 3

Values of MDL versus order (a) and NMSE versus zero-time delay (b) corresponding to the FAD fluorescence decay deconvolution. MDL values decreased drastically after order 2 and remained low afterward; an optimal order of 3 was selected. The NMSE reached a global minimum at $0.8\text{-}\mathrm{ns}$ delay.

3.2.

Validation on Time-Resolved Fluorescence Microscopy Ex Vivo Measurement from Human Arteries

Finally, the proposed automated Laguerre deconvolution method was also validated on TRFS data obtained from human arteries. Deconvolution results from a sample artery fluorescence decay are shown in Fig. 4 . For this particular case, optimal values at 400 nm for $\text{alpha}=0.88$ and order $L=4$ were obtained. Two exponential terms were used for the LSIR deconvolution of all artery TRFS data. Based on the low residuals level $(<5\%)$ and flat autocorrelation function, it can be observed that both the Laguerre and multiexponential methods successfully deconvolved the artery TRFS data. The lifetime values obtained from the arterial TRFS data ranged from $\sim 1.5\text{to}3\mathrm{ns}$ , which is consistent with values reported in the literature.^{8, 9, 11}

Fig. 4

Deconvolution results from a sample of human artery TRFS data using the automated Laguerre (a) and multiexponential (b) methods. The fluorescence IRFs were accurately estimated by both methods.

In addition, we wanted to assess ranges of values of the parameter $\alpha$ , the order $L$ , and the delay that were optimal for deconvolving human arterial TRFS data. For this purpose, we estimated the average values of these parameters $(\text{mean}\pm \mathrm{SE})$ as a function of wavelengths from all datasets collected. Results of this analysis are shown in Fig. 5 . It was observed that optimal $\alpha$ values for accurate deconvolution of arterial TRFS data were between 0.90 to 0.94, and there was no obvious variation between $\alpha$ values with emission wavelengths [Fig. 5a]. The optimal orders were mostly between 4 and 5, and there was no apparent correlation between them and the emission wavelengths [Fig. 5b]. The zero-shift delay, on the other hand, was highly correlated with emission wavelength, going from $\sim 0.5\mathrm{ns}$ at $400\mathrm{nm}\text{to}\sim 0.95\mathrm{ns}$ at $550\mathrm{nm}$ , with an almost linear increase of delay with the emission wavelength [Fig. 5c]. We also estimated the average NMSE as a function of emission wavelength, which indicates excellent estimation performance for all wavelengths $(\mathrm{NMSE}<2\%)$ [Fig. 5d].

Fig. 5

Average ( $\text{mean}\pm \text{SE}$ , from 47 datasets) values for (a) the Laguerre parameter $\alpha$ , (b) expansion order, (c) zero-time delay, and (d) NMSE as a function of emission wavelengths from all human artery TRFS datasets. Values of $\alpha$ were between 0.90 and 0.94 and did not change with wavelength. Expansion orders were between 3 and 5 and did not change with wavelength either. Zero-time delay increases with wavelength. NMSE was less than 5% and increased above $500\mathrm{nm}$ due to lower SNR at higher wavelengths.

4.

Discussion

We introduce a fully automated deconvolution method for online TRFS data analysis based on an iterative implementation of the Laguerre deconvolution technique. In this new implementation, the optimization of the expansion parameters (i.e., Laguerre $\alpha$ parameter and expansion order) is performed as part of the deconvolution routine. In addition, this method provides a means for compensating for the well known wavelength-dependence delay between the instrument response and the experimentally measured fluorescence decay. Since optimal estimation of the fluorescence IRF is performed automatically, this method offers the opportunity to realize online analysis of TRFS data and will thus facilitate the use of time-resolved fluorescence spectroscopy in real-time applications. The following is a discussion of the results validating this method, and its advantages over previous implementa-tions of the Laguerre deconvolution technique and the more standard LSIR multiexponential method.

The Laguerre parameter $\alpha$ defines an orthonormal set of DLFs by specifying their rate of exponential decline: fast or slow converging DLFs can be defined by choosing a small or large $\alpha$ value, respectively. Thus, fast or slow decaying fluorescence IRFs can be optimally expanded on a set of DLFs defined by a small or large $\alpha$ parameter value, respectively. The proposed automated Laguerre deconvolution method performs a nonlinear least-square search for the optimal values of $\alpha$ to guarantee that the fluorescence IRF is properly expanded and estimated. Our validation results on TRFS data from fluorescence lifetime standards and tissue endogenous fluorophores showed that our proposed algorithm is capable of converging first to an optimal value of $\alpha$ using nonlinear least-square optimization, and then to an optimal fluorescence IRF expansion using linear least-square estimation of the expansion coefficients. This is a significant improvement over previous implementation of the Laguerre deconvolution method, in which the values of $\alpha$ are determined either by trial and error or together with the expansion coefficients using nonlinear optimization as well. The first approach is computationally efficient; however, since the user is required to evaluate different values of $\alpha$ until finding an optimal one, this approach is unsuitable for online applications.¹⁵ The second implementation, proposed by Maarek,¹⁴ provides a method for estimating both the parameter $\alpha$ and the IRF expansion; however, it is computationally very expensive (nonlinear optimization is applied not only for the $\alpha$ parameter but for the expansion coefficients, too), thus is not suitable for online analysis either. By separating the optimization of $\alpha$ from the estimation of the expansion coefficients, our method provides a faster alternative for performing fully automated deconvolution, which is suitable for online TRFS data analysis.

Model complexity quantification is a very well studied problem in the field of mathematical modeling and system identification.²⁰ A number of model complexity figures have been proposed to determine the level of complexity needed by a given model to properly account for the dynamics being investigated. Among the numerous complexity figures being proposed, MDL is one of the most robust and widely used.²⁰ The Laguerre deconvolution method constitutes a modeling problem, in which the fluorescence IRF is modeled as a linear expansion of Laguerre functions. As such, the model complexity is directly linked to the expansion order $L$ . Our results clearly demonstrated that the MDL can be successfully used for quantifying the expansion complexity, therefore determining the optimal expansion order needed for accurate fluorescence IRF expansion. As shown in Fig. 3a as an example, the MDL value usually decreases as the order expansion increases, until an optimal order is reached, after which the MDL value does not further decrease significantly. Since the MDL criterion sacrifices model estimation for level of complexity, in some instances it is convenient to choose a higher order than the one determined by MDL to warranty good fluorescence IRF estimation. Such heuristic rules can be easily implemented on the proposed algorithm.

During TRFS measurement, the instrument response cannot always be acquired at the same emission wavelength as the recorded fluorescence decay.¹⁹ Since the refractive index of the various optical materials present in a TRFS instrument depend on the photon wavelengths, the instrument response and fluorescence pulses may travel at different speeds throughout the TRFS instrumentation. Because of this speed difference, the instrument response and the measured fluorescence pulse may arrive with a relative time delay between each other. Thus, during the deconvolution process, it may be necessary to account for this zero-time shift or delay. Making the assumption that the optimal fluorescence IRF estimation is achieved when both the instrument response and fluorescence pulse are aligned in time (zero delay), we decided to use a goodness of the fit index, such as the NMSE, to estimate the experimental zero-time delay. As shown in Fig. 3b as an example, the NMSE values usually show a global minimum at the optimal delay, allowing a systematic determination of the zero-time shift.

For the past several years, TRFS has been evaluated as a potential clinical tool for minimally invasive and nondestructive tissue diagnosis. Any successful clinical diagnostic instrument based on TRFS has to be capable of acquiring and processing the fluorescence data from the tissue in quasi-real-time to provide some type of clinical feedback to the clinician. To this regard, an automated method for processing TRFS data online will facilitate the translation of this technology into the clinical arena. Our validation results on human artery TRFS data presented here demonstrate the capability of our method to automatically deconvolve time-resolved fluorescence data from tissue. Our method can thus be eventually implemented on a clinical TRFS instrument to perform the online processing of the fluorescence data needed for diagnosis evaluation.

Our group is particularly interested in using time-resolved fluorescence measurements (spectroscopy and imaging) for quantifying the biochemical composition of atherosclerotic plaques, which is clinically relevant for detecting those patients with high risk of heart attacks and strokes. To this regard, we were also interested in investigating what values of the parameter $\alpha$ were optimal for accurate deconvolution of human artery TRFS data. Our results showed that optimal values for $\alpha$ were concentrated within a range of 0.90 to 0.94 with no apparent correlation between $\alpha$ and the emission wavelength [Fig. 5a]. These findings indicate that it may be possible to speed up the optimization of the parameter $\alpha$ , for instance, by assigning an initial condition within that interval (e.g., 0.92) for the nonlinear least-square algorithm. In an extreme case, it may be even feasible to use a fixed suboptimal valued for $\alpha$ (e.g., $\alpha =0.92$ ) for deconvolving any arterial TRFS. This hypothesis will be tested as part of our future research.

Similarly, we also investigated how many Laguerre functions were necessary for accurate estimation of arterial-fluorescence IRF. Our results showed that optimal expansion orders were mostly between 4 and 5, regardless of the emission wavelength [Fig. 5b]. These findings suggest that we can limit our search for the optimal order in a very reduced interval (e.g. 4 to 5), which would significantly increase the computational speed of our method. As in the case of the parameter $\alpha$ , in the extreme case, it may be even suitable to fix the expansion order to an upper bound of the optimal value (e.g., $L=5$ ) to ensure the goodness of the fit by slightly sacrificing model complexity. Again, this hypothesis will be tested as part of our future research.

As previously discussed, due to the wavelength-dependent refractive index of the optical materials traveled by light pulses in the TRFS instruments, it is common to observe a zero-time delay between the recorded instrument response and measured fluorescence pulses. In our particular TRFS instrument, both the scattering light of the excitation laser (instrument response) and the sample fluorescence emission are collected through a silica/silica step index optical fiber bundle. Since the refractive index of this fiber is higher for shorter than for longer wavelengths, fluorescence pulses at longer wavelengths than the excitation wavelength travel faster than the measured scattered instrument response. Thus, the instrument response is usually delayed with respect to the fluorescence pulse, and the relative delay is expected to increase for longer fluorescence emission wavelengths. Our results from arteries are in agreement with this phenomenon, as it was observed from how the relative delay of the instrument response increases with the emission wavelength [Fig. 5c]. These results underscore the need for taking into consideration the zero-time delay for proper estimation of the fluorescence IRF, especially for TRFS instruments with high temporal resolution.

The validation results on the human artery, taken together, are of special relevance for the application of TRFS in clinical tissue diagnosis. The findings discussed suggest the possibility of customizing our proposed automated Laguerre deconvolution method for online analysis of TRFS data from a particular type of tissue (e.g., coronary arteries or epithelial tissue). For instance, a TRFS instrument could be characterized and calibrated to determine suboptimal values for zero-time delay at specific wavelengths, and would allow the algorithm to search around these for optimal delay values. In addition and as discussed before, a fixed expansion order could be predefined as well. Finally, since optimal $\alpha$ values for deconvolving TRFS data of a particular tissue are expected to be limited to a narrow range, a customized and faster version of our algorithm can be implemented to search only in that relevant interval. Similar approaches customizing the proposed automated Laguerre deconvolution method for a particular application will also be useful for facilitating online applications of TRFS in numerous other areas (e.g., chemistry, biochemistry, and drug development).

Perhaps the only disadvantage of the proposed automated Laguerre deconvolution over the previous heuristic implementation extensively discussed in a previous publication¹⁵ is its slower performance, as can be observed in Table 1. This is understandable, since a full model parameter optimization is performed in the automated version: the parameter $\alpha$ is optimized by nonlinear least-square minimization, and both the order and delay values are optimized by exhaustive searching of values minimizing complexity (MDL) and performance (NMSE) figures of merit, respectively. On the other hand, no user intervention is required for tuning up the model parameters, allowing online analysis of TRFS data. Furthermore, and as discussed before, the method can be customized for a given sample type and TRFS instrumentation to speed up its performance. Also interesting was the fact that the proposed automated Laguerre deconvolution performed at least as fast as the standard LSIR method. This result indicates that all the intrinsic advantages of the Laguerre deconvolution with respect to the LSIR method can still be delivered without significant increase in computational time.

Our future efforts include extending the automated Laguerre deconvolution method for fluorescence lifetime imaging microscopy (FLIM) data analysis. In the context of FLIM, the analysis of a single image requires deconvolution of the instrument response from the fluorescence decay of each pixel within the field of view. We have recently introduced and validated a Laguerre expansion-based approach for FLIM data deconvolution that performs at least 2 orders of magnitude faster than standard FLIM deconvolution algorithms.²⁵ As for the case of the original Laguerre TRFS deconvolution, the current Laguerre FLIM deconvolution also requires to select the expansion parameters (i.e., $\alpha$ and order) by trial and error. Thus, an automated version of the Laguerre FLIM deconvolution method would facilitate the use of this novel technique in various applications, including our current research on FLIM-based atherosclerosis and cancer diagnosis.

In summary, we present a fully automated deconvolution method for TRFS data analysis based on an iterative Laguerre expansion approach. The method has been extensively validated on TRFS data from fluorescence lifetime standards, tissue constituents, and human tissue ex vivo. The main advantage of the proposed method is that it does not require any user intervention for tuning up the deconvolution process. We believe that our method will facilitate the use of TRFS in applications where online data analysis is required.