Resolution of 4 components in the same pixel in FLIM images using the phasor approach

In several cellular systems, the phasor FLIM approach has shown the existence of more than 2 components in the same pixel, a typical example being free and bound NADH. In order to properly quantify the concentrations and the spatial distributions of fluorescence components associated with different molecular species we developed a general method to resolve 3 and 4 components in the same pixel using the phasor approach. The method is based on the law of linear combination of components valid after transformation of the decay curves to phasors for each pixel in the image. In principle, the linear combination rule is valid for an arbitrary number of components. For 3 components we use only the phasor position for the first harmonic, which has a small error, while for 4 components we need the phasor location at higher harmonics that have intrinsically more noise. As a result of the noise in the higher harmonics, caused by limited photon statistics, we are able to use linear algebra to resolve 4 components given the position of the phasors of 4 independent components in mixtures of dyes and 3 components for dyes in cellular systems.


Introduction
Fluorescence lifetime imaging microscopy (FLIM) is a popular imaging technique increasingly applied in diverse areas of biological research [1][2][3][4][5]. However, quantitative multicomponent analysis of fluorescence lifetime data in the same pixel is not widely used. Yet this analysis could enable one to separate different components in the same pixel and quantify their amounts as chemical species change with time. In this case, the time range available will depend on the time needed to acquire a FLIM image with sufficient precision to perform analysis of multiple components in the same pixel. This concept is different from the common method based on the phasor approach to identify components that are spatially separated in the FLIM image. Here we focus on a general method to quantify 3 and 4 components in the same pixel. In other words, we aim at producing an image of the different components rather than a spatially resolved image of different lifetime values.
In many research labs, FLIM imaging is done using the popular time-domain instrumentation which is mostly limited to analysis of two components-either by exponential fitting of fluorescence decay or by phasor analysis [3][4][5][6][7][8][9][10]. Resolving multiple components in the same pixel from data collected in the time correlated single photon counting (TCSPC) approach involves fitting the intensity decays with multiple exponentials and then associating the exponential coefficients to their molecular fractions. This fitting approach is difficult to perform when the photon statistics are poor and when the lifetime components are in a similar time range. However, in the phasor method which can be applied to time-domain as well as frequency-domain data, this particular difficulty of the time-domain approach can be, in part, avoided. In this work we consider the case in which about 50 to 100 photons are collected in a typical pixel and show that the phasor approach can provide a natural and robust way to resolve multi-components (3 and 4) in a Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. single pixel [11]. To demonstrate our point, we will make use of the linear combination property of phasors. This approach simplifies the problem when the decay of a molecular species is not a single exponential, which occurs frequently with fluorescent proteins or for autofluorescence in which molecular species can decay as a sum of multiple exponentials. In such cases, a molecular species decays with a multi-exponential decay; but still it has a unique phasor that can be used to identify that molecular species.
For this analysis we assume that we know the phasor positions of individual molecular species either because we can measure them in independent samples or because there are identifiable pixels in an image where we have a pure chemical species (though not necessarily a single exponential decay). When such is the case, the application of the rule of linear combination of phasors reduces the problem to finding the parameters of the linear combination that represent the fractional intensity of each component at any given pixel. This rule assumes that the molecular species are not interacting. If they interact, then the linear combination rule will not apply and molecular species with different phasor values should appear in the analysis and their positions can be identified in the phasor plot and then in the image [12].
In this article we first describe the rule of linear combination of phasors. Then we use simulations to determine how the amount of noise in a pixel affects the resolution of the phasor in the same pixel for 3 and 4 components. Next we perform measurements of mixture of dyes of different chemical compositions and lifetimes to verify that the results of the simulation are reflected in the actual measurements. This comparison will show whether there are important sources of noise in the data that have been neglected in the simulations. Finally, we analyze biological samples which presumably could contain pixels with 3 components contributing to the overall fluorescence intensity in a pixel [8]. In our experiments we use dyes that paint the DNA of cells with the idea that in the same pixel we may observe the emission of several dyes. In this case we use 3 different dyes of different lifetimes and calculate, in each pixel, the relative contribution of the dyes to the overall fluorescence. We tried using 4 dyes to paint the nucleus of cells, but we have no evidence that we could have 4 different lifetimes in the same pixel. Instead, in the case of solutions we were sure to have 4 lifetime components in the same pixel.

