Simplified approach to least-square fitting of fluorescence lifetime ophthalmoscopy (FLIO) data by fixating lifetimes

: Fluorescence lifetime imaging ophthalmoscopy (FLIO) is a new imaging modality in ophthalmology. For clinical investigations, the amplitude-weighted mean of two or three lifetime components is usually analyzed. In this study, we investigated the eﬀects of ﬁxation of lifetime components. This resulted in slightly higher ﬁt errors but mean lifetimes were highly correlated to those from ﬁts with variable individual lifetimes. Furthermore, this approach resulted in a similarly good discrimination of diabetic retinopathy patients from controls, a reduction of the computational workload, a de-noising of the mean lifetime images and allows higher local resolution. Thus, ﬁxation of lifetimes in the ﬁt of FLIO data could be superior for clinical routine analysis of FLIO data.


Introduction
Fluorescence lifetime imaging ophthalmoscopy (FLIO) is a new imaging modality in Ophthalmology [1,2]. Currently, five prototype devices worldwide are in use for experimental diagnostic trials in various retinal diseases [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] as well as in animal studies [18,19]. With FLIO, it is possible to investigate intrinsic fluorophores of the retina and the retinal pigment epithelium (RPE) by analyzing their fluorescence lifetimes, i.e. the mean time, which fluorescent molecules remain in an excited electronic state. Fluorescence lifetimes are characteristic for single fluorophores as well as for their embedding matrix. Fluorophores, known to contribute to the autofluorescence of the ocular fundus, are lipids, lipoproteins and retinal-derived compounds of the age-pigment lipofuscin, collagens, elastin, glycated proteins, FAD, and NADH as minor contributor [20]. Furthermore, retinol fluorescence was investigated by Sharma et al. by two-photon imaging [21]. As, however, retinol and NADH have their excitation maximum in the UV, their contribution to FLIO might be limited [22]. The mixture of fluorophores in the retinal tissue results in a multi-exponential fluorescence decay over time. This can be approximated by a least-square fit of a series of exponential functions, representing individual lifetime components according to where I(t) is the fluorescence intensity at time t, I 0 is the intensity at t = 0, τ i are the time constants of fluorescence decay, a i describe their abundance, and n is the number of exponentials. IRF describes the temporal instrument response function and the asterisk denotes a convolution integral. A least-square fit procedure, seeking for the minimum of the goodness-of-fit parameter χ 2 , is typically used for the approximation of experimental decay data by adapting a i and τ i [23].
In clinical FLIO, the ocular fundus is scanned by a pulsed laser for fluorescence excitation and the emitted photons are recorded by time-correlated single photon counting (TCSPC). The number of photons is limited by the laser power that can be applied to the eye without hazard risk, the aperture of the pupil of the eye, the quantum yield of naturally occurring fluorophores, and the total investigation time, which should be convenient for the patient (typically 2-3 minutes). Thus, typical FLIO recordings have a signal-to-noise ratio (SNR) which allows for determining no more than three exponentials. For clinical routine investigations, a mean lifetime τ m is calculated as the amplitude weighted mean of single decay times (2) Although FLIO, this way, does not identify single retinal-derived or other flurophores, τ m is a useful additional marker for clinical diagnostics of age-related macular degeneration [4,6,11,17], diabetic retinopathy [9,24], Stargard's disease [14], macular teleangiectasia type 2 [5], retinal artery occlusion [16], central serous chorioretinopathy [10], choroideremia [15], retinitis pigmentosa [13,25], and Alzheimer's disease [26]. Pathologic changes in fluorescence lifetime may be due to a change in the abundance of the fluorophores mentioned above, their chemical alteration, or a change in their cellular micro-environment. However, we have to consider more fluorophores in the retina than the number of lifetime components that we can obtain by approximating the decay data. Thus, τ m is used to describe the overall fluorescence decay. But determining up to 6 parameters (a 1...3 and τ 1...3 ) from a fit of the measured decays by the model, given in Eq. 1, and subsequently averaging them according to Eq. 2, has several disadvantages: First, a sufficient SNR is needed to determine 6 parameters with sufficient accuracy. This is achieved by averaging the signal over several pixels (so called pixel binning) which, however, results in a loss of lateral resolution. Second, there is considerable correlation between the model parameters a 1...3 and τ 1...3 . Thus, they cannot be assumed as independent or even orthogonal. Third, fitting the model in Eq. 1 to the data needs a non-linear procedure which is time-consuming. In order to overcome these shortcomings, here we present a simplified model of fluorescence decay. This model assumes three fluorescence components with fixed lifetime components but variable abundance, i.e. in the fit of Eq. 1 to the data only a 1...3 will be modified. The results of this model were compared to a standard fit with the variables a 1...3 and τ 1...3 .

