Developing and Testing a Bayesian Analysis of Fluorescence Lifetime Measurements

FRET measurements can provide dynamic spatial information on length scales smaller than the diffraction limit of light. Several methods exist to measure FRET between fluorophores, including Fluorescence Lifetime Imaging Microscopy (FLIM), which relies on the reduction of fluorescence lifetime when a fluorophore is undergoing FRET. FLIM measurements take the form of histograms of photon arrival times, containing contributions from a mixed population of fluorophores both undergoing and not undergoing FRET, with the measured distribution being a mixture of exponentials of different lifetimes. Here, we present an analysis method based on Bayesian inference that rigorously takes into account several experimental complications. We test the precision and accuracy of our analysis on controlled experimental data and verify that we can faithfully extract model parameters, both in the low-photon and low-fraction regimes.


Introduction
Förster resonance energy transfer, or FRET, is a fluorescence technique commonly used to access spatial information on length scales smaller than the diffraction limit of light [1]. In standard fluorescence, illuminating light is used to excite a fluorophore into a higher energy state, and the fluorophore subsequently relaxes into its ground state either by emitting a photon or through a non-radiative decay pathway. If another fluorophore is near, typically within % 10 nm, the two fluorophores can interact through dipole-dipole interactions termed FRET. FRET confers an additional decay path where the excited florophore, termed the donor, can transfer its energy to the nearby, unexcited fluorophore, termed the acceptor, which can then release the energy as a photon or through non-radiative decay. As the emission spectra of commonly used donor and acceptor pairs are spectrally distinct, one common method of measuring the the average FRET efficiency is to compare the relative intensities collected from the two channels. However, this method has drawbacks including spectral bleed-through and a sensitivity to changes in fluorophore concentration and excitation light intensity [2]. As an alternative to using fluorescence intensity to quantify FRET, fluorescence lifetime imaging microscopy, or FLIM, can be used [3][4][5][6]. FLIM is a general technique that allows changes in a flurophore's local environment to be probed. While use of FLIM is not limited to measuring changes in FRET, it can be used in this context without some of the drawbacks of an intensity based measurement. In time-domain FLIM, a narrow pulse of light is used to excite fluorophores into an excited state. Fluorophores that decay from their excited states can do so by releasing a photon. A subset of the released photons are detected, and for each detected photon, the arrival time is measured relative to the excitation pulse. The amount of time fluorophores spend in their excited state depends on the number of decay paths available. Donor fluorophores are chosen such that when they decay from their excited states, they do so at a constant rate, leading to photon emission time distributions that are exponential with a single characteristic decay time. This characteristic decay time is known as the fluorescence lifetime and is typically on the order of nanosconds. When donor fluorophores are undergoing FRET, they will spend, on average, a shorter amount of time in their excited states, leading to a reduced lifetime and quantum efficiency [7]. In a sample where only a fraction of donor fluorophores are undergoing FRET, the photon emission time distribution will be the sum of two exponentials with different lifetimes. By comparing the amplitudes of these two exponentials, the relative fraction of donors undergoing FRET can be measured. In practice, additional complications are present, including photons collected from spurious background and time delays introduced by the collection system itself. These effects must be accounted for to infer the relative amplitudes and lifetimes of the emitted photon distributions from the measured photon arrival time histograms. Several approaches have been used in order to estimate these parameters, including least-squares fitting [8], rapid lifetime determination [9], phasor methods [10][11][12], and Bayesian approaches [13], each with their own advantages and disadvantages.
Here we utilize and extend the Bayesian approach previously described [13] to take into account biexponential decays and additional experimental factors and we test the performance of our method using experimental data.