Algorithms and mathematics
This following section describes the algorithms and equations used in this work. These equations have been published previously and are repeated here for the reader's convenience [11]. In the phasor approach, the fluorescence and harmonic signals collected from each point of the image are transformed to the phasor space using the following transformations [6,8].
where, I i,,j (t) is the intensity decay, g i,,j (ω) and s i,,j (ω) are the X and Y coordinates in the phasor plot, respectively; ω is the angular repetition associated with the frequency ( f ) of the laser (ω=2πf ), T is the period of the laser repetition frequency and n is the harmonic frequency; i, j are the X and Y coordinates of a pixel of the original microscopy image, respectively. In frequency-domain measurements, the transformation to the phasor plot uses the following relations, where, g i,j (ω) and s i,j (ω) are again the X and Y coordinates of the phasor plot, respectively, and m i,,j and f i,,j are the modulation and phase at the image pixel i, j. The longer lifetime is represented by the larger phase angle in the phasor plot. The distribution of phasor points originating from FLIM measurements appears on (for mono-exponential decays, points 1, 2, 3 and 4 in figure 1) or inside (for multi-exponential decays) the universal circle (figure 1). The linear combinations are shown by the blue lines ( figure 1). If each component has a multi-exponential decay, its location will be represented by phasors not in the universal circle, but the principle of linear combination remains valid. For decays described by sum of exponentials (2 in the derivation below), the intensity decay can be expressed as, The calculated phasor coordinates from this decay are given by, These expressions can be simplified using the definition of fractional intensity: The phasor coordinates of the mixture are given by, For multiple exponential decays the law of phasor addition holds and thus, The fractions are normalized using the following equation: Equations (11)-(13) are the system of equations that are used in this paper to solve for the unknown fractions f n .