Subjects
In this study, we re-analyzed data from a previous study in addition to newly acquired data of patients with diabetic retinopathy and healthy controls [9]. All investigations in humans were carried out after informed consent was obtained from the subjects. The study adhered to the tenets of the declaration of Helsinki and was approved by the local ethics committee. Reference values for the lifetime components τ 1...3 were determined from 92 healthy control subjects (control cohort 1) in the age range of 19 through 81 years (mean: 46.0 ± 22.2 years). In the patient group, 34 subjects (mean age: 63.7 ± 8.5 years, range: 46-79 years) were included. Nine patients had mild and nine patients had moderate non-proliferative retinopathy, and 16 patients presented with severe non-proliferative retinopathy. A control group of 48 subjects (control cohort 2) was selected out of control cohort 1with a mean age of 66.2 ± 8.7 years (range: 41-81 years), which was not significantly different from that of the patients (p = 0.190). As lifetimes are known to be slightly different in phakic and pseudophakic subjects, only phakic subjects were included.

FLIO procedure and data processing
FLIO is described in detail elsewhere [8,12,20]. Briefly, a picosecond laser diode with a repetition rate of 80 Mhz was coupled into a laser scanning ophthalmoscope (Spectralis, Heidelberg Engineering, Heidelberg, Germany) for fluorescence excitation at 473 nm. The maximum radiant power through the pupil for a pulsed exposure for FLIO is 200 µW. Therefore, FLIO is using an average laser power well below the limit of 2.8 mW for the specified wavelength, recommended by the American National Standards Institute (ANSI Z136.   [27]. A detailed consideration of laser safety aspects of the device is given elsewhere [8]. Fluorescence photons were detected by TCSPC (SPC-150, Becker&Hickl GmbH, Berlin, Germany) in a short-wavelength (SSC: 498-560 nm) and a long-wavelength spectral channel (LSC: 560-720 nm). The detector count rates are monitored and the laser power can be reduced in order to avoid detector overflow. FLIO provides 30-degree images with a frame rate of 9 frames per second and a resolution of 256 × 256 pixels. Photon histograms over time, describing the fluorescence decay, were least-square fitted with a series of three exponential functions (Eq. 1) using the software SPCImage 5.1 (Becker&Hickl GmbH, Berlin, Germany). This way, a 1...3 and τ 1...3 were determined for all patients and controls (fit 1). In a second step, the lifetime components τ 1...3 were fixed as determined from the control cohort 1 and the fit according to Eq. 1 was repeated (fit 2). For both fits, τ m was calculated from Eq. 2 and averaged over the central area, the inner and the outer ring of the standardized Early Treatment of Diabetic Retinopathy Study (ETDRS) -grid ( Fig. 1) using the software FLIMX, which is documented and freely available for download online under the open source BSD license ((http://www.flimx.de). In the same way, a 1...3 and τ 1...3 from fit 1 as well as a 1...3 from fit 2 and the respective χ 2 values were averaged over the grid areas. All these procedures were performed for both spectral channels independently.

Statistical analysis
SPSS 24 (SPSS, Inc., Chicago, IL, USA) was used throughout this study, a p < 0.05 was considered significant. First, all samples of parameters, measured at the central area, inner, and outer ring of the ETDRS-grid as well as for the cohorts of controls and patients were checked for normal distribution using Kolmogorov-Smirnov test. As most of the parameters did not show normally distribution, averages were given as median values and non-parametric tests were used. Respective parameters were compared between fit 1 and fit 2 using the Wilcoxon-test. To compare results for patients and controls, Mann-Whitney U-test was used. Pearson's correlation coefficient was used to describe the mutual dependence of parameters. As a general measure of the interaction of all parameters per spectral channel, Cronbach's α was determined as from the number of variables N and the mean of their Pearson correlation coefficientsr. The discrimination of patient and age-matched controls (control cohort 2) by single parameters from fit 1 and fit 2 was investigated by generating the receiver-operating characteristics (ROC) and described by the area-under-curve (AUC). Finally, a generalized linear mixed model with a logit link function with random intercept was established for each of the fit procedures. The predicted probability of diabetic retinopathy was calculated from sufficiently independent fit parameters. Therefore, parameters with an AUC over 0.5 were considered. The Pearson correlation coefficient of these parameters was determined pairwise. If this was above 0.8, only the parameter with the higher AUC value was included. Subsequently, in a ROC analysis for this predicted probability the AUC as well as sensitivity and specificity were determined.

Analysis of control cohort 1
Fit 1 of these 92 healthy subjects revealed a clear and highly significant age-dependence of lifetimes (Table 1). The correlation of lifetime component τ 1 with age is shown as an example in Fig. 2. Thus, it was not possible to assign fix values to the lifetimes for fit 2. These had to be calculated from a linear regression with age according to equations given in Table 2. Using these lifetimes and fixing them in fit 2 resulted in qualitatively similar results than fit 1. Figure 3 (top) shows the τ m from both fits for the SSC of a healthy subject as a typical example.
Qualitatively, however, there were significant differences between the τ m from fits 1 and 2 for all areas of the ETDRS-grid except for the outer ring in the LSC. Data is given as median as well Table 1. Determination coefficient R 2 and significance p for the correlation of lifetimes from fit 1, averaged over the whole image excluding the optic disk, with subject age    Table 6 in the appendix. Furthermore, χ 2 increased by fixation of the decay times τ 1...3 (Fig. 3 bottom,  Table 6). Cronbach's α, given in Table 3, indicates that parameters of fit 2 are (except for the central area in the SSC) somewhat more independent from each other than that from fit 1. Table 3. Cronbach's a describing the dependence of parameters a 1...3 and τ 1...3 from fit 1

Discrimination of patients with diabetic retinopathy from healthy controls
Patients suffering from diabetic retinopathy showed generally longer fluorescence lifetimes than the subjects from the control cohort 2 with the same age. This holds for τ m from both fit procedures for all investigated areas and in both spectral channels ( Table 4). The Mann-Whitney U-test showed significant differences between patients and controls for both fits and all areas in the SSC as well as for fit 1 in the LSC. τ m from fit 2 in the LSC were significantly different only in the central area. The AUC values of the ROC curve show a somewhat better discrimination between patients and controls by τ m in fit model 1 for all regions and both spectral channels (Table 5). Four parameter from fit 1 met our criteria for inclusion into the predicted probability of diabetic retinopathy in a generalized linear mixed model (AUC > 0.5, Pearson's coefficient for the correlation with other parameters <0.8): a 2 (AUC = 0.670) as well as a 3 (AUC = 0.822) from the central area in the SSC, a 1 of the outer ring (AUC = 0.556), and a 3 from the central area (AUC = 0.661) in the LSC. The ROC-analysis with the predicted probability, established from these fit parameters, resulted in an AUC value of 0.865, a sensitivity of 89.53%, and a specificity of 73.53% at the optimal cut-off according to the Youden-index. From fit 2, we got a superior discrimination (AUC = 0.98, sensitivity = 91.17%, specificity = 97.92%) by the use of three fit parameters only: a 2 of the central area (AUC = 0.792) and a 3 of the outer ring (AUC = 0.903) from the SSC as well as a 1 of the outer ring from the LSC (AUC = 0.601).

Lifetime image quality and local resolution
It turned out that the images of lifetime components from fit 1 were much noisier than that from fit 2. Figure 4 gives an example. This noise is much less pronounced in the τ m images, but still persists as shown in Fig. 5. Furthermore, fitting three free parameters per spectral channel (fit 2) instead of 6 (fit 1) reduces the requirements with respect to SNR of the primary data, i.e. the decay curves per pixel. Thus, for fit 2 it is possible to reduce the binning from 5 × 5 to 3 × 3 pixels. This results in sharper τ m images as shown in Fig. 5 as well.

Discussion
In this work, we explored the possibility of fitting ocular fundus fluorescence decays, measured by FLIO, with three fixed lifetimes. This approach was already pursued by Murashova et al. for slices of murine retina, RPE, choroid, and sclera using lifetimes from the literature [28]. From a theoretical point of view, this approach assumes a mixture of three fluorophores which may vary in their concentration only. Obviously, this is an oversimplification. This also follows from the age dependence of these lifetimes, found here, which is in agreement with Dysli at el. [12]. On the other hand, the assumption of just three fluorescence components in the complex anatomical structure of the ocular fundus, as is made by the conventional FLIO fits with three variable lifetimes per spectral channel (fit 1), is also a simplified model as we know about a variety of fluorophores within the retina. Eldred and Katz isolated ten different fluorescent compounds only from the age pigment lipofuscin [29] and Schweitzer et al. measured the lifetimes of nine other naturally occurring fluorophores [20]. However, the SNR of FLIO measurements does not allow for fitting of more than three components. Furthermore, the high values, found for Cronbach's α (Table 5), show a strong interaction of all fluorescence parameters, obtained from a three-exponential fit. This vindicates the use of the amplitude-weighted mean τ m value for a generalized description of the fluorescence decay as was done in clinical studies. This is also possible if fixed lifetime components (as in fit 2) are used. Although this approach results in slightly different τ m than fit 1, the values show a very high correlation (R = 0.936 to 0.987, Table 6) and, thus, seem to give the same clinical information. This is shown here by the discrimination of diabetic retinopathy patients and controls.
However, higher χ 2 values for fit 2 compared to fit 1 indicate that the goodness of fit is inferior when fixed lifetime components are used. This is logically consistent as every additional fit parameter allows for a better approximation of the measurement. The χ 2 values ≤1.33, which we found for the inner and outer ring, however, indicate still sufficiently good fits and, thus, fit 2 can be considered as an appropriate model for our data [23]. This statement, however, has to be restricted to the extra-foveal retina/RPE. Since lifetimes, calculated according to the equations in Table 2, are determined from measurements in healthy subjects averaged over a 30°field of the retina excluding the optic disk, they result in an underestimation of the influence of the macular pigment, which is known to have very short lifetimes [3,8]. This results in considerably longer mean lifetimes τ m and higher χ 2 values of fit 2 compared to fit 1 in the central grid areas (Fig. 3). For the same reason, fit 2 may go wrong at the optic disc where very long lifetimes are abundant [20].
Cronbach's α showed only a marginal reduction upon fixation of the lifetime components (Table 3). This might seem surprising at a first glance. However, as a 3 , the abundance of the longest decay component τ 3 , turned out to be below 2% (Table 6), a 1 and a 2 must add up to almost 100% and, thus, must be highly dependent, resulting in high values of Cronbach's α.
Patients suffering from diabetic retinopathy have significantly longer lifetimes τ m from both fits in the SSC compared to subjects from a control cohort with the same mean age. Whereas fit 1 showed significantly longer τ m for the LSC as well, in fit 2 significance was found for the macula only. Nevertheless, generally both fits confirmed the findings by Schmidt et al. [9]. The discrimination of patients and controls was slightly better for τ m from fit 1 ( Table 5). Using the predicted probability of diabetic retinopathy, constructed from abundance values a1. . . 3 in a generalized linear mixed model, fit 2 resulted in higher AUC, sensitivity, and specificity values. The reason might be that a1. . . 3 are somewhat more independent in fit 2. Although, with respect to patient discrimination, a higher AUC for fitting with fixed lifetimes does not prove its superiority, it indicates that this approach can be considered as not inferior to fitting with variable lifetimes. It, however, has to be pointed out that FLIO, basically, is an imaging modality showing disease-specific molecular signature of fluorophores and their environment. It is not intended for use as an automated discrimination or diagnostic tool.
An advantage of the use of fixed lifetimes (fit 2) is that it needs a linear approximation only in contrast to non-linear procedure to fit exponential functions to the decay (fit 1). This decreases the computational load dramatically and yields in an immediate availability of the result. The greatest advantage, however, seems to be the conspicuous reduction of noise in the images. This is best seen in the abundance a 2 as demonstrated in Fig. 4. This noise reduction is in agreement with Walsh et al. [30], who found a dramatic reduction of the error of the abundance value upon fixation of the lifetime components. As FLIO is an imaging technique, noise reduction is essential in order to make otherwise hidden structures visible. Although it has to be pointed out that this noise reduction is less pronounced in the τ m images, it is still sufficient to reduce the pixel binning resulting in improved local resolution (Fig. 5). This can help in detecting small details in the images, which possibly indicate pathologic alterations, and thus, may improve the diagnostic capabilities of FLIO. The reason for the lower noise probably is the structure of the χ 2 error plane. For six free parameters in the fit (a 1...3 and τ 1...3 ), this may have multiple local minima. According to Knutson et al., already a fit with 2 free lifetimes may result in a valley in the χ 2 plane which makes it difficult to find the global minimum of the square sum of errors [31].
In real world noisy data, the least-square fit can easily stuck in a local χ 2 minimum, which can be a different one for neighboring pixels and, thus, result in noise of the fit parameters (Fig. 4) and, finally, τ m . Fixating the lifetime components, however, greatly reduces the complexity of the χ 2 error plane and the fit algorithm is able to find the global minimum of the fit error. This is reflected by less noisy τ m images. Knutson et al. suggested a global analysis of decays, measured in two spectral channels, using identical lifetime components in both channels to overcome the problem. [31] We tried this approach on our data by fitting with the mean values of fixed parameters, determined for SSC and LSC (data not shown). This, however, resulted in unacceptably high χ 2 values and a poor discrimination of diabetic patients and controls. Thus, this approach was not further pursued.
The following limitations of our study have to be mentioned: First, we tested the discrimination of patients and healthy subjects for the specific cohorts, included here, only. Thus, we cannot say how discrimination of both fit procedures work for other pathologies. The non-inferiority of fit 2, using fixed lifetime components, cannot be generalized without further studies. Second, we did not repeated measurements on the same subjects. Thus, we cannot comment on the reproducibility of parameters from either fit. Finally, the confocality of FLIO greatly reduces the contribution of lens fluorescence to measurements at the ocular fundus. However, as the lens of the human eye shows a very strong fluorescence, an influence of the lens on our measurements cannot be ruled out completely. Furthermore, fluorescence lifetimes of lens may be age dependent and, particularly, may change with beginning cataracts. As neither patients with severe cataracts nor pseudophakic patients were included in this study, equations in Table 2 may be different for cataract or pseudophakic patients.
In conclusion, we found that it is possible to fit FLIO data with fixed lifetime components if accounted for age-dependence. Although this approach results in a somewhat reduced goodness of fit based on a higher χ 2 , it gives τ m images which are qualitatively similar to that of fits with free individual lifetime components. Quantitatively, a high correlation between τ m from both fit methods was found. Furthermore, a similarly good discrimination between patients suffering from diabetic retinopathy and controls was found for both methods. The fit with fixed lifetime components results in a reduction of computational effort as well as a noise reduction in the τ m images and allows a higher spatial resolution. As a reduced number of parameters in the fit also reduces the requirements at the SNR, it seems possible to acquire FLIO images with shorter exposure times. A reduced imaging time would be advantageous for patient comfort as well as application of the technique in clinical routine. This however has to be proven in further studies using different exposure times in one cohort of subjects. Taken together, our results indicate that fit of FLIO data with fixed lifetime components should be considered for further testing in subsequent studies.

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