Translator Disclaimer
1 March 2009 Fully automated deconvolution method for on-line analysis of time-resolved fluorescence spectroscopy data based on an iterative Laguerre expansion technique
Author Affiliations +
Abstract
Time-resolved fluorescence spectroscopy (TRFS) is a powerful analytical tool for quantifying the biochemical composition of organic and inorganic materials. The potential of TRFS for tissue diagnosis has been recently demonstrated. To facilitate the translation of TRFS to the clinical arena, algorithms for online TRFS data analysis are essential. A fast model-free TRFS deconvolution algorithm based on the Laguerre expansion method has previously been introduced. One limitation of this method, however, is the need to heuristically select two parameters that are crucial for the accurate estimation of the fluorescence decay: the Laguerre parameter α and the expansion order. Here, a new implementation of the Laguerre deconvolution method is introduced, in which a nonlinear least-square optimization of the Laguerre parameter α is performed, and the optimal expansion order is selected based on a minimum description length criterion (MDL). In addition, estimation of the zero-time delay between the recorded instrument response and fluorescence decay is also performed based on normalized mean square error criterion (NMSE). The method is validated on experimental data from fluorescence lifetime standards, endogenous tissue fluorophores, and human tissue. The proposed automated Laguerre deconvolution method will facilitate online applications of TRFS, such as real-time clinical tissue diagnosis.

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.3 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).13 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.3 In the particular case of tissue TRFS data, fluorescence emission typically originates from several endogenous fluorophores, each one potentially presenting complex decay dynamics.3 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.3 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.16 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,15 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 α to guaranty accurate estimations of the fluorescence IRF.14, 15 The Laguerre parameter (0<α<1) determines the rate of exponential decline of the DLFs.16 Thus, small α values are usually needed for expanding short-lived fluorescence decays, while large α 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 α 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 α , in which simultaneous nonlinear optimization of both α and the Laguerre expansion coefficients is performed.14 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.15 Thus, we propose a more computationally efficient approach for optimizing the Laguerre parameter, in which nonlinear optimization is only applied to estimate α , 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.19 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 α , 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(n) is given by the convolution of the fluorescence IRF h(n) with the instrument response x(n) :

Eq. 1

y(n)=Tm=0K1h(m)x(nm),n=0,,N1.
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(n) and x(n) . The Laguerre deconvolution technique expands the fluorescence IRF on an orthonormal set of discrete time Laguerre functions (DLF) bjα(n) :

Eq. 2

h(n)=T.j=0L1cjαbjα(n).
In Eq. 3, cjα are the unknown Laguerre expansion coefficients (LEC), which are to be estimated from the input-output data; bjα(n) denotes the j ’th order orthonormal DLF;16 L is the number of DLFs used to model the IRF, thus defining the order of the expansion. The Laguerre parameter α (0<α<1) determines the rate of asymptotic decline of the DLFs. Thus, larger α values imply longer convergence time to zero. By inserting Eq. 2 into Eq. 1, the convolution Eq. 1 becomes:

3.

y(n)=j=0L1cjαvjα(n)
vjα(n)=i=0K1bjα(m)x(nm).
The functions vjα(n) , representing the digital convolution of the input x(n) 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

024030_1_039902jbod1.jpg
Finally, the expansion coefficients ĉα can be estimated by least-square estimation as:

Eq. 5

ĉα=(VαTVα)1VαTy¯.
Once the expansion coefficient ĉα has been estimated, the fluorescence IRF h(n) can be computed from Eq. 2 and the average lifetimes can be calculated as:3

Eq. 6

τ=T.n=0Nn.h(n)n=0Nh(n).

2.1.2.

Nonlinear least square optimization of the Laguerre parameter

The choice of the Laguerre parameter α is critical in achieving accurate fluorescence IRF expansions. As a first approach, the parameter α 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 α 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 ( α in this case) is chosen so that an objective function is minimized. If we write Eq. 4 as:

Eq. 7

y¯=Vαc¯α,
then we know that the least-Squares solution of this equation is given by:

Eq. 8

ŷ=Vαĉα=Vα(VαTVα)1VαTy¯.
We now define the error in our estimation as:

Eq. 9

εN=y¯ŷ=y¯Vαĉα=y¯Vα(VαTVα)1VαTy¯,
hence our cost function can be defined as:

Eq. 10

F(α)=εNTεN.
Our task is to find a value of α that minimizes the prior cost function:

Eq. 11

α̂=argminαF(α).
Since the cost function is a nonlinear function of α , minimization requires an iterative scheme for determining α̂ . Thus, an iterative scheme for α is defined as:

Eq. 12

α̂(k+1)=α̂(k)γ(k)[dF(α)dα](k),
where γ(k) is a scaling factor chosen in a way such that the value of α̂ is incremented by a fixed amount after each iteration. Note that the superscript k denotes the iteration number. Estimation of the Jacobian dF(α)dα can be achieved by numerical methods. Alternatively, closed form recursive formulation of dF(α)dα can also be derived as follows.

Substituting εN from Eq. 9 into Eq. 10, we get:

Eq. 13

F(α)=[y¯Vα(VαTVα)1VαTy¯]T[y¯Vα(VαTVα)1VαTy¯].
This can be further simplified to:

Eq. 14

F(α)=y¯T[INVα(VαTVα)1VαT]y¯=y¯T[y¯Vα(VαTVα)1VαTy¯],
where IN is the identity matrix of order N . By identifying the second term on the RHS as εN , it can be further simplified to:

Eq. 15

F(α)=y¯TεN.
Thus, the Jacobian can be simply expressed as:

Eq. 16

dF(α)dα=y¯Td(εN)dα.
Using the definition of εN as in Eq. 9, its derivative with respect to α becomes:

Eq. 17

d(εN)dα=d(y¯Vαĉα)dα=d(Vαĉα)dα=[d(Vα)dαĉα+Vαd(ĉα)dα].
The key variables vjα(n) forming the matrix Vα can be computed recursively using the following formulation:16

18.

vjα(n)=αvjα(n1)+αvj1α(n)vj1α(n1),
v0α(n)=αv0α(n1)+1αx(n).
From here, we can derive a recursive formulation for the terms in d(Vα)dα :

19.

dvjα(n)dα=12α[vjα(n1)+vj1α(n)]+α[dvjα(n1)dα+dvj1α(n)dα]dvj1α(n1)dα,
dv0α(n)dα=12α[v0α(n1)]+α[dv0α(n1)dα]+121αx(n).
Similarly, by taking the derivative with respect to α to the least-square solution of the expansion coefficients as defined in Eq. 5, we obtain:

Eq. 20

dĉαdα=d[(VαTVα)1VαTy¯]dα ={d[(VαTVα)1]dαVαTy¯+(VαTVα)1d[VαTy¯]dα}.
To calculate d[(VαTVα)1]dα , let us first consider that:
d[(VαTVα)1(VαTVα)]dα=d[(VαTVα)1]dα(VαTVα)+(VαTVα)1d[(VαTVα)]dα.
Since (VαTVα)1(VαTVα)=IN , the LHS becomes zero, thus:

Eq. 21

d[(VαTVα)1]dα=(VαTVα)1d[(VαTVα)]dα(VαTVα)1.
Substituting Eqs. 5, 21 into Eq. 20, we get:

Eq. 22

dĉαdα=(VαTVα)1{d(VαT)dαy¯[d(VαT)dαVα+VαTd(Vα)dα]ĉα}.
By inserting Eq. 22 into Eq. 17, the Jacobian can be expressed as:

Eq. 23

dF(α)dα=y¯T(d(Vα)dαĉα+Vα(VαTVα)1{d(VαT)dαy¯[d(VαT)dαVα+VαTd(Vα)dα]ĉα}).
Finally, our last iterative relation for calculating α̂ becomes:

Eq. 24

α̂(k+1)=α̂(k)γ(k)y¯T(d(Vα)dαĉα+Vα(VαTVα)1{d(VαT)dαy¯[d(VαT)dαVα+VαTd(Vα)dα]ĉα})(k),
where the terms of d(Vα)dα and d(VαT)dα can be calculated using the recursive relations in Eq. 19. Iterative updating of the Laguerre parameter α̂ is performed until the cost function F(α) 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),20 which is defined as:

Eq. 25

MDL(L)=F(α)Nexp[ln(N)LN].
Here, N is the number of measured samples, F(α) 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 MDL(L) is computed. The value of L that yields the minimum 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.19 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.

2.2.

Method Validation on Experimental Time-Resolved Fluorescence Spectroscopy Data

The performance of the automated Laguerre deconvolution technique was assessed with experimental TRFS data from: 1. fluorescence lifetime standards, 2. a set of well characterized tissue endogenous fluorphores, and 3. human tissue. A brief description of the TRFS measurement is presented here. The method was validated in terms of lifetime values retrieved by our method, and their agreement with the corresponding values obtained using the LSIR approach and the ones reported in the literature. In addition, the estimated residuals normalized with respect to the peak amplitude of the measured fluorescence decay were observed to assess estimation accuracy. Finally, the whiteness of the residuals was also evaluated by analyzing the autocorrelation function of the residuals [i.e., checking whether or not its values were within the 95% confidence level for random correlation, computed as ±1.96*C(0)N , where C(0) is the autocorrelation values for lag zero].21

2.2.1.

Instrumentation for time-resolved fluorescence spectroscopy measurements

The experiments were conducted with a TRFS system, built up on a first prototype,22 currently being evaluated for clinical diagnostic of cancer and cardiovascular diseases.8, 9, 10, 11 Briefly, the autofluorescence of the sample was induced with a Q-switched UV laser ( 355nm , 1-ns pulse width, 10-kHz repetition rate) using a bifurcated fiber optic probe. The collected autofluorescence was dispersed by a monochromator (MicroHR, Horiba), and detected with a gated multichannel plate photomultiplier tube (transit-time spread of 90ps ). The autofluorescence was temporally resolved using a digital oscilloscope (bandwidth 2.5GHz , sampling rate 10Gsampless ) coupled to a preamplifier (bandwidth 1.5GHz ).

2.2.2.

Fluorescence lifetime standards and tissue endogenous fluorphores

The fluorescence lifetime standards used for our method validation were selected to cover a broad range of emission wavelengths (360to650nm) and radiative lifetimes (0.4to12ns) , and included rose Bengal (33,000, Sigma-Aldrich, Saint Louis, Missouri), rhodamin B (25,242, Sigma-Aldrich), and 9-cyanoanthracene (15,276, Sigma-Aldrich) dissolved in either ethanol, methanol, or distilled water to a concentration of 106M . In addition, our deconvolution method was also validated in three tissue endogenous fluorophores, relevant for potential fluorescence-based tissue diagnosis:1, 2 collagen type 1 (C3511, Sigma-Aldrich), NADH (N8129, Sigma-Aldrich), and flavin adenine dinucleotide (FAD) disodium salt dihydrate (F6625, Sigma-Aldrich). The collagen was tested in dry form, while NADH and FAD were tested in 106-M phosphate buffered saline (PBS) solutions. For each sample, the fluorescence decay corresponding to the peak emission wavelength of that fluorophore was recorded. After the fluorescence measurement, the scattered laser pulse temporal profile was measured at a wavelength slightly below the excitation laser line. The laser pulse energy at the tip of the excitation fiber probe was adjusted to 5μJpulse .

2.2.3.

Tissue autofluorescence

Finally, the automated Laguerre deconvolution method was also tested on TRFS measurements obtained from human tissue. We used five human postmortem coronary arteries for this validation study, taking advantage of our ongoing research on optical diagnosis of atherosclerosis. These five specimens were a combination of normal arteries and atherosclerotic arteries of various types. In this case, the time-resolved fluorescence spectra were measured from the lumen side of the arteries, covering the spectral range between 400to600nm . The scattered laser pulse temporal profile was also measured right after each time-resolved fluorescence spectrum acquisition. The laser pulse energy at the tip of the excitation fiber probe was also adjusted to 5μJpulse .

2.2.4.

Comparison with standard deconvolution method

To compare the deconvolution performance of our method with more standard approaches, we also analyzed the same TRFS data using the classical multiexponential least-square iterative reconvolution (LSIR) approach.3 LSIR applies nonlinear least-square optimization methods to estimate the parameters of a multiexponential IRF that would fit best its convolution with the instrument response with the fluorescence decay data. In this study, one or two exponential components were considered to model the fluorescence IRFs.

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 α=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 α=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 α=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.

024030_1_039902jbo1.jpg

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.

SampleSolventWavelength(nm)OptimalalphaOptimalorder Lifetime (ns) Computation time (ms)
Opt. LAGLSIRLiteratureLAGOpt. LAGLSIR
9-CAEthanol450 0.96±0.02 4 12.61±0.58 12.18±0.73 11.7 to 12.287266153
Rhodamin BEthanol580 0.84±0.01 4 2.52±0.10 2.44±0.14 2.60 to 3.01525867
Rhodamin B H2O 580 0.74±0.01 4 1.55±0.03 1.58±0.05 1.48 to 1.67519368
Rose BengalEthanol570 0.62±0.01 3 0.65±0.02 0.61±0.01 0.8436945
Rose BengalMethanol570 0.57±0.01 3 0.48±0.01 0.48±0.01 0.5433273
CollagenN/A400 0.91±0.02 3 3.49±0.14 3.61±0.15 2 to 43175366
NADHPBS450 0.67±0.01 4 0.45±0.03 0.36±0.01 0.30 to 0.50454139
FADPBS530 0.93±0.02 3 3.11±0.13 3.17±0.12 3 to 43157779

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 400nm , 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 α=0.91 and order L=3 were obtained, and two exponential components were needed for proper multiexponential deconvolution ( tau1=4.82 , tau2=1.21 , a1a2=0.34 ). The NADH fluorescence IRF and lifetime value at the peak emission of 450nm were also properly estimated by both methods, as reflected in Figs. 2c and 2d, and Table 1. For this fluorophore, optimal values of α=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 530nm . Optimal values for α=0.93 and order L=3 were obtained, and two exponential terms were used ( tau1=3.32 , tau2=0.81 , a1a2=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.

024030_1_039902jbo2.jpg

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.8ns , 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-ns delay.

024030_1_039902jbo3.jpg

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 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 1.5to3ns , 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.

024030_1_039902jbo4.jpg

In addition, we wanted to assess ranges of values of the parameter α , 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 (mean±SE) as a function of wavelengths from all datasets collected. Results of this analysis are shown in Fig. 5 . It was observed that optimal α values for accurate deconvolution of arterial TRFS data were between 0.90 to 0.94, and there was no obvious variation between α 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 0.5ns at 400nmto0.95ns at 550nm , 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 (NMSE<2%) [Fig. 5d].

Fig. 5

Average ( mean±SE , from 47 datasets) values for (a) the Laguerre parameter α , (b) expansion order, (c) zero-time delay, and (d) NMSE as a function of emission wavelengths from all human artery TRFS datasets. Values of α 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 500nm due to lower SNR at higher wavelengths.

024030_1_039902jbo5.jpg

3.3.

Computation Speed Assessment of the Automated Laguerre Deconvolution

To compare the computation speed of the proposed Laguerre deconvolution with those of both its original implementation (in which the Laguerre parameter is fixed a priori) and the standard LSIR method, the computation times for deconvolving each of the TRFS datasets from the fluorescence standards and endogenous fluorophores were reported, as shown in Table 1. As expected, the original Laguerre deconvolution implementation was 1 to 2 orders of magnitude faster than both the LSIR and the automated Laguerre deconvolution. However, the computation speed of the proposed automated Laguerre deconvolution was comparable to that of the multiexponential LSRI method, and it was even faster for those fluorophores showing more than one exponential decay component. The algorithms were implemented in Matlab (version R2007a, MathWorks, Natick, Massachusetts) and ran on a PC computer having a dual core 2.4-GHz processor.

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 α 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 α 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 α value, respectively. Thus, fast or slow decaying fluorescence IRFs can be optimally expanded on a set of DLFs defined by a small or large α parameter value, respectively. The proposed automated Laguerre deconvolution method performs a nonlinear least-square search for the optimal values of α 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 α 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 α 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 α until finding an optimal one, this approach is unsuitable for online applications.15 The second implementation, proposed by Maarek,14 provides a method for estimating both the parameter α and the IRF expansion; however, it is computationally very expensive (nonlinear optimization is applied not only for the α parameter but for the expansion coefficients, too), thus is not suitable for online analysis either. By separating the optimization of α 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.20 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.20 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.19 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 α were optimal for accurate deconvolution of human artery TRFS data. Our results showed that optimal values for α were concentrated within a range of 0.90 to 0.94 with no apparent correlation between α and the emission wavelength [Fig. 5a]. These findings indicate that it may be possible to speed up the optimization of the parameter α , 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 α (e.g., α=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 α , 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 α 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 publication15 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 α 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.25 As for the case of the original Laguerre TRFS deconvolution, the current Laguerre FLIM deconvolution also requires to select the expansion parameters (i.e., α 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.

Acknowledgments

This work was supported by the American Heart Association, Texas Affiliate, Beginning Grant-in-Aid grant 0765102Y and NIH grant 1-R21-CA132433.

References

1. 

R. Richards-Kortum and E. Sevick-Muraca, “Quantitative optical spectroscopy for tissue diagnosis,” Annu. Rev. Phys. Chem., 47 555 –606 (1996). https://doi.org/10.1146/annurev.physchem.47.1.555 0066-426X Google Scholar

2. 

N. Ramanujam, “Fluorescence spectroscopy data of neoplastic and non-neoplastic tissues,” Neoplasia, 2 (1–2), 89 –117 (2000). https://doi.org/10.1038/sj.neo.7900077 1522-8002 Google Scholar

3. 

J. R. Lakowicz, Principles of Fluorescence Spectroscopy, 3rd ed.Springer, Berlin (2006). Google Scholar

4. 

M. A. Mycek, K. T. Schomacker, and N. S. Nishioka, “Colonic polyp differentiation using time-resolved autofluorescnence spectroscopy,” Gastrointest. Endosc., 48 (4), 390 –394 (1998). https://doi.org/10.1016/S0016-5107(98)70009-4 0016-5107 Google Scholar

5. 

L. Marcu, M. C. Fishbein, J. M. Maarek, and W. S. Grundfest, “Discrimination of human coronary artery atherosclerotic lipid-rich lesions by time-resolved laser-induced fluorescence spectroscopy,” Arterioscler., Thromb., Vasc. Biol., 21 (7), 1244 –1250 (2001). https://doi.org/10.1161/hq0701.092091 1079-5642 Google Scholar

6. 

L. Marcu, J. A. Jo, P. V. Butte, W. H. Yong, B. K. Pikul, K. L. Black, and R. C. Thompson, “Fluorescence lifetime spectroscopy of glioblastoma multiforme,” Photochem. Photobiol., 80 (1), 98 –103 (2004). https://doi.org/10.1562/2003-12-09-RA-023.1 0031-8655 Google Scholar

7. 

P. V. Butte, B. K. Pikul, A. Hever, W. H. Yong, K. L. Black, and L. Marcu, “Diagnosis of meningioma by time-resolved fluorescence spectroscopy,” J. Biomed. Opt., 10 (6), 064026 (2005). https://doi.org/10.1117/1.2141624 1083-3668 Google Scholar

8. 

L. Marcu, Q. Y. Fang, J. A. Jo, T. Papaioannou, A. Dorafshar, T. Reil, J. H. Qiao, J. D. Baker, J. A. Freischlag, and M. C. Fishbein, “In vivo detection of macrophages in a rabbit atherosclerotic model by time-resolved laser-induced fluorescence spectroscopy,” Atherosclerosis, 181 (2), 295 –303 (2005). https://doi.org/10.1016/j.atherosclerosis.2005.02.010 0021-9150 Google Scholar

9. 

J. A. Jo, Q. Fang, T. Papaioannou, J. D. Baker, A. H. Dorafshar, T. Reil, J. H. Qiao, M. C. Fishbein, J. A. Freischlag, and L. Marcu, “Laguerre-based method for analysis of time-resolved fluorescence data: application to in vivo characterization and diagnosis of atherosclerotic lesions,” J. Biomed. Opt., 11 (2), 021004 (2006). https://doi.org/10.1117/1.2186045 1083-3668 Google Scholar

10. 

W. H. Yong, P. V. Butte, B. K. Pikul, J. A. Jo, Q. Fang, T. Papaioannou, K. Black, and L. Marcu, “Distinction of brain tissue, low grade and high grade glioma with time-resolved fluorescence spectroscopy,” Front. Biosci., 11 1255 –1263 (2006). https://doi.org/10.2741/1878 1093-4715 Google Scholar

11. 

L. Marcu, J. A. Jo, Q. Fang, T. Papaioannou, T. Reil, J.—H. Qiao, J. D. Baker, J. A. Freischlag, and M. C. Fishbein, “Detection of rupture-prone atherosclerotic plaques by time-resolved laser-induced fluorescence spectroscopy,” Atherosclerosis, (0021-9150) https://doi.org/10.1016/j.atherosclerosis.2008.08.035 Google Scholar

12. 

D. V. O'Connor, W. R. Ware, and J. C. Andre, “Deconvolution of fluorescence decay curves—critical comparison of techniques,” J. Phys. Chem., 83 (10), 1333 –1343 (1979). https://doi.org/10.1021/j100473a019 0022-3654 Google Scholar

13. 

W. R. Ware, L. J. Doemeny, and T. L. Nemzek, “Deconvolution of fluorescence and phosphorescence decay curves—least-squares method,” J. Phys. Chem., 77 (17), 2038 –2048 (1973). https://doi.org/10.1021/j100636a003 0022-3654 Google Scholar

14. 

J. M. I. Maarek, L. Marcu, W. J. Snyder, and W. S. Grundfest, “Time-resolved fluorescence spectra of arterial fluorescent compounds: reconstruction with the laguerre expansion technique,” Photochem. Photobiol., 71 (2), 178 –187 (2000). https://doi.org/10.1562/0031-8655(2000)071<0178:TRFSOA>2.0.CO;2 0031-8655 Google Scholar

15. 

J. A. Jo, Q. Fang, T. Papaionnou, and L. Marcu, “Fast model-free deconvolution of fluorescence decay for analysis of biological systems,” J. Biomed. Opt., 9 (4), 743 –752 (2004). https://doi.org/10.1117/1.1752919 1083-3668 Google Scholar

16. 

V. Z. Marmarelis, “Identification of nonlinear biological-systems using Laguerre expansions of kernels,” Ann. Biomed. Eng., 21 (6), 573 –589 (1993). https://doi.org/10.1007/BF02368639 0090-6964 Google Scholar

17. 

J. A. Jo, A. Blasi, E. M. Valladares, R. Juarez, A. Baydur, and M. C. Khoo, “A nonlinear model of cardiac autonomic control in obstructive sleep apnea syndrome,” Ann. Biomed. Eng., 35 (8), 1425 –1443 (2007). https://doi.org/10.1007/s10439-007-9299-5 0090-6964 Google Scholar

18. 

G. D. Mitsis, R. Zhang, B. D. Levine, and V. Z. Marmarelis, “Modeling of nonlinear physiological systems with fast and slow dynamics. II. Application to cerebral autoregulation,” Ann. Biomed. Eng., 30 (4), 555 –565 (2002). https://doi.org/10.1114/1.1477448 0090-6964 Google Scholar

19. 

M. Van Den Zegel, N. Boens, D. Daems, and F. C. De Schryver, “Possibilities and limitations of the time-correlated single photon counting technique: a comparative study of correction methods for the wavelength dependence of the instrument response function,” Chem. Phys., 101 (2), 311 –335 (1986). https://doi.org/10.1016/0301-0104(86)85096-0 0301-0104 Google Scholar

20. 

L. Ljung, System Identification, Theory for the User, 2nd ed.Prentice Hall PTR, Englewood Cliffs, NJ (1999). Google Scholar

21. 

D. R. James and W. R. Ware, “Recovery of underlying distributions of lifetimes from fluorescence decay data,” Chem. Phys. Lett., 126 (1), 7 –11 (1986). https://doi.org/10.1016/0009-2614(86)85107-7 0009-2614 Google Scholar

22. 

Q. Y. Fang, T. Papaioannou, J. A. Jo, R. Vaitha, K. Shastry, and L. Marcu, “Time-domain laser-induced fluorescence spectroscopy apparatus for clinical diagnostics,” Rev. Sci. Instrum., 75 (1), 151 –162 (2004). https://doi.org/10.1063/1.1634354 0034-6748 Google Scholar

23. 

N. Boens, W. Qin, N. Basaric, J. Hofkens, M. Ameloot, J. Pouget, J. P. Lefevre, B. Valeur, E. Gratton, M. vandeVen, N. D. Silva, Y. Engelborghs, K. Willaert, A. Sillen, G. Rumbles, D. Philips, A. J. W G. Visser, A. vanHoek, J. R. Lakowicz, H. Malak, I. Gryczynski, A. G. Szabo, D. T. Krajcarski, N. Tamai, and A. Miura, “Fluorescence lifetime standards for time and frequency domain fluorescence spectroscopy,” Anal. Chem., 79 (5), 2137 –2149 (2007). https://doi.org/10.1021/ac062160k 0003-2700 Google Scholar

24. 

D. P. Jonathan and M. Mary-Ann, “Design and development of a rapid acquisition laser-based fluorometer with simultaneous spectral and temporal resolution,” Rev. Sci. Instrum., 72 (7), 3061 –3072 (2001). https://doi.org/10.1063/1.1379957 0034-6748 Google Scholar

25. 

J. A. Jo, Q. Y. Fang, and L. Marcu, “Ultrafast method for the analysis of fluorescence lifetime imaging microscopy data based on the Laguerre expansion technique,” IEEE J. Sel. Top. Quantum Electron., 11 (4), 835 –845 (2005). https://doi.org/10.1109/JSTQE.2005.857685 1077-260X Google Scholar
©(2009) Society of Photo-Optical Instrumentation Engineers (SPIE)
Aditi S. Dabir, Chintan Trivedi, Yeontack Ryu, Paritosh Pande, and Javier A. Jo "Fully automated deconvolution method for on-line analysis of time-resolved fluorescence spectroscopy data based on an iterative Laguerre expansion technique," Journal of Biomedical Optics 14(2), 024030 (1 March 2009). https://doi.org/10.1117/1.3103342
Published: 1 March 2009
JOURNAL ARTICLE
13 PAGES


SHARE
Advertisement
Advertisement
Back to Top