Brief introduction to phasors
The final results (equations (11)-(13)) give the law of phasor addition, which was originally written by G.
Weber [13]. Although Weber did not use the term phasor nor did he propose a graphical-geometrical interpretation, his work inspired the modern phasor approach. The modern phasor plot was first described in the review by Jameson et al [14] and the vector addition property of phasors was first used to correct phase and modulation (frequency domain) data for background signals [15]. The field was dormant for some time, although Brand's group implemented the phasor approach to following calcium binding [16]; they did not, however, use the term 'phasor' to describe their method. The new millennium witnessed a resurgence of the phasor approach. In 2004, Clayton, Hanley and Vermeer [17] discussed the application of the graphical phasor approach (not termed as such) using a single light modulation frequency in FLIM systems. In this work the plots shown had the axes labeled as A and B. In 2005 Hanley and Clayton followed up this work with a paper formally naming their graphs AB plots [18] They applied this approach to cell data, using FLIM, to study the phosphorylation of the epidermal growth factor. In an important paper in 2005 [19], Redford and Clegg carefully reviewed the mathematics underlying the phasor approach, which they termed 'polar plots', and pointed out the advantages of such plots in FLIM (although they did not work with cells but rather carried out work in vitro in a fast mixing chamber). In 2007, Colyer et al [20], described a mathematical model and physical implementation for a low cost digital frequency domain FLIM (DFD-FLIM) system. They also presented an error analysis that showed the precise parameters for maximizing the quality of lifetime acquisition, and presented phasor plots on some systems. In 2008, Digman et al [6] applied the phasor approach to FLIM studies on live cells with particular attention to its use in dealing with autofluorescence as well as FRET. Since these early works, dozens of FLIM studies on cells have utilized the phasor approach (see for example [1,3,8,11,[21][22][23][24], to list just a few). In addition to FLIM work on cells, numerous papers have appeared using phasors to study various systems in cuvettes including studies on pyrene [25], fluorescent beads [24], protein systems [26][27][28], nucleic acids [29] and membranes [30,31].
2.3. The phasor approach to resolving multiple species in the same pixel In the case of 3 and 4 components, there exists a unique exact algebraic solution to the system of equations (11)-(13) which provides the three (or four) fractions ( f 1 , f 2 and f 3 ) for every point in the phasor plot under reasonable experimental error commonly obtained in FLIM [20] analysis. Hence, we focus here on the experimental resolution of 4 components in the same pixel. If the phasor point is displaced from the expected solution due to noise, we will still get the fractional intensity contribution of the components that give the specific point in the phasor plot. Notice that this observation simply means that we are propagating the error of the phasor location (due to statistics) to the components fractional intensities. This procedure can be done by error propagation rules in the linear system we are considering here, no matter the number of components, provided that these components are independent. In any case this approach is a simple algebra problem. The 4-component case will give a quadrilateral shape in the phasor plot (figure 1, identified by phasor points 1, 2, 3 and 4). In this case the unknowns are the 4 fractional components. Equations (11)-(13) and the system of equations cannot be resolved since we have 3 equations and 4 unknowns.
To obtain more than 4 parameters we make use of the phasor transformation into higher harmonics, with each harmonic producing a new pair of equations (the equations for G and S for the selected harmonics). We will use only the first harmonic for 3 component analysis. Using the second harmonic for 4 component analysis, the system of equation is determined for 4 fractional intensities and can be solve by linear algebra methods. The 4 equations for the linear combination of harmonics are written below: Where f 1 , f 2 , f 3 and f 4 are the unknown fractional intensities to be retrieved, g n h1 and s n h1 are the S and G components of the first harmonic for the 4 phasor positions of the pure components and g n h2 and s n h2 are the S and G components of the second harmonic, respectively. The terms in the right side of equations (14)- (17), G h1 , S h1 , G h2 and S h2 are the harmonics experimentally measured first and second harmonics at the pixel where we are attempting to find the 4 fractions f n . Note that this simple linear algebra system could not be solved if the harmonics were not orthogonal to one another. Since phasor coordinates of higher orders are orthogonal, the new pair of equations for each harmonics provides a solvable system, since it would be a linear combination and the determinant of the system would be different from zero. In other words, the system of equations (14)- (17), is determined using linear algebra methods to give the 4 fractional intensities f n . The algorithm is applied for each pixel of the image. Of course we assume that the harmonics of the pure components can be measured with a given precision. In the following we study, using simulations, the effect of noise in a pixel and its' effect on the precision of the fractional intensities at the pixel.
The algorithm we describe is super-fast. An image 256×256 can be resolved in 4 components in less than 1 s. Data acquisition and analysis was performed using the Globals for Images-SimFCS software (G-Soft Inc. Champaign, IL-USA) and the DIVER microscope described in [30]. The microscope is connected to a 320 Fast-FLIMbox (ISS, Inc., Champaign-Illinois, USA) to allow FLIM measurements and phasor analysis. The FLIMbox utilizes Spartan 6 FPGA and allows multifrequency data analysis up to 8-th harmonics with a fundamental frequency of 80 MHz. The description of the electronics used in the FLIMbox was previously presented (Colyer et al [20] 20). This electronics makes use of the photon counting frequencydomain method. The output consists of a series of harmonics of the decay. In this work we use the first and second harmonics, but higher harmonics are provided by the electronics up to about 7 harmonics. As shown in Colyer et al (Ref [20]), the error for each of the coordinate for each of the harmonics is proportional to the inverse square-root of the number of photons in each phasor harmonics.

Cell culture and staining with nuclear probes NIH3T3 cells were cultured in Dulbecco's Modified
Eagle Medium (DMEM, Thermo Fisher Scientific) supplemented with 10% Fetal Bovine Serum (FBS, Sigma-Aldrich), penicillin (100 U ml −1 ) and streptomycin (100 μg ml −1 ) at 37°C in a 5% CO 2 incubator. Twenty four hours prior to FLIM imaging, the cells were plated in fibronectin (Sigma-Aldrich) covered glass bottom dishes (Mattek). To improve dye staining in live cells, we permeabilized the cells with 0.1% TritonX-100 for 10 min at 4°C, washing twice before and after the treatment with 1X Dulbecco's Phosphate Buffer Solution (DPBS) [31]. Triple staining of the cells was performed sequentially. First, cells were incubated 5 min at 37°C with Rose Bengal (Sigma-Aldrich, dye content 80%, working solution 10 μgml −1 ), and washed twice with DPBS. After, cells were incubated at least 5 min with Acridine Orange (Sigma-Aldrich, 13 ng ml −1 ) and NucBlue™ Live ReadyProbes™ Reagent (Thermo Fisher Scientific, working solution was a 1/30th dilution of manufacturer´s recommended 1 drop/1.5 ml). All solutions were prepared in DMEM phenol-free, high glucose, HEPES medium (Thermo Fisher Scientific). We simulated images of 256×256 pixels with an average intensity of 50 counts at each point of the image ( figure 2(A)). We simulated a total of 21 images (256×256 pixels) with the fractional intensities, equal across the whole image (table 1). The values shown in the table are before normalization.