Bayesian Framework
Our framework is based on a previously described Bayesian analysis approach for measuring lifetimes from FLIM data [13]. For an introductory overview of Bayesian analysis, we direct the reader to [14]. Bayes' Law states that given a set of data, t, and a set of model parameters θ, then, where p(θ|t), the probability of the model parameters given the measured data, is referred to as the posterior distribution, p(t|θ) is referred to as the likelihood function, and p(θ) is referred to as the prior distribution. The aim of Bayesian inference approaches is to find the posterior distribution for the given model and data, and hence what the probability is for each possible set of model parameters θ.
In time-domain FLIM measurements, a narrow laser pulse is used to excite fluorophores in the sample, and the arrival times of photons emitted from the fluorophores are recorded. Fluorophores undergoing FRET will have a shorter florescence lifetime compared with fluorophores not undergoing FRET. When only a fraction of fluorophores in the sample are undergoing FRET, the resulting distribution of photon emission will be a sum of exponentials, where each exponential has a different lifetime, and each exponential is weighted by the number of photons collected from the respective source. In addition, there exists a constant background of photons due to noise in the detector and stray light, taken to be from a uniform distribution. In the following, we consider photons from each of these sources separately and construct the likelihood function as follows, where t is the arrival time of a photon relative to the excitation pulse, τ S and τ L are respectively the short and long fluorescence lifetimes, f S and f L are the fractions of photons from the short and long lifetime distributions respectively, f B is the fraction of photons from the uniform background given by Here p i (t|f i ) is the probability of the photon arriving at time t given that the photon originates from fraction f i . Eq (2) represents the likelihood model when time is taken to be continuous. However, in practice, photon arrival times collected with TCSPC are discretized into bins, and this discretization must be taken into account. If the bins are numbered sequentially and of width Δt, such that b i represents the bin containing photons with arrival time, (i − 1)Δt t iΔt, then the likelihood function becomes, Thus, Eq (3) serves as the discrete form of the likelihood function, Eq (2).

Instrument Response Function
One complexity in experimental TCSPC measurements is that a delay is introduced to photon arrival times, termed the Instrument Response Function (IRF). In order to account for this effect, the IRF was experimentally measured (see FLIM Measurements). The measured IRF is then convolved with the idealized probability density functions for the exponential distributions in order to construct the likelihood function. Taking this effect into account leads to, where p em,j is the idealized exponential distribution, taken to be / e À t=t j , where j 2 {S, L} is an index labeling the exponential distribution and IRF(t) is the experimentally measured instrument response function.

Posterior Distribution
Using Eq (4) in Eq (3) leads to the final form of our likelihood function, For comparison with experiments using control dyes where the lifetimes of the two molecules are well characterized, we choose a prior distribution such that the distribution is uniform for the fractions in the domain f j 2 [0, 1], and τ S and τ L are set to the measured values for Coumarin 153 and Erythrosin B respectively. With this choice of prior, Eq (1) becomes, pðyjtÞ / pðtjyÞ ð6Þ and hence our posterior distribution is proportional to our likelihood function in the constrained parameter space. To build the posterior distribution, parameter space is searched by evaluating the likelihood function on a grid of uniform spacing. Alternatively, parameter space can be searched stochastically using the Markov chain Monte Carlo method, yielding equivalent results (S1 Fig).

Effects of Periodic Excitation
For a single exponential decay, the probability of measuring a photon at time t, given a decay lifetime, τ, is given by, Where τ is the lifetime of the fluorophore. In practice, many sequential excitation pulses are used, and it's possible that a fluorophore excited by a given pulse doesn't emit a photon until after a future pulse. Taking this effect into account for a single exponential decay leads to [13], where T is the excitation pulse period and k is an index counting previous pulses. The sum is a geometric series, which converges to, Thus, accounting for periodic excitations leads to a prefactor 1 1À e À T t , which for a given T and τ is constant. As we treat exponentials from populations with short and long lifetimes separately, this factor can safely be absorbed into the normalization constant, leaving the probability distribution unchanged.

FLIM Measurements
FLIM measurements were carried out on a Nikon Eclipse Ti microscope using two-photon excitation from a Ti:sapphire pulsed laser (Mai-Tai, Spectra-Physics, 865 nm or 950 nm wavelength, 80 MHz repetition rate, % 70 fs pulse width), a commercial scanning system (DCS-120, Becker & Hickl), and hybrid detectors (HPM-100-40, Becker & Hickl). The excitation laser was collimated by a telescope assembly to avoid power loss at the XY galvanometric mirror scanner and to fully utilize the numerical aperture of a water-immersion objective (CFI Apo 40x WI, NA 1.25, Nikon). Fluorescence was imaged with a non-descanned detection scheme with a dichroic mirror (705 LP, Semrock) that was used to allow the excitation laser beam to excite the sample while allowing fluorescent light to pass into the detector path. A short-pass filter was used to further block the excitation laser beam (720 SP, Semrock) followed by an emission filter appropriate for Coumarin and Erythrosin B (550/88nm BP, Semrock, or 552/ 27nm BP, Semrock). A Becker & Hickl Simple-Tau 150 FLIM system was used for time correlated single photon counting [15]. The instrument response function was acquired using second harmonic generation of a urea crystal [15].
For the data shown in Fig 1, the TAC range was set to 7 × 10 −8 with a Gain of 5, corresponding to a 14 ns maximum arrival time. The TAC offset was set to 6.27%. The TAC limit high and limit low were set to 5.88% and 77.25%, respectively, resulting in a 10 ns recording interval. Erythrosin B and Coumarin 153 samples were prepared at 10 mM and 15 mM, respectively. Lifetimes were measured and fixed at values of 3.921 ns and 0.453 ns for Coumarin 153 and Erythrosin B respectively.
For the data shown in Fig 2, the TAC range was set to 5 × 10 −8 with a Gain of 5, corresponding to a 10 ns maximum arrival time. The TAC limit high and limit low were set to 95.29% and 5.88%, respectively, resulting in a 10 ns recording interval. Illumination intensity was set such that %2.5 × In vivo FLIM-FRET measurement U2OS cell line was maintained in Dulbecco's modified Eagle's medium (DMEM, Gibco), supplemented with 10% Fetal Bovine Serum (FBS, Gibco), and 50 IU/mL penicillin and 50 mg/ mL streptomycin (Gibco) at 37˚C in a humidified atmosphere with 5% CO2. Cells were seeded on a 25-mm diameter, #1.5-thickness, round coverglass coated with poly-D-lysine (GG-25-1.5-pdl, neuVitro). Transient transfection of pCMV-mTurquoise2-GFP plasmid was done with TransIT-2020 (Mirus), and cells were imaged 24 hours later. During imaging, the cells were maintained at 37˚C on a custom-built temperature controlled microscope chamber, while covered with 1.5 ml of imaging media and 2 ml of white mineral oil (VWR) to prevent evaporation. The excitation wavelength was 850 nm, and the emission filter was 470/40 (Chroma). The excitation laser power was adjusted to 4 mW. Becker and Hickl SPCM acquisition parameters were set to 10x zoom, 256 × 256 image pixels, 5 second integration, and 256 ADC resolution.

Software Implementation
All algorithms were implemented in MATLAB. The code used is freely available on Github at https://github.com/bryankaye1/bayesian-analysis-of-fluorescent-lifetime-data. Posterior  distributions were generated by evaluating the likelihood function in a grid space of parameter values and were marginalized before estimation of the mode and mean for each parameter.

Results
In a sample where only a subset of fluorophores are undergoing FRET, photon emission distributions take the form of a biexponential distribution, with some fraction of the distribution consisting of photons from a short-lifetime exponential, another fraction consisting of photons from a long-lifetime exponential, and some fraction coming from a spurrious background distribution. The goal of FLIM analysis is to infer the relative weights of these distributions, along  with the lifetimes of the two exponential distributions, from the measured histogram of photon arrival times (Fig 1A). Here we apply an analysis based on Bayesian inference in order to infer the most likely set of parameters from experimentally measured data. The output of our algorithm is a posterior distribution, which gives the relative probability of measuring a given set of parameters (Fig 1B). To characterize our approach, we test our analysis in both the low-photon and low-fraction regimes, representing two extremes where data may be collected.

Low-Photon Regime
While the biexponential nature of FLIM histograms is apparent when the histogram is constructed using a large number of photons (Fig 1A), the histogram's underlying distribution is less obvious when the photon count is low (Fig 2A). Previous work has estimated the minimum number of photons necessary to achieve a certain accuracy in determining fluorescence lifetimes from TCSPC measurements [16]. In this regime it can be difficult to extract accurate estimates of the fraction of short-lifetime photons through methods that rely on histogram fitting. This low-photon regime is relevant in many applications of FLIM, due to the fundamental tradeoff between the number of photons collected and both the spatial-temporal precision of the measurement and the light dose received by the sample. Thus, methods that can improve the precision and accuracy of parameter estimation in the low-photon count regime could potentially lead to a practical increase in spatial-temporal resolution and lower light doses.
In order to test the accuracy and sensitivity of our analysis, fluorescence lifetime measurements were taken using Erythrosin B and Coumarin 153, two reference dyes with well characterized lifetimes of 0.47 ± 0.02 ns and 4.3 ± 0.2 ns respectively [17]. These dyes were mixed at a fixed ratio, and fluorescence lifetime measurements were taken (Fig 2A, Materials and Methods) in order to generate a master list of photon arrival times. A fixed number of photons were randomly sampled from the master list in order to construct a histogram of photon arrival times, and analyzed to infer an estimate of the fraction of short-lifetime fluorophores, f S , taken as either the mean or the mode of the posterior distribution, while the known lifetimes were held fixed (Materials and Methods). This process was repeated 300 times in order to produce an error estimate for each given photon count, and was repeated for total photon counts spanning % 3 orders of magnitude (Fig 2B).
We find good agreement between the estimates of the fraction of short-lifetime photons for total photons counts larger than % 200 photons, using either the posterior mean or posterior mode as a fraction estimate (Fig 2B). Slight discrepancies between estimates using the posterior These histograms were then analyzed in order to estimate the fraction of short-lifetime photons, f S . (B) The estimated short-lifetime fraction, f S , varies linearly with the constructed short-lifetime fraction for three different total photon numbers, with a small offset. Squares: estimate from posterior mode. Dots: estimate from posterior mean. Dashed lines: linear fits with slopes 0.9933 ± 0.0026, 1.0085 ± 0.0024, and 1.0106 ± 0.0031, and offsets of 0.002000 ± 0.0484 × 10 −4 , 0.004400 ± 0.0814 × 10 −4 , and 0.009500 ± 0.1992 × 10 −4 for low, medium, and high intensities respectively (95% confidence interval). Intensities correspond to data collected at %1.5 × 10 5 , 1.2 × 10 6 , and 4.8 × 10 6 counts per second, for low, medium, and high intensity respectively. Inset: mean and posterior mode are apparent due to truncation and the fact that the posterior distribution is skewed (Fig 1B), and thus in general the mode and the mean of the distribution are not equal. As a measure of the error in our parameter estimation, we compute the standard deviation of the estimates from the 300 numerical replicates (Fig 2C) for each photon count. Fitting a power law to all data points except for the four smallest photon counts yields an exponent of −0.48 ± 0.04 (95% confidence interval), consistent with the exponent of −0.5 predicted from the central limit theorem in the limit of high n photon . This ffiffiffiffiffiffiffiffiffiffiffi n photon p scaling is also evident for other analysis methods, including the rapid lifetime determination method [9], which has comparable error for high n photon .

