uFLIM -- Unsupervised analysis of FLIM-FRET microscopy data
Francesco Masia, Paola Borri, Walter Dewitte, Wolfgang Langbein
uuFLIM – Unsupervised analysis of FLIM-FRET microscopy data
Francesco Masia,
1, 2, ∗ Walter Dewitte, Paola Borri, and Wolfgang Langbein School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK School of Biosciences, Cardiff University, Museum Avenue, Cardiff CF10 3AX, UK
Despite their widespread use in cell biology, fluorescence lifetime imaging microscopy (FLIM)data-sets are challenging to analyse. This is because, in each spatial position, there is typically asuperposition of multiple fluorescent components. Here, we present a data analysis method makingthe most out of the information in the available photon budget, while being fast and unbiased.The method, called uFLIM, determines spatial distributions and temporal dynamics of multiplefluorescent components with no prior knowledge. It goes significantly beyond current approacheswhich either assume the functional dependence of the dynamics, e.g. an exponential decay, orrequire dynamics to be known. Its efficient non-negative matrix factorization algorithm allows forreal-time data processing. We show that uFLIM is capable to disentangle the spatial distributionand spectral properties of 5 fluorescing probes, from only two excitation and detection channels anda photon budget of 100 detected photons per pixel. By adapting the method to data exhibitingF¨orster resonant energy transfer (FRET), we retrieve the spatial and transfer rate distribution ofthe bound species, without constrains on donor and acceptor dynamics. The deployment of thispowerful unsupervised method allows to extract the maximum information from FLIM-FRET data.
I. INTRODUCTION
Fluorescence microscopy is a widely used tool to studythe distribution of biomolecules in living cells and tis-sues, with high contrast, specificity, and spatial resolu-tion. The decay dynamics of the fluorescence intensityfollowing pulsed excitation can reveal information on thelocal environment of the emitting fluorophore. This con-cept is used in fluorescence lifetime imaging microscopy(FLIM), where spatially-resolved emission dynamics arerecorded [1]. Typically, the emission intensity is mea-sured as function of the delay after an excitation pulse,but there are also frequency-domain implementations [2].Spatially-resolved fluorescence dynamics have been usedto sense local variation of temperature [3], pH [4, 5], andionic concentration [6]. FLIM can also be used to dis-tinguish multiple spectrally overlapping fluorophores viatheir different decay dynamics [7].Among the various processes which alter the lifetimeof an emitter, F¨orster resonant energy transfer (FRET)offers the possibility of studying protein-protein inter-action [8, 9]. Here, two proteins of interest are taggedwith different fluorophores, called donor and acceptor.The emission spectrum of the donor spectrally overlapswith the absorption spectrum of the acceptor, and thelight excitation is chosen resonant to the absorption ofthe donor. If the distance of the two fluorophores issmall enough, typically in the nanometre range, signif-icant non-radiative transfer of the excitation occurs fromthe donor to the acceptor. Such energy transfer can bedetected by a quenching of the donor emission and a cor-responding enhancement of the acceptor emission. Forhigh accuracy and sensitivity, a method to detect thetransfer not relying on absolute intensities is preferable, ∗ Electronic address: [email protected] and this can be achieved by measuring the change ofthe fluorescence dynamics in FLIM. The energy transferprovides an additional relaxation channel for the donor,increasing its decay rate, and a corresponding delayedexcitation of the acceptor. Notably, FLIM-FRET is notaffected by absolute intensity changes originating fromphotobleaching, illumination inhomogeneity and/or con-centration distributions.To analyze FLIM, a common approach is to fit the sig-nal decay assuming a mono- or bi-exponential decay be-haviour. Recently, global analysis methods offering fasteralgorithms compared to pixel by pixel fitting have beenreported [10, 11]. These methods assume exponential de-cay dynamics, and the instrument response function intime-domain needs to be known to extract the exponen-tial time constants. FRET is observed as and additionaldecay rate and can be extracted from the fit parame-ters [10, 11]. However, while being a convenient mathe-matical function to use, and the simplest solution of rateequation models, an exponential decay is not necessarilyrepresenting the physical behaviour of a fluorophore.Phasor analysis is an alternative analysis approach[12–14]. In this method, each FLIM pixel is representedby two quantities, namely the real and imaginary partof the Fourier coefficient of the first harmonic (typicallyreferring to the excitation repetition rate) normalized tothe amplitude of the zeroth harmonic. These values arethen interpreted as coordinates in the resulting ”phasorplot” in the complex plane. Pure exponential decay dy-namics of varying decay times are forming a semi-circlein this plot. Due to the linearity of the transform, mixedexponential decay dynamics are resulting in averages ofpure component phasors, and thus have amplitudes in-side the circle. The phasor analysis provides a usefulqualitative tool when applied to FLIM-FRET data. Theoccurrence of energy transfer can be identified as a devia-tion of the phasors from the values obtained in regions ofthe sample occupied only by unbound donor molecules. a r X i v : . [ phy s i c s . b i o - ph ] F e b More quantitative information on the FRET process canbe obtained if some assumptions are imposed, e.g. amono-exponential decay of the unbound donor fluores-cence.FLIM data can also be analyzed by linear unmixingof the intensity decay in the basis of selected referencepatterns [15] measured a priori in samples with similarproperties as the sample under investigation. In thismethod, each fluorescence decay is approximated as alinear combination of reference decay curves by minimiz-ing the Kullback-Leibler discrepancy (KLD) [16], whichmaximizes the likelihood of the model for data showingPoisson noise, i.e. the noise expected from the photonemission statistics. Typically, a gradient descent methodusing multiplicative update rules is used to find the non-negative fractional concentrations of the reference pat-terns [16]. The reference patterns are either extractedfrom singly labelled control samples or by selecting re-gions of the image which are assumed to show the indi-vidual components. This approach has been applied tothe analysis of multispectral time-domain FLIM, and upto nine different fluorophores could be visualised [7]. Thesupervised determination of the components and theirdynamics is generally complicating the analysis, and in-troduces a bias. Specifically in the analysis of FRET-FLIM data, where the different components interact, re-liably determining the individual component dynamics ischallenging.Recently, Smith and co-authors [17] have developed adeep learning network to analyse FLIM and FLIM-FRETdatasets. However, the benefits of the free-fit approachare superseded by the typical shortcoming of deep learn-ing, i.e. the long computational time required to gen-erate the training set and to train the neural network.Moreover, the training dataset is generated assuming abi-exponential functional dependence of the fluorescencedynamics and require the knowledge of the instrumentresponse.In this work, we demonstrate an unsupervised FLIManalysis (uFLIM) method, using a fast non-negative ma-trix factorization (NMF) algorithm [18] and random ini-tial guesses for both the spatial distribution and decaytraces of the factorization components. Similarly to thepattern-based unmixing, the NMF method decomposesthe data into a linear combination of few components,but in this case the method does not require the knowl-edge of the component patterns which can be deducedas part of the factorization. The method offers all of theadvantages of pattern unmixing, i.e. there is no specificassumption on the analytical dependence of the fluores-cence dynamics, and as a consequence there is no needto know the instrument response function. Notably, itprovides these benefits at higher speed and dropping therequirement of prior knowledge of reference patterns. Wedemonstrate the performance of uFLIM in distinguish-ing multiple spectrally overlapping fluorescing proteins,showing that the method can retrieve the spatial dis-tribution and dynamics of 5 fluorescent protein probes while using a simple FLIM set-up with only two excita-tion lasers and two detection channels.Building on this method, we introduce a FRET anal-ysis, which uses the donor and acceptor dynamics de-termined by uFLIM from samples or sample regions notshowing FRET. Energy transfer is quantitatively char-acterised by the quantum efficiencies of donor and ac-ceptor emission into the detection channels, as well asthe mean and variance of the transfer rate distribution,here assumed to be log-normal [19]. The analysis deter-mines the values of these quantities, and thus the FRETpair dynamics, together with the spatial distribution ofthe FRET pairs, by minimizing the NMF factorizationerror. Other components, such as autofluorescence, canbe retrieved at the same time without prior knowledge.Notably, uFLIM-FRET does not assume a functional de-pendence of the fluorescence dynamics, but calculates theaverage donor and acceptor dynamics in the FRET pairfrom their unbound dynamics using the distribution ofFRET rates.
II. RESULTSA. uFLIM analysis method
Firstly, measured FLIM data are reshaped as a ( N s × N t ) matrix D where N s and N t indicate the number ofspatial and temporal points, respectively. Then, we intro-duce N c components to represent the data, and use NMFto determine the spatial distribution matrix S of N s × N c elements, and the dynamics matrix T of N c × N t ele-ments. Multiple spectral channels, if present, are stackedin N t sequentially, and multiple data can be analysed to-gether by stacking them in N s . We assume in the follow-ing that D is given as number of detected photons, whichhas Poissonian noise with a standard devation given by (cid:112) D ij . We utilise a fast NMF algorithm which mini-mizes the residual || D − ST || , where || . || indicates theFrobenius norm [18]. This method provides the decompo-sition of maximum likelihood in case of Gaussian whitenoise in the data, i.e. a noise independent of the datavalue. To be able to use this algorithm, which is 2-3 or-ders of magnitude faster than gradient descent methodsaccounting for non-white noise, we partially whiten thedata before factorization (see Methods). The partiallywhitened data D w is then factorized by NMF, minimiz-ing E = || D w − S w T w || , and the resulting decomposition S w and T w is de-whitened to recover the factorization ofthe original data. This treatment of noise is providingequivalent results to minimizing the KLD for the dataconsidered (see SI). B. uFLIM application I: Single spectral channeldatasets
Here we demonstrate the uFLIM analysis of experi-mental data reported in Ref. [20, 21], in which the fluo-rescence lifetime of a dye changes due to variations in theenvironmental conditions, specifically the T2-AMPKARconstruct in presence of the 991 activator resulting inFRET. In these measurements, a single channel detectsthe dynamics of the T2-AMPKAR compound using time-correlated single photon counting (TCSPC) with 50 pstime bins. We have analysed the data in Fig. 4 ofRef. [20], with 4 × µ s/pixel for a singleuFLIM step on an Intel i7-8700 CPU. We found thatconvergence (error change below 1 ‰ per iteration) wasreached within 5 iterations. Further computational timesreported below were measured on the same CPU.The dynamics of the three components are given inthe bottom plot of Fig. 1. The first component T isinconsistent with fluorescence, and is attributed to sys-tematic artefacts in the data, for example due to non-linearity of the TCSPC at high count rates. The sec-ond component T is slowest, and is thus attributedto the emission of T2-AMPKAR without 991. Thedynamics of the third component T is faster, and isattributed to the T2-AMPKAR - 991 pair. The rel-ative integrated counts of the components, given by r i = (cid:16)(cid:80) s,t S s,i T i,t (cid:17) / (cid:16)(cid:80) i,s,t S s,i T i,t (cid:17) , are r = 0 . r = 0 . r = 0 . T and T .For visualization, the spatial distributions S and S are encoded using a hue-saturation-value (HSV) colourmapping at maximum saturation. The value (V), whichis the brightness, is taken as square root of S + S , nor-malised for each image. The hue (H) is given by thepoint-wise contrast ( S − S ) / ( S + S ), offset and scaledas indicated. We observe a change of colour of the HSVmaps from green to violet with increasing concentrationof 991, showing increasing fraction of T2-AMPKAR with991 attached. To compare with the global fitting expo-nential decay analysis in Ref. [20], the average lifetime (cid:104) τ (cid:105) = ( { S } τ + { S } τ ) / ( { S } + { S } ) is given in the in-set of the bottom panel in Fig. 1, where { . } indicates the1-norm, and τ and τ are the lifetimes of the individualcomponents gven by the first moments of their dynamics.The resulting (cid:104) τ (cid:105) exhibits a dependence on the activatorconcentration consistent with Ref. [20].This example shows that uFLIM is able to analyzeFLIM experiments with resulting lifetimes similar toglobal exponential fitting, yet providing the dynamics ofthe components not constrained to an exponential decay.Notably, measurement artefacts and autofluorescence canbe separated into additional components, and their influ- -2 0 2 4 6 80.000.020.040.06 T T T T delay (ns) á t ñ ( n s ) concentration ( m M) r m M0.1 m M1 m M10 m M25 m M50 m M -0.60.2 FIG. 1: Results of the uFLIM algorithm applied on a datasetfrom Ref. [20] of TCSPC FLIM on HepG2 expressing the T2-AMPKAR compound as a function of the concentration ofthe 991 activator. In these measurements, a single channeldetects the dynamics of the T2-AMPKAR compound. Thedata have been factorised into three components which showdifferent dynamics. Top: Concentrations S and S and dis-played with a HSV mapping as discussed in the text for dif-ferent concentrations of 991 as indicated. The contrast rangeis encoded as hue as given by the colour scale. Bottom: Tem-poral dynamics T of the three retrieved components. Inset:Weighted average lifetime (cid:104) τ (cid:105) (black symbols) and fraction of S in each image r = { S } / ( { S } + { S } ) (red symbols) forthe different fields of view versus activator concentration. ence on the analysis is accordingly reduced. As additionalexample, we shown in the SI Sec. S2 the uFLIM analysisof time-gated FLIM images of mixtures of two differentdyes, and its ability to recover their dynamics and distri-bution. C. uFLIM application II: Multiple spectralchannels and unmixing of many fluorescent proteins
Imaging living cells which are expressing multiple flu-orescent proteins (FPs) is crucial when disentangling theprotein interaction network. Here, we explore the ca-pability of uFLIM to extract the spatial distributionof a large number of FPs, by unmixing their spec-tral and temporal profiles. A similar question wasasked in Ref. [7] using the pattern-matching algorithmon spectrally-resolved fluorescence lifetime imaging mi-croscopy (sFLIM) data. sFLIM was employed with se-quential excitation at three wavelengths and detectionover 32 spectral channels. Up to nine fluorescent probescould be separated, for data having a photon budgetof around 1000 photons per pixel in the bright regions.However, this result required prior knowledge of fluores-cence decay and spectral signature patterns, a constrainthat can be lifted with uFLIM.To test the performance of uFLIM on sFLIM datasets,we generated synthetic data combining several FPs.Since the large number of excitation and detection chan-nels used in Ref. [7] are not available in most FLIMsystems, we simulate here a much simpler system withonly two excitation lasers (at wavelenghs of 460 nm and490 nm) and two detection channels (over wavelengthranges of 500–550 nm and 550–700 nm). We use anexcitation repetition rate of r = 40 MHz, a detectionrange from -1 ns to 24 ns with l = 1000 temporal chan-nels, and a gaussian instrument response function (IRF)exp( − t /w ) with w = 141 . f (see SI Table S1), with spatial distributions givenby selected paintings [22], which were cropped and re-sized to 256 ×
256 pixels, converted to greyscale usinga gamma of 1.5, and normalized to have unity mean,yielding the distribution matrix F f . For each combi-nation of excitation wavelength (index e ) and detectionchannel (index d ), we define a scaled spatial distribution F def = c def F f , where c def accounts for the quantum ef-ficiency and the extinction coefficient of the FP, and thefraction of photon emission by FP f detected by channel d , see SI. We also define the fraction of photons detectedin a given channel as ˆ c def = c def / (cid:80) d,e c def , and the frac-tion of detected photons contributed by a given FP asˆ c f = (cid:80) d,e c def / (cid:80) d,e,f c def . The measured FP dynam-ics over the time t , represented by the matrix T f , arecalculated as the convolution between the Gaussian IRFand a mono-exponential decay with decay rate γ f , see SI.The noiseless sFLIM synthetic data are then obtained bymultiplying the paintings with the dynamics, and adding their emission, assuming equal spatially-integrated num-bers of each FP, yielding D s ed = A (cid:88) f F def T f , (1)where the normalization A is ensuring { D s } = N s I tot ,and the average number of photons I tot per spatial pointwas chosen to be 100 or 10 in the results shown. To sim-ulate photon counting detection and corresponding noise,the integer values of a random variable following Poissonstatistics with a mean value given by the noiseless sFLIMdata are taken as photon counting sFLIM data. Compu-tational time was reduced by partially binning the 1000time channels in T f according to the method describedin the SI Sec. S3.This synthetic data is then analysed by uFLIM ac-cording to the method described in Sec. II A. As a firsttest, for direct comparison with Ref. [7], we retrieved thespatial distribution S w in a single step NMF, with thedynamics T w fixed by T f , i.e. assuming prior knowledgeon the dynamics. The resulting S are shown in Fig. 2for I tot = 10 . The spatial distributions of all 8 FPs arewell retrieved, with a root-mean-square error of about300 and about 20% relative error. We emphasise thatthis was achieved using a two-channels in excitation anddetection, compared to 3 and channels in Ref. [7]. FPswith properties differing significantly from each other arewell recovered, while more error is visible for FPs withsimilar properties, for example for mEos2 and mVenus,and for FPs with weak emission, such as LSS-mKate2.Even for a much smaller photon budget I tot = 100 (seeFig. S3), spatial distributions are recovered, albeit withaccordingly larger noise.To evaluate if the retrieval could be improved by max-imizing the likelihood for the Poisson statistics, we haveimplemented a gradient descent minimising the KLD. Wehave used a multiplicative update rule [16], and, as ini-tial guess of S , either the result of the NMF, or the so-lution of the linear system D = ST (see SI Sec. S4). Inboth cases, we did not observe a relevant improvementof the results compared to the fast NMF algorithm (seeFig. S5 and Fig. S4 for I tot = 10 and Fig. S7 and Fig. S6for I tot = 100), despite a computational time about 15–50 times longer. Using the fast NMF, the uFLIM com-putational time was 5 µ s/pixel. This indicates that thewhitening transformation applied by us, combined withfast NMF algorithm, is a suitable alternative to the com-putationally expensive gradient descent method.Next, we applied uFLIM to retrieve the spatial dis-tribution and the FP spectral and dynamic propertiesdirectly from the measured data, with no prior knowl-edge. Since determining the FP properties additionallyto the spatial distribution is more taxing on the data in-formation content, we have removed the three FPs withthe smallest differences in their properties from the 8previously used (see Table I). To simulate the presenceof unknown variations from the nominal properties, thesFLIM data are generated by introducing a ±
20% rela-
WasCFP BrUSLEE mBeRFP Dendra2(Red)MiCy mEos2(Green) mVenus LSS-mKate2
M=6.5e3, a=2.1e3
M=6.7e3, a=1.3e3 M=8.1e3, a=1.1e3
M=9.6e3, a=8.5e2
M=1.5e3, a=6.2e2 M=1.1e4, a=2.0e3 M=5.2e3, a=1.6e3 M=2.0e3, a=4.3e2
WasCFP
BrUSLEE mBeRFP
Dendra2(Red)MiCy mEos2(Green) mVenus LSS-mKate2
M=1.6e3, m=-1.7e3r=3.2e2 M=4.6e2, m=-4.9e2r=1.0e2 M=1.1e3, m= -1.4e3r=2.1e2 M=9.1e2, m=-6.9e2r=1.1e2M=5.2e2, m=-6.0e2r=1.1e2 M=4.0e3, m=-3.5e3r=6.6e2 M=2.7e3, m=-3.0e3r=5.1e2 M=1.2e3, m=-9.5e2r=2.0e2
FIG. 2: Spatial distributions obtained by applying uFLIM tosFLIM synthetic data generated with 8 fluorescent proteins(see labels), having spatial patterns given by selected paint-ings, and detected by a two-channel FLIM set-up (see text).uFLIM in this cases assumes prior knowledge on the FP dy-namics, i.e. uses fixed T f . Greyscale is from m to M . Toprows: Retrieved S with m = 0 and the maximum ( M ) as indi-cated. The spatially averaged pixel values ( a ) are also given,for comparison with the nominal values ˆ c f I tot , see Table S1.Bottom rows: Difference between the retrieved and nominaldistributions corresponding to the top rows, with M and m as given. r indicates the root-mean-square. tive variation in c def and 1 /γ f , and an additional ± c f , taken at random from a uniform dis-tribution. We use the iterative uFLIM method whereboth S and T are calculated. The nominal FP proper-ties, before parameter variation, are used to generate theinitial value of T , while the guesses for S are obtainedby solving the system D = ST and then setting negativevalues to zero. We constrain the dynamics of a given FPto be the same for all excitation and detection channels,by replacing at each NMF iteration step the dynamicscalculated for the different channels with their average.The iteration is stopped if the factorisation error has notimproved for three consecutive steps, allowing for a max-imum of 100 iterations. Here, a single iteration step tookabout 2 µ s/pixel, and typically 10-25 steps were used.Fig. 3 shows the retrieved spatial distributions and FPproperties obtained for I tot = 10 . The spectral and tem-poral properties extracted from the retrieved quantitiesare given in red in Table I, showing a good agreementbetween the retrieved and original S and T . Even for I tot = 10 (see Fig. S9), the retrieval works reasonably,showing only some cross-talk between the FPs with mostsimilar properties, mBeRFP and Dendra2(Red). Results Name of FP / f ˆ c f ˆ c f ˆ c f ˆ c f ˆ c f τ (ns) painting WasCFP / 1 0.27 0.07 0.53 0.13 0.36 5.05
The creation of
Adam
The Hay Wain
The ambassadors
Old woman and boy with candles
The great wave .
00 0.07 3.73 off Kanagawa I tot = 10 are given in red, where the lifetimes are the firstmoment of the retrieved dynamics for positive times. Thestandard deviations of the retrieved parameters are given ingreen. can be slightly improved by subsequently minimizing theKLD (see Fig. S14 and Fig. S15), however this takes 2 to5 times longer than the fast NMF depending on I tot andthe choice of initial guesses.We note that the number of FPs retrievable within acertain error depends in a complex way on their proper-ties, especially on their differences, as well as the signalstrength I tot , and the FP spatial distributions. There-fore, for a given experiment, a reliable determination ofthe retrieval error should be obtained via repeated re-trievals using new realizations of the photon counts D from probability distributions determined by the mea-sured counts. To give an example, for the parametersshown in Table I, by evaluating over 10 realisations ofthe photon shot noise, we found that the absolute devia-tions for ˆ c def and ˆ c f and the relative deviation for τ arebelow 1%.To exemplify the benefits of using retrieved propertiesversus fixed properties, we show in Fig. S17 the FP dis-tributions obtained from the data of Fig. 3 fixing the FPproperties to the nominal ones, not including the varia-tions introduced. Significant systematic errors are foundfor weak FPs, e.g. Dendra2(Red) and MiCy. With de-creasing I tot , the noise in the data is increasing and therelative importance of the systematic error decreases, sothat for I tot = 100 (see Fig. S18), these systematic errorsare less relevant.We stress that retrieving both the spatial distributionand the FP spectral and dynamic properties from the -5 -4 -3 -2 WasCFP BrUSLEE mBeRFP Dendra2(Red) MiCy
M=10347, a=4091 M=10927, a=2046 M=12735, a=1903 M=15077, a=1333 M=1270, a=649Mm M=1444, m=-147, r=593 M=238, m=-834, r=186 M=551, m=-1168, r=221 M=1038, m=-653, r=186 M=211, m=-930, r=409
FIG. 3: Spatial distributions and properties of 5 FPs retrieved by uFLIM from sFLIM synthetic data with I tot = 10 , generatedas described in the text, on a greyscale from m to M . uFLIM in this cases is applied with no prior knowledge on the FPproperties. Top row: Retrieved S with m = 0 and the maximum M as indicated. The spatially averaged pixel values ( a ) aregiven, having the nominal values ˆ c f I tot , see Table I. Middle row: Retrieved dynamics T f (red), and corresponding originaldynamics (black). Bottom row: Difference between the retrieved and original distributions, using M and m as given. r is theroot-mean-squares of the distribution differences. measured data directly, so that the separate measure-ments on reference samples with individual FPs are nolonger required, is a significant advantage of uFLIM. No-tably, the spectral and dynamic properties of FPs arevarying with their environment, and thus can be differ-ent between pure solutions and cellular samples. Fur-thermore, long-term drift of the instrument response canadd to systematic deviations between the FP propertiesused and the ones present in the sample of interest. Re-moving the need for such prior knowledge is therefore amajor advantage. D. uFLIM-FRET analysis method
Beyond the unsupervised analysis of FLIM data, wehave extended the algorithm to retrieve the spatial distri-bution of FRET pairs. In the literature, FRET efficien-cies are often derived from fitting the measured dynamicsby exponential decays and comparing the resulting decaytimes with the decay constant measured in samples whereonly the donor is present. These methods are limited bythe assumption of exponential decay dynamics, and re- quire the knowledge of the instrument response function.As demonstrated in the previous sections, in uFLIM,the temporal dynamics is retrieved without prior knowl-edge nor assumption of an exponential decay. ThereforeuFLIM can be applied to data showing pure donor andacceptor dynamics as components, providing the normal-ized pure donor and acceptor dynamics which we call T d and T a , with { T d } = { T a } = 1. FRET occurswhen donor and acceptor are in close proximity, forminga donor-acceptor pair (DAP). The FRET process intro-duces a non-radiative decay of the donor with rate γ (asketch of the energy diagram is shown in Fig. S19). Wethen describe the effect of FRET on the donor excitation,by subtracting the FRET transfer at points in the past,propagated to the present. More details on the calcu-lation of the modified donor and acceptor dynamics forDAPs are given in Sec. IV B.In the uFLIM-FRET method, a minimum of threeNMF components in T are used, given by T d , T a andthe calculated dynamics of a donor-acceptor pair ˜ T f . Thelatter is a function of the FRET rate γ and the ratio q ofthe detection probability of an excitation decay of accep-tor to donor. The spatial distributions S d , S a , and S f ,of these components are determined from the data byNMF. Additional components can be added to the NMFanalysis, for example to take into account autofluores-cence, determining their temporal dynamics and spatialdistribution without additional constrain starting fromrandom initial values, as shown for example in Sec. S15.The most likely values of γ and q given the data arethe ones minimizing the residuals of the NMF. The modelcan be expanded to several FRET components with dif-ferent rates, which is a typical situation in FRET due tothe variation in distance and relative orientation of thedonor and acceptor transition dipoles [23]. Such a vari-ation can be efficiently rationalized using a distributionof rates P ( γ ; ¯ γ, σ ) of mean value ¯ γ and relative standarddeviation σ , resulting in the FRET dynamics T f (¯ γ, σ, q ) = (cid:90) P ( γ ; ¯ γ, σ ) ˜ T f ( γ, q ) dγ . (2)Again, the most likely values of the parameters ¯ γ , σ ,and q are minimizing the residual of the NMF, which arefound using a computationally efficient method detailedin the SI Sec. S8. E. uFLIM-FRET application I: Analysis ofsynthetic data
To verify the method, we first use synthetic data. Weconsider two detectors which are mostly detecting thedonor and acceptor emission, respectively, given by R d = R a = 0 .
9. The FLIM system is the same as discussed inthe previous section. We consider that donor and accep-tor fluorescence have a single-exponential dynamics, withdecay rates of γ D =0.33/ns and γ A =0.385/ns, respec-tively, corresponding to the decay lifetimes of mNeon-Green and mCherry. We vary the spatially averagedtime-integrated photon counts of the donor emission I d ,proportional to the one of the acceptor emission, I a , us-ing I a = 0 . I d throughout. The dynamics of the DAPdetected by the two channels is calculated according toEq.(m9). In typical experimental conditions, we expectthat in the focal volume the donor-acceptor distance andorientation are not fixed to single values, but rather fol-low a distribution, leading to a distribution of FRETrates. Here, we have considered that the FRET ratesfollow a log-normal distribution [19] given by P ln ( γ ) = 1 γξ √ π exp (cid:18) − ( µ − ln γ ) ξ (cid:19) , (3)where µ and ξ are determined by ¯ γ and σ as ξ = (cid:112) ln ( σ + 1), and µ = ln(¯ γ ) − ξ /
2. Such a distributionis expected for a product of many statistically indepen-dent positive random variables. We generated data with¯ γ taking values of ¯ γ s = 0 . σ given by σ s = 0 .
5. The relative detection efficiencybetween donor and acceptor was taken to be q s = 1. Inthe following, symbols without subscript refer to the pa-rameter values changed by the algorithm, while symbols with the subscript s refer to the values used to generatethe data, and symbols with the subscript r are valuesresulting from the algorithm.Various relative strength of DAP and donor emission, I f /I d , are considered, where I f are the spatially aver-aged time-integrated photon counts of the DAP emission.As spatial distributions of donor, acceptor, and DAPwe used Monet’s Nymph´eas , Van Gogh’s
Starry Night ,and Leonardo’s
La Gioconda , respectively. The drawings[22] were cropped, resized to 256 ×
256 pixels, and con-verted to greyscale. The synthetic data D s = S s T s + b are then created by multiplying each pixel of the imageswith the corresponding decay curve, and adding the darkcounts b , which we characterize by their equivalent inten-sity I b = bN t . For the data shown, we have considered b = 0, 0.001, and 0.01, corresponding to I b = 0, 2, and20.The photon counting data D is generated from D s us-ing Poissonian statistics as before, and we repeated theanalysis for ten realizations of D . To reduce the analy-sis time, we apply some binning commencing 0.5 ns afterexcitation (see SI Sec. S3). The data are then factorisedusing the donor ( T d ), acceptor ( T a ), and FRET ( T f )components over a grid of the FRET parameters ¯ γ , σ ,and q . Donor and acceptor dynamics without FRET aretaken as known – in experiments these would have beenmeasured and retrieved by uFLIM. No free componentsare used, so that the factorization is a single step NMF forthe spatial distributions S w which minimise the residual E . The initial guesses for S w are random. The depen-dence of the factorisation error over the parameter spaceis shown in Fig. S25.Fig. 4 shows uFLIM-FRET results for the values of ¯ γ r , σ r , and q r minimizing the residual, for a specific data re-alisation with ¯ γ s = 0 . I b = 2 and κ = 1. Resultsfor small intensity and strong FRET ( I d = I f = 100)are given on the left, for large intensity and weak FRET( I d = 16 I f = 10000) in the middle, and for large inten-sity and strong FRET ( I d = I f = 10000) on the right.The first row shows the data summed over the tempo-ral channels, N t ¯ S , where the images of donor and ac-ceptor are visible, and the FRET image is discerniblefor strong FRET. The synthetic data dynamics T s , d , T s , a , and T s , f , are given as solid lines in Fig. 4 (bot-tom). The second, third and fourth rows from the topshow the spatial distributions S d , S a , S f retrieved byNMF, recovering the corresponding images well, alsoin conditions of small intensity (left) and weak FRET(middle). The difference between the original and re-trieved data is quantified using the relative error (cid:15) = || D s − ST || / || D s || , and similarly the reconstructionof the individual components is quantified by the rela-tive errors (cid:15) i = || S s , i T s , i − S i T i || / || S s , i T s , i || , where i ∈ { d , a , f } . The mean values ( (cid:104) . (cid:105) ) and the standard de-viations (] . [) of the reconstruction errors calculated overthe data realizations are shown in Fig. S27. As expected,the reconstruction error decreases with increasing inten-sities. We find that the error scales approximately as S f S a S d S t I d =100, I a =80, I f =100 M=619m=32M=321.79a=102.18M=339.80a=83.51M=408.82a=95.37 time (ns) F l uo r e s ce n ce dyn a m i c s I d =10 , I a =8 10 , I f =625 M=42088m=1724M=27403.46a=9983.42M=18845.74a=7945.66M=3969.42a=696.47 time (ns) I d =10 , I a =8 10 , I f =10 M=56770m=3852M=27297.56a=10000.49M=19626.75a=7999.06M=26261.68a=9999.55 time (ns)
FIG. 4: Results of uFLIM-FRET for synthetic data using ¯ γ s = 0 . σ s = 0 . q s = 1, I b = 2 and κ = 1. The three columnsrefer to different intensities as indicated. Top row: the time summed data N t ¯ S , on a linear grey scale from a minimum m (black) to a maximum M (white) as indicated. The second to fourth rows show the retrieved spatial distributions of the donor S d , acceptor S a , and FRET S f . Here m = 0 and a is the average pixel value over the image. The bottom panels show thesynthetic original dynamics of donor T d (black), acceptor T a (green), and DAPs undergoing FRET T f (blue), with the retrievedFRET dynamics given as red dashed lines. The signal of the donor (acceptor) detector are given in the top (bottom) panel,respectively. The dynamics are normalized to have a sum of unity over the 2000 temporal points of both detectors. Theretrieved FRET rate distribution parameters are ¯ γ r = (480 . ± . /µ s, σ r = 0 . ± . q r = 0 . ± . γ r = (502 . ± . /µ s, σ r = 0 . ± . q r = 0 . ± . γ r = (499 . ± . /µ s, σ r = 0 . ± . q r = 1 . ± . / √ I d (see Fig. S28). In general, (cid:15) , (cid:15) d , and (cid:15) a dependmostly on I d , while (cid:15) f is affected by both I d and I f . Wenote that all errors are below 10% for the high intensitycase, and that they are always much larger than the pa-rameter retrieval error, since they are dominated by theshot noise in the realizations.The uFLIM-FRET analysis is largely superior to thephasor analysis approach, as we show in the SI Sec. S10using the same data. Specifically, to extract quantitativeinformation, a phasor analysis needs to assume a sim-ple model of the dynamics, and the abundance of thedonor-acceptor pairs undergoing FRET and the FRETefficiency are hardly disentangled. The spatial distribu-tions of the donor-only and DAPs obtained with the pha-sor analysis poorly reflect the original distributions (seeSI Fig. S26).The retrieved dynamics of the FRET component T f (dashed lines in Fig. 4) agree well with the original syn-thetic dynamics, which is confirmed by the close matchof original and retrieved values of the parameters ¯ γ r , σ r ,and q r given in the caption. Their mean values andstandard deviations over the ensemble of realizations aregiven in Fig. S39. The errors decrease as the intensitiesincrease, showing that the method is correctly retrievingthe FRET parameters. The standard deviation, whichis due to the photon shot noise in each realization, israther similar for the different parameters, with σ beingretrieved with somewhat less accuracy as its influenceon the dynamics T f is less. However, we see also somesystematics for low intensities, in particular σ is underes-timated. To verify if this could be due to the remainingnon-whiteness of the noise in the analysed data D w wehave repeated the factorisation using the gradient de-scent minimizing the KLD, with the fast NMF results asinitial guesses. The retrieved spatial distributions andFRET parameters obtained with the two methods aregenerally very similar (see Fig. S72), confirming the suit-ability of the fast NMF algorithm on whitened data forthe analysis of data showing Poisson noise. Additionally,the gradient descent comes with two orders of magnitudelonger computational time of about 50 µ s/pixel for givenFRET parameters, and of the order of 5000 evaluationsare used to find the parameters minimizing the error,which makes it unsuitable for real time analysis.The difference in the accuracy among the different pa-rameters can be understood by looking at the curvatureof the reconstruction error along the directions definedby the parameters. The curvature is much smaller alongthe σ direction, resulting in a lower accuracy in the de-termination of this parameter (see Sec. S14).Further results for different ¯ γ s and I b are given in theSI Sec. S13. In the case of a small FRET rate ¯ γ s = 0 . S d component, T f differs from T s , f , and the value of σ is underestimated.For a higher rate ¯ γ s = 0 . γ = 0 . γ = 0 . I b = 0), the method isable to retrieve the correct parameters of the FRET dis-tribution with smaller error than for I b = 2 (see Fig. S37and Fig. S38), and retrieval is possible down to lower in-tensities I d = 32 and mean FRET rates of 0.5 and 0.9/ns.Conversely for large dark rate ( I b = 20), higher I d and I f are required for retrieval, see SI Fig. S40 and Fig. S41.The dependence of the reconstruction and FRET pa-rameter retrieval errors on the image size is analysed inthe SI Fig. S30. The systematic errors of the mean FRETparameters are not significantly affected by the numberof pixels N s . The standard deviation instead is scaling as1 /N s , steeper than the 1 / √ N s dependence expected forthe shot noise. We note that each pixel comes with itsown concentration in S , so that the number of photonsper retrieved information is independent of N s , as longas the number of spatial points is much larger than thenumber of FRET parameters.We have repeated the analysis in the case of negligi-ble direct excitation of the acceptor molecules, choosing κ = 0 in Eq.(m7), and accordingly do not include a pureacceptor component with a dynamics T a in the NMF.The corresponding increase in contrast, and reductionin free parameters, results in smaller errors of both re-construction and retrieved parameters, as shown in theSI Fig. S48 to Fig. S64. Finally, we show the ability ofuFLIM-FRET to retrieve the FRET parameters and theDAP spatial distribution in the presence of an additionalcomponent, such as autofluorescence, in Sec. S15. F. uFLIM-FRET application II: Analysis ofexperimental data
To show that uFLIM-FRET works well also with ex-perimental data, we analyse FLIM-FRET experimentson Arabidopsis roots co-expressing two tagged inter-acting transcription factors, SHORT-ROOT (SHR) andSCARECROW (SCR), with levels of both proteins ele-vated in the endodermis by expression under control ofthe SCR promoter (pSCR) [24]. In this system the SCRfactor is tagged with YFP (donor), while the SHR proteinis tagged with the acceptor RFP. Details on the experi-mental settings can be found in [24]. Only one channel,centered at the donor emission, has been acquired in theFLIM measurements. The data was binned both spa-tially (2 ×
2) and temporally (to a bin size of 100 ps).We have used uFLIM on images of roots expressing only pSCR::SCR:YFP to retrieve the donor and autofluores-cence dynamics. Then we have applied uFLIM-FRET onthe images of roots co-expressing pSCR::SCR:YFP and pSCR::RFP:SHR , using a logarithmic binning starting at3 ns and the time zero set to the peak of the autofluores-cence component. We therefore have R d = 1 and R a = 00 -1.5 0 1.500.050.100.050.1 5 10 20 M=9441.9 M=739.2M=651.2 time (ns) F l uo r e s ce n ce dyn a m i c s FIG. 5: uFLIM-FRET analysis of pSCR expressed SCRand SHR in the Arabidopsis root endodermis. The im-ages show the distribution of the three components used inthe uFLIM-FRET analysis on a grayscale as in Fig. 3 with m = 0. Red: donor ( pSCR::SCR:YFP ), green: autofluores-cence, Blue: DAP. The corresponding dynamics T d (red), T f (blue), and T u (green) are shown in the graph. and the DAP dynamics simplifies to T f (¯ γ, σ ) = (cid:90) P ( γ ; ¯ γ, σ ) ˜ T d ( γ ) dγ. (4)Fig. 5 shows the results of uFLIM-FRET, yielding ¯ γ r =0 . / ns and σ r = 0 .
01. This corresponds to a quencheddonor decay time of τ f = τ d / (1 + ¯ γ r τ d ) = 1 . E = 1 − (¯ γ r τ d + 1) − = 0 .
66, wherea donor lifetime of τ d = 3 .
57 ns has been estimated fromthe first moment of T d .The analysis reveals an accumulation of DAP in theendodermis of the root, where both donor and acceptorare expressed, in line with reported single-pixel lifetimeanalysis of these data [24] (see also the distribution of theretrieved average lifetime in Sec. S16). As with the syn-thetic data, using the gradient descent minimizing theKLD instead of fast NMF does not lead to significantchanges in the parameters (¯ γ r = 0 . / ns and σ r ∼ µ s/pixel for a single evaluation, similar to what ob-served for the synthetic data, and ∼
700 iterations arerequired to find the parameters minimizing the error.
III. DISCUSSION
We have demonstrated a data analysis method, whichwe call uFLIM, to analyze FLIM data in an unsuper- vised way. It employs a fast non-negative factorizationalgorithm on partially whitened data to infer the emis-sion dynamics and the spatial distribution of emittingmolecules. The method offers several advantages com-pared to other approaches in the analysis of FLIM dataavailable in literature. Firstly, it does not make assump-tions on the shape of the dynamics, which is instead thestarting point of standard fitting techniques. Secondly,the algorithm does not require reference patterns, whichare the component dynamics, as input. Importantly, themethod is fast, capable of analysing data in real timeon desktop computers. It can unmix spectrally resolvedFLIM images where several spectrally overlapping fluo-rescing probes are present, extending the multiplexingcapabilities of FLIM.Based on uFLIM, we developed uFLIM-FRET whichextracts FRET rates and spatial distributions. Here, theindividual donor (and acceptor if detected) emission dy-namics, which can be determined by uFLIM, is used tocalculate the DAP dynamics for a distribution of FRETrates. uFLIM-FRET determines the values of the dis-tribution parameters which minimize the residual of theNMF, at the same time as determining the spatial distri-bution of donor, acceptor, and DAPs. The distributionparameters characterize the fluctuations in the separa-tion and orientation of the donor and acceptor in theDAP, going beyond the approximation of a single FRETrate. Additional known or unknown components canbe added to the retrieval. uFLIM-FRET can estimatethe FRET parameters even in presence of a unknownautofluorescence signal yielding half of the total signal.The method can be adapted to retrieve donor, acceptor,and DAP dynamics without separate donor and acceptordata. Generally, the more information is available, themore parameters can be retrieved. The precision of re-trieval depends on the corresponding effect on the data –the larger the difference between, for example, donor, ac-ceptor, and FRET dynamics over the detected channels,the higher the precision.The fast NMF algorithm used allows analysis in realtime for typical FLIM microscopes. This speed comeswith an approximate treatment of the noise in the data,possibly leading to systematic errors in the retrieval. Wehave demonstrated that these errors are typically not sig-nificant by comparing the results with a gradient descentalgorithm taking into account the exact noise model ofthe data, and thus providing a correct maximum like-lihood estimation, at the cost of orders of magnitudelonger computational time.Notably, the method also offers the possibility to com-press the data of FLIM experiments into the spatial dis-tributions of few components, which facilitates the usageof FLIM-FRET as a high-throughput tool for cell biology.In order to enable widespread adoption of uFLIM-FRET as method of choice to analyze FLIM data,the corresponding software is provided. Informa-tion on the data underpinning the results pre-sented here, including how to access them, can be1found in the Cardiff University data catalogue athttp://doi.org/10.17035/d.2020.0115661402.
Acknowledgments
This work was supported by the Cardiff UniversityData Innovation URI Seedcorn Fund. F.M. acknowledgesthe Ser Cymru II programme (Case ID 80762-CU-148)which is part-funded by Cardiff University and the Eu-ropean Regional Development Fund through the WelshGovernment. P.B. acknowledges the Royal Society forher Wolfson research merit award (Grant WM140077).The authors thank Christopher Dunsby (Imperial Col- lege London) and Ikram Blilou (King Abdullah Univer-sity of Science and Technology) for kindly providing theexperimental data used in the manuscript. Discussionswith Peter Watson and Camille Blakebrough-Fairbairnare gratefully acknowledged.
Contribution statement
F.M. and W.L. developed the method, with input fromP.B. and W.D.. F.M. implemented the method and an-alyzed the data, with input from W.L.. All authors con-tributed to the manuscript writing. [1] M. Y. Berezin and S. Achilefu, Chem. Rev. , 2641(2010).[2] M. Raspe, K. M. Kedziora, B. van den Broek, Q. Zhao,S. de Jong, J. Herz, M. Mastop, J. Goedhart, T. W. J.Gadella, I. T. Young, et al., Nat. Meth. , 501 (2016).[3] K. Okabe, N. Inada, C. Gota, Y. Harada, T. Funatsu,and S. Uchiyama, Nat. Comms. , 705 (2012).[4] A. Orte, J. M. Alvarez-Pez, and M. J. Ruedas-Rama,ACS Nano , 6387 (2013).[5] F.-J. Schmitt, B. Thaa, C. Junghans, M. Vitali, M. Veit,and T. Friedrich, Biochim. Biophys. Acta , 1581(2014).[6] H. C. G. Alexandra V. Agronskaia, L. Tertoolen, J.Biomed. Opt. , 9 (2004).[7] T. Nieh¨orster, A. L¨oschberger, I. Gregor, B. Kr¨amer, H.-J. Rahn, M. Patting, F. Koberling, J. Enderlein, andM. Sauer, Nat. Meth. , 257 (2016).[8] Y. Sun, R. N. Day, and A. Periasamy, Nat. Protoc. ,1324 (2011).[9] A. Margineanu, J. J. Chan, D. J. Kelly, S. C. Warren,D. Flatters, S. Kumar, M. Katan, C. W. Dunsby, andP. M. W. French, Sci. Rep. , 28186 (2016).[10] S. Laptenok, K. Mullen, J. Borst, I. van Stokkum,V. Apanasovich, and A. Visser, J. Stat. Softw. , 1(2007).[11] S. C. Warren, A. Margineanu, D. Alibhai, D. J.Kelly, C. Talbot, Y. Alexandrov, I. Munro, M. Katan,C. Dunsby, and P. M. W. French, PLOS One , 1 (2013).[12] A. H. A. Clayton, Q. S. Hanley, and P. J. Verveer, J.Microsc. , 1 (2004).[13] M. A. Digman, V. R. Caiolfa, M. Zamai, and E. Gratton,Biophys. J. , L14 (2008).[14] C. Stringari, A. Cinquin, O. Cinquin, M. A. Digman,P. J. Donovan, and E. Gratton, Proc. Natl. Acad. Sci. U.S. A. , 13582 (2011).[15] I. Gregor and M. Patting, Pattern-Based Linear Unmix-ing for Efficient and Reliable Analysis of MulticomponentTCSPC Data (Springer International Publishing, Cham,2015), pp. 241–263, ISBN 978-3-319-15636-1.[16] D. D. Lee and H. S. Seung, in
Advances in Neural In-formation Processing Systems 13 , edited by T. K. Leen,T. G. Dietterich, and V. Tresp (MIT Press, 2001), pp.556–562.[17] J. T. Smith, R. Yao, N. Sinsuebphon, A. Rudkouskaya, N. Un, J. Mazurkiewicz, M. Barroso, P. Yan, andX. Intes, Proceedings of the National Academy of Sci-ences , 24019 (2019), ISSN 0027-8424.[18] J. Kim and H. Park, SIAM Journal on Scientific Com-puting , 3261 (2011).[19] N. Balakrishnan, S. Kotz, and N. L. Johnson, ContinuousUnivariate Distributions: Volume 1 (Wiley-Blackwell;,1994).[20] G. Chennell, J. R. Willows, C. S. Warren, D. Carling,M. P. French, C. Dunsby, and A. Sardini, Sensors ,E1312 (2016).[21] Data can be found from https://omero.bioinformatics.ic.ac.uk/webclient/?show=project-4553 .[22] The pictures used in this manuscript have beendownloaded from Wikipedia, https://commons.wikimedia.org/wiki/File:Michelangelo_-_Creation_of_Adam_(cropped).jpg , https://commons.wikimedia.org/wiki/File:John_Constable_-_The_Hay_Wain_(1821).jpg , https://commons.wikimedia.org/wiki/File:Hans_Holbein_the_Younger_-_The_Ambassadors_-_Google_Art_Project.jpg , https://commons.wikimedia.org/wiki/File:Peter_Paul_Rubens_-_Old_Woman_and_Boy_with_Candles.jpg , https://commons.wikimedia.org/wiki/File:Tsunami_by_hokusai_19th_century.jpg , https://commons.wikimedia.org/wiki/File:Meisje_met_de_parel.jpg , https://commons.wikimedia.org/wiki/File:Venus_botticelli_detail.jpg , https://commons.wikimedia.org/wiki/File:Georges_Seurat_-_A_Sunday_on_La_Grande_Jatte_--_1884_-_Google_Art_Project.jpg , https://commons.wikimedia.org/wiki/File:Monet_Water_Lilies_1916.jpg , https://commons.wikimedia.org/wiki/File:Mona_Lisa,_by_Leonardo_da_Vinci,_from_C2RMF_retouched.jpg , https://commons.wikimedia.org/wiki/File:Van_Gogh_-_Starry_Night_-_Google_Art_Project.jpg , https://commons.wikimedia.org/wiki/File:Agnolo_bronzino-_Ritratto_del_Nano_Morgante_come_bacco_(front)_,_1552,_galleria_Palatina.jpg .[23] I. V. Gopich and A. Szabo, Proc. Natl. Acad. Sci. U. S.A. , 7747 (2012).[24] Y. Long, Y. Stahl, S. Weidtkamp-Peters, M. Postma,W. Zhou, J. Goedhart, M.-I. S´anchez-P´erez, T. W. J. Gadella, R. Simon, B. Scheres, et al., Nature , 97(2017), ISSN 1476-4687.
IV. METHODSA. uFLIM
As stated in the main text, we assume that the FLIMdata D is given as number of detected photons, andhas Poissonian noise with a standard devation given by (cid:112) D ij , with the first index corresponding to the spa-tial and the second to the temporal position. To adaptthis data for use with the fast NMF method, which pro-vides the decomposition of maximum likelihood in caseof Gaussian white noise in the data, we partially whitenit before factorization as follows. We generate the time-averaged image ¯ S and the spatially averaged dynamics¯ T by averaging D along the temporal and spatial points,respectively,¯ S i = 1 N t (cid:88) j D ij , ¯ T j = 1 N s (cid:88) i D ij , (m1)The background-subtracted, partially whitened data D w are then defined as D w ij = D ij − b (cid:112) ¯ S i ¯ T j , (m2)with the average dark counts b , which can be measuredindependently. We assume here that b is independent ofposition, time, and spectral channel, as it is typically thecase for scanning time-correlated single photon count-ing, and note that inhomogeneous dark counts can besubtracted in the same fashion. In D w the data hasbeen divided by the expected standard deviation of thedata when factorized into the average spatial and tem-poral dependence. This method whitens spatially de-pendent time-integrated intensities, as well as spatiallyintegrated time-dependent intensities. D w is then fac-torized by NMF, minimizing E = || D w − S w T w || , andthe resulting decomposition S w and T w is de-whitenedto recover the factorization of the original data S ij = S w ij (cid:112) ¯ S i , T ij = T w ij (cid:113) ¯ T j , (m3)so that D ≈ ST + b . B. uFLIM-FRET
The modified donor dynamics ˜ T d in the DAP is calcu-lated using ˜ T d i ( γ ) = T d i − i − (cid:88) j =1 f j ( γ ) ˆ T d i − j +1 (m4) iterating along the temporal channel i = 1 , , .., l . Thisexpression contains the normalized and zero-centereddonor dynamicsˆ T d k = (cid:40) T d m + k − /T d m for k ≤ l − m + 1ˆ T d l − m T d l /T d2 l − m − k for k > l − m + 1 (m5)where m is the temporal channel at which T d is maxi-mum. In Eq.(m5), we have extrapolated the donor exci-tation decay beyond the last measured point l using thedecay observed over the extrapolation time interval priorto l . The FRET transfer f j ( γ ) at time j is given by themodified occupation of the donor excited state at thattime, the time-step ∆, and the transfer rate γ , f j ( γ ) = ˜ T d j γ ∆ . (m6)These equations determine the effect of FRET on thedonor excitation, by subtracting the FRET transfer atpoints in the past, propagated to the present using thedonor excitation response function ˆ T d . A sketch visualiz-ing Eq.(m5) is shown in Sec. S7. We approximate ˆ T d bythe measured donor emission dynamics, normalized to itsmaximum and starting from its maximum as time zeroof the response. This is adequate for FRET rates smallerthan the inverse time resolution of the measurements,and is consistent with the finite resolution of the datafor which it is used. To support this statement we havecompared the resulting dynamics with the analytical so-lution of the donor excitation modified by a single FRETprocess in the simple condition of a mono-exponential de-cay for the pure donor and a Gaussian IRF, as shown inthe supplementary information (SI) Sec. S6. The time-resolution limitation can be controlled by refining thesystem dynamics, for example by deconvolution of a re-sponse function before analysis. Note however that de-convolution is modifying the noise of the data from thesimple Poisson distribution of photon counts.The modified acceptor excitation dynamics is calcu-lated using the same approach as˜ T a i ( γ ) = κT a i + i − (cid:88) j =1 f j ( γ ) ˆ T a i − j +1 . (m7)This expression contains the normalized and zero-centered acceptor dynamicsˆ T a k = (cid:40) T a m + k − /T a m for k ≤ l − m + 1ˆ T a l − m T a l /T a2 l − m − k for k > l − m + 1 . (m8)The normalisation of T d and T a ensures the conserva-tion of the number of excitations by the transfer fromdonor to acceptor. The parameter κ quantifies the di-rect excitation of the acceptor molecules by the laser (seeFig. S19). While κ can be included in the parameters tobe retrieved by the method, we assume in the followingthat κ is known a priori , noting that it is given by therelative absorption cross-section of acceptor and donor3at the excitation wavelength, and can be determined in-dependently. We assume to have two spectral channels,and that the donor and acceptor emission is detecteddominantly by the respective channels, with some cross-talk. We therefore introduce the fraction of donor R d (acceptor R a ) emission detected by the donor (acceptor)channel, respectively. The dynamics ˜ T D ( ˜ T A ) detectedin the donor (acceptor) channel for a donor-acceptor pairundergoing FRET with rate γ is then given by˜ T D ( γ, q ) = R d ˜ T d ( γ ) + q (1 − R a ) ˜ T a ( γ )˜ T A ( γ, q ) = (cid:0) − R d (cid:1) ˜ T d ( γ ) + qR a ˜ T a ( γ ) , (m9)respectively. Here, we have introduced the ratio q ofacceptor to donor detection probability of an excitationdecay over both channels, to take into account the dif-ferent quantum efficiency of the acceptor and donor, andthe different probability of detecting an emitted photonin the two channels, including detector efficiency and fil-ter performance. The values of R d and R a can be simplymeasured using samples of only donor or acceptor fromthe ratio of the number of detected photons in donor versus acceptor channel. Determining q instead requiresto additionally determine the relative excitation rates ofthe donor versus acceptor molecules in the two samples,requires knowledge of relative molar concentration andrelative absorption κ . We have considered here the situ-ation where R d and R a have been measured, while q isretrieved as part of the retrieval process.Typically, the acceptor has a small absorption at theexcitation wavelength, so κ (cid:28)
1, and is hardly detectedby the donor channel, so 1 − R a (cid:28)
1. In the mainmanuscript we show the more challenging condition of κ = 1, where three components (donor, acceptor, andDAP) need to be included, while the simpler case κ = 0,showing an improved retrieval for a given photon budget,is given in the SI.In the NMF, the temporal points of donor and ac-ceptor channel are concatenated into the temporal di-mension of D . Three NMF components in T are used,given by the donor, (cid:2) R d T d , (cid:0) − R d (cid:1) T d (cid:3) , the acceptor[(1 − R a ) T a , R a T a ] and the calculated dynamics of adonor-acceptor pair ˜ T f ( γ, q ) = (cid:104) ˜ T D ( γ, q ) , ˜ T A ( γ, q ) (cid:105)(cid:105)