Simulations of a mixture of 4 components
For each point generated in the simulations we get values of G and S at two harmonics (80 and 160 MHz in our example). For the simulations we know the corresponding fractional intensities of the components that have generated the point. We then added noise to the G and S phasor coordinates of the point according to a Gaussian distribution with a given standard deviation, gradually increasing the amount of noise as shown in figure 2(B). As shown in Colyer et al [20], this is the standard deviation of the noise affecting each pixel of each harmonics evaluated. Hence, each of the pixels in the 256×256 image has random noise added to the phasor coordinates G and S. We then apply the linear algebra algorithm that recovers the fractional intensities. Our expectation is that increasing the noise on the phasor components will decrease the precision of the recovered fractional intensities as shown in figure 2(B). For each simulated point of the image we measure the mean deviation of the recovered value of the fractional intensities with respect to the input values of the fractional intensities used to calculate a given point. The simulation was done for a system composed of 3 and 4 components. The results are shown in figure 2(B), where the mean standard deviation of the recovered fractional  intensities is plotted against the noise added in units of G and S. The blue points are for 3-component simulation and the red points are for the 4-component simulation. For a relatively large range of errors the average error of the recovered fractional intensities is linear and, of course, the error in the fraction recovered is smaller for 3 components than for 4 components. An example of simulated data with noise is shown in figure 2(A) where the noise level was 0.1 (Gaussian with unity standard deviation) in the G and S coordinates for each harmonics for a 4 components simulation.

Simulation of linear combination of 4 components
The simulated data according the above linear combinations are shown in figure 3(A). For each image and each component, we show the difference between the input data and the recovered data in figure 3(A). All data are recovered relatively well. We note that this is a 4 component analysis and that the original image contains very few counts. For the simulated mixture of four components noise was added to the value of G and S as shown in figure 3(A). In each panel of figure 3(B), the X axis is the simulated fractional intensity of each component and the Y axis is the recovered fractional intensity using the proposed algorithm. The correlation coefficient of each plot is always above 0.94 and the maximum deviation is less than 10%. Note the maximal deviation is obtained for the fraction F3 as well as the lowest correlation coefficient (0.94). The correlation coefficient as well as the maximum deviation is largely independent on the data noise when noise is added to the data from 0.02 (in terms of coordinates g and s) to 0.16 as shown in figure 3(A). The lack of apparent dependence on noise is due to the way data are simulated; in this case 256×256 pixels are simulated and then the phasor is averaged after the calculation for the whole image. If we had simulated a smaller number of pixels, such as ROI of 16×16, the noise will have become more important.

Experimental confirmation of the simulations using mixtures of pure components
We then constructed a system of these 4 dyes that we mixed to reproduce the conditions used for the simulations. The results are shown in figure 4. The mixture at the center of the 4 green lines was obtained adding 4 components of equal fractional intensity. The individual components were added to the mixture until the solution contained one pure component. These points, which are on the universal semicircle, are indicated with a colored circle in in figure 4.

Result of components analysis for the biological samples
We measured a total of 47 cells. The phasor plot of all cells together is shown in figure 5

Conclusions
In this paper we present a simple algebraic algorithm to resolve the fraction of four components in mixtures of dyes and of 3 components for dyes in cells. The   algorithm is fast. All 47 images were analyzed at the single pixel level in about 50 s. The quantification of relative amount of the dye or any other fluorescent molecules in cells could be used to classify diseases or the response of cells to different stimuli. The distribution of the pixels histograms for each cell is another parameter that can be used to classify cell. Multiplexing several dyes in a sample and determine quantitatively their abundance is a crucial capability for modern diagnostics techniques. The phasor approach in this work adds to the already powerful methods allowed by the phasor analysis in which parts of a sample can be recognized by different lifetimes. This work opens another level of analysis by determining the number and abundance of molecular species in a single pixel.
For the algorithm to work for 4 components one needs data at 2 harmonics, which is a simple matter using the high frequency FLIMbox that has a bandwidth of about 1 Ghz. Harmonics 2 and 3 have very little noise as shown by the resolution of mixtures of 4 components. We note that in addition to the lifetime phasor, one may acquire the spectral phasor for the same sample at the same time, which can also be used to resolve multiple components [32].

Author contributions
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.