Low-Fraction Regime
We next tested our results in the regime where a relatively large number of photons are collected, but the fraction of photons originating from the short-lifetime component is low. This regime is relevant in systems where a large number of donor molecules are present, but interactions leading to FRET are relatively rare. In order to test the performance of our algorithm in this regime, fluorescence lifetime measurements were taken of Erythrosin B and Coumarin 153 as representative short-and long-lifetime dyes respectively. Unlike the measurements taken in the low-photon regime, separate fluorescence lifetime measurements were taken for each dye, generating separate master photon histograms (Fig 3A). A fixed number of photons could then be numerically sampled from each master histogram in order to create test histograms containing a prescribed fraction of photons originating from the short-lifetime dye, which were then analyzed in order to estimate the short-lifetime fraction while the lifetimes were held fixed at their previously measured values (Materials and Methods). Data was collected at %1.5 × 10 5 , 1.2 × 10 6 , and 4.8 × 10 6 counts per second, corresponding to low, medium, and high intensity respectively, and histograms from each intensity were analyzed separately.
Photons were sampled from master curves such that the total number of photons was fixed at 5 × 10 7 , with a prescribed fraction of photons originating from the short-lifetime distribution. This process was repeated 100 times for each condition. Across orders of magnitude, the short-lifetime fraction estimated from our algorithm varies linearly with the prescribed shortlifetime fraction (Fig 3B), with linear fits giving slopes of 0.9933 ± 0.0026, 1.0085 ± 0.0024, and 1.0106 ± 0.0031, and offsets of 0.002000 ± 0.0484 × 10 −4 , 0.004400 ± 0.0814 × 10 −4 , and 0.009500 ± 0.1992 × 10 −4 for low, medium, and high intensities respectively (95% confidence interval). The estimated short-lifetime fraction differs from the known short-lifetime fraction by a small bias factor, evident by the small positive offsets in the linear fits (Fig 3B, Inset). We hypothesize that this offset may be due to a number of factors, including non-monoexponential photon emission from the dyes, slight mischaracterization of the lifetimes or the instrument response function, or an intensity dependence of the FLIM measurement system. While the magnitude of the bias varies with intensity, the magnitude of the bias is relatively small, overestimating the fraction by less than one percent for the highest intensity tested.
In many applications, the changes in FRET fraction are more relevant than the actual fraction values themselves. Thus, we next considered the accuracy of measuring changes in the short-lifetime fraction, which were derived from the results in Fig 3B by subtracting values adjacent on the short-lifetime fraction axis. While the estimated short-lifetime fractions contain a small bias (Fig 3B), the bias is largely removed when changes in short-lifetime fraction are considered (Fig 3C). Consistent with this removal of bias, fitting linear equations to the estimated changes in short lifetime fraction vs. prescribed short lifetime fraction gives slopes of 0.9573 ± 0.1895, 1.0013 ± 0.1445, and 0.9580 ± 0.1717, and offsets of (0.2332 ± 0.3712) × 10 −5 , (0.3215 ± 0.3088) × 10 −5 , and (0.4617 ± 0.4286) × 10 −5 for low, medium, and high intensities respectively (95% confidence interval) (Fig 3C). These results demonstrate the accuracy and precision of our method for measuring changes in short-lifetime fraction across many orders of magnitude. For a short-lifetime fraction of 2 −7 , the sample standard deviation decays with increasing photon number. Fitting a power law yields an exponent of −0.4764 ± 0.0471 (95% confidence interval), consistent with the exponent of −0.5 predicted from the central limit theorem and as was the case for the low-photon regime measurements (Fig 2C).

In vivo Testing and Method Comparison
As FLIM is commonly used to measure FRET in living systems between biological fluorophores, which may contain complications not accounted for in our model, we next tested the applicability of our method in living cells. FLIM measurements were carried out on U2OS cells transfected with a plasmid carrying mTurquoise2-4AA-GFP, a fusion between the FRET pair of mTurquoise2 and GFP ( Fig 4A). As these two fluorophores are physically attached to each other in close proximity, a fraction of the donor mTurquoise2 molecules undergo FRET, and thus have a short lifetime. However, as these fluorophores must undergo maturation before being functional, some fraction of mTurquoise2 molecules will be attached to GFP that are not fully mature and thus will not undergo FRET, leading to a long-lifetime fraction. Photons were pooled from pixels grouped using boxcar windowing into groups of either 3 × 3, 7 × 7, or 11 × 11 pixels and analyzed using either the Bayesian analysis presented here, or least-squares fitting. (B) Histograms showing the probability density of the long-lifetime fraction from images in (A). The probability density functions from Bayesian analysis were found to have mean values of 0.648 ± 0.069, 0.642 ± 0.037, and 0.640 ± 0.026 (mean ± s.d.) for 3 × 3, 7 × 7, and 11 × 11 binning respectively, while the means values were found to be 0.558 ± 0.137, 0.648 ± 0.050, and 0.652 ± 0.034 (mean ± s.d.) for 3 × 3, 7 × 7, and 11 × 11 binning respectively using least-squares-fitting. To test the performance of our method as a function of the photon number, photon arrival times were pooled within the cell using boxcar windowing, from areas of either 3 × 3, 7 × 7, or 11 × 11 pixels corresponding to average photon counts of 1,087 ± 260, 5,826 ± 1,451, and 14,105 ± 3,673, respectively (Fig 4B). In order to more readily make comparisons with the least-squares method, we here infer the relative amplitudes of the biexponential decay, instead of the relative photon populations previously considered, and thus consider the long-lifetime amplitude fraction instead of the long-lifetime photon fraction considered above. Using the Bayesian method, the distributions of the long-lifetime fraction were found to have mean and standard deviation values of 0.648 ± 0.069, 0.642 ± 0.037, and 0.640 ± 0.026 (mean ± s.d.) for 3 × 3, 7 × 7, and 11 × 11 binning respectively, showing little change for all conditions.
In order to compare the results from the Bayesian method presented here with the more commonly used least-squares fitting method, we repeated our analysis using fitting routines built into the Becker & Hickl software (Fig 4A and 4B). Mean and standard deviation values of the long-lifetime fraction distributions were found to be 0.558 ± 0.137, 0.648 ± 0.050, and 0.652 ± 0.034 (mean ± s.d.) for 3 × 3, 7 × 7, and 11 × 11 binning respectively. While the mean value for the long-lifetime fraction is similar for higher photon counts, there is significant discrepancy for 3 × 3 binning, where an asymmetric distribution is evident. Furthermore, for the highest photon counts, the mean values agree within error for the Bayesian method and leastsquares fitting, indicating a convergence between the two methods in the limit of high photon number. Thus, for low-photon counts, the Bayesian method presented here provides long-lifetime fraction estimates that are more accurate and precise than the nonlinear least-squares method.

Discussion
Here we presented an extension of previous Bayesian inference approaches to FLIM data analysis that takes into account additional experimental complexities. Using controlled experimental data as a test case, we show that this analysis performs remarkably well in both the lowphoton and low-fraction regimes.
In the low-photon regime, we can estimate the low-lifetime fraction, f S , with a precision of 0.003 and a bias of 0.017 using only 200 photons. At a photon collection rate of 2 × 10 5 photons per second, this number of photons corresponds to an acquisition time of only 1 millisecond. As the precision scales as / n À 1=2 photon (Fig 2C), if one instead wanted a higher precision of 0.001, one could instead collect data for 9 milliseconds. In the low-fraction regime, using 5 × 10 7 photons, for a short-lifetime fraction, f S , of 0.0156, we find a precision of 0.000096 and a bias of 0.0046. With an acquisition rate of 1.5 × 10 6 photons per second, this corresponds to % 33 seconds of acquisition time. As the precision in this regime also scales / n À 1=2 photon (Fig 3C), if one requires a higher precision of 0.000032, this could be obtained by acquiring data for nine times as long, or 300 seconds. Thus, in both the low-photon and low-fraction regimes, our results show the required number of photons, and hence the acquisition time, necessary to achieve a given level of precision.
One limitation of our implementation is that we evaluate the posterior distribution at equally spaced points. A large parameter space must be searched, and the analysis presented here is relatively slow compared to other parameter searching techniques. For example, when 4 parameters are searched using a Markov chain Monte Carlo approach to stochastically optimize our likelihood, the computation time is reduced by a factor of %10-20 with no loss of accuracy (S1 Fig). Here we have focused on the use of FLIM to measure changes in FRET, yet it has wider applications, including in metabolic imaging [18] and in measuring local changes in environment, including pH [19] as well as oxygen [20] and Zn 2+ [21] concentrations. The analysis presented here is general, and should be applicable to FLIM measurements in these other systems as well.