Single-Cell Biochemical Multiplexing by Multidimensional Phasor Demixing and Spectral Fluorescence Lifetime Imaging Microscopy

Revealing mechanisms underpinning cell function requires understanding the relationship between different biochemical reactions in living cells. However, our capabilities to monitor more than two biochemical reactions in living cells are limited. Therefore, the development of methods for real-time biochemical multiplexing is of fundamental importance. Here, we show that data acquired with multicolor (mcFLIM) or spectrally resolved (sFLIM) fluorescence lifetime imaging can be conveniently described with multidimensional phasor transforms. We demonstrate a computational framework capable of demixing three Forster resonance energy transfer (FRET) probes and quantifying multiplexed biochemical activities in single living cells. We provide a comparison between mcFLIM and sFLIM suggesting that sFLIM might be advantageous for the future development of heavily multiplexed assays. However, mcFLIM—more readily available with commercial systems—can be applied for the concomitant monitoring of three enzymes in living cells without significant losses.


INTRODUCTION
The fluorescence lifetime (τ) is the average time a fluorescent molecule spends in the excited state before returning to the ground state with the emission of a photon [5,44]. Often, τ depends on the physicochemical characteristics of the environment surrounding the fluorophore but does not depend on the fluorophore concentration; thus, fluorescence lifetime sensing has been applied successfully to probe cell biochemistry (e.g., pH, analyte concentration, enzymatic activities, protein-protein interactions, and conformational changes). Fluorescence lifetime imaging microscopy (FLIM) is commonly used to map cell biochemistry in a quantitative and low-invasive way [18,28,45,46].
Genetically expressed (e.g., using fluorescent proteins) biosensors permit researchers to probe diverse biochemical reactions exploiting the use of Forster resonance energy transfer (FRET). FRET is the non-radiative transfer of energy from a donor fluorophore to an acceptor molecule that can occur when the two molecules are sufficiently close, typically <10 nm [13,37]. FRET causes a loss in fluorescence intensity (quenching) and a decrease of the donor's fluorescence lifetime, both of which are directly proportional to the fraction of energy transferred from the donor to the acceptor. Thus, FRET can be exploited to quantify interactions at the nanometer scale via using a diffraction-limited microscope [75]. This property makes FRET an excellent tool for minimally invasive assays to probe molecular interactions, signal transduction, and biochemical activity in living cells [2,48].
Historically, the estimation of fluorescence lifetime relied on the iterative fitting of the experimental decays [4,31] or simple analytical relations between signals integrated in a few timewindows (e.g., [65]). Over the last decade, however, nonparametric data analysis using frequency-domain techniques [20,33,41], phasor-based representation [10,18,25,35,68], and extended phasors [8] have become very popular in the community. This strategy is simple, computationally efficient, and, more importantly, does not require model assumptions (e.g., mono or bi-exponential) that in complex, multiplexing assays might easily break down (e.g., because of background, cross-talks, and spurious signals). This is especially important at the low photon budget we must operate instrumentations to minimize phototoxicity.
With a large palette of fluorescent proteins that can be used to construct FRET biosensors of different colors (from the UV/ blue to near infrared [3,16,27,40,53,54,66]), we can monitor multiple signaling events in living cells. Dual FRET biosensing was demonstrated with different experimental configurations and demixing analysis frameworks: multichannel ratiometric detection [1,30,57,67], simultaneous homo-FRET and hetero-FRET detected with time-resolved anisotropy and global analysis [74], and dual-color FRET-FLIM to follow two biosensors using time-domain analysis [17]. However, the simultaneous detection of three and more FRET pairs remains challenging. Previously three-protein interaction was measured using "triple fluorophore" three-way FRET sensing [7,29,64,69].
Recently, we have shown that the rational design of FRET pairs aimed to optimize the utilization of the visible spectrum enables multi-color FLIM (mcFLIM) to multiplex at least three FRET probes with a single excitation laser [27] at the speed and resolution necessary for live cell imaging [71]. The utilization of additional excitation wavelengths and hyperspectral detection might extend this strategy to the integration of a higher number of markers, biosensors, or optogenetic tools [39,67,71].
However, computational frameworks for the robust and sensitive multiplexing of FRET are at their infancy, and further work is necessary to improve our multiplexing capabilities. Phasor analysis of time decays has often been applied to quantify FRET for single probes by separating the two states of a typical sensor (low/high FRET) [10,14,18,26,35,36,38,41,47,49,50,55,60,70,72,73]. Similarly, the integration of spectrally resolved FLIM (sFLIM) [34] and multidimensional phasor analysis has been successfully applied to the quantification of single FRET probes [25]. In Fries et al. [27], we have illustrated how to utilize multicolor FLIM and multidimensional phasors for demixing three FRET probes. However, to our knowledge, the demixing of multiple fluorescent species (donor, acceptors, and FRET for several probes simultaneously) has not been sufficiently explored.
In this work, first we briefly illustrate the development of a photon-efficient spectrally resolved FLIM based on off-the-shelf components. Second, we demonstrate the generalization of the NyxSense computational framework which we had first introduced for mcFLIM applications [27]. Most importantly, we provide a detailed comparison of various algorithms aimed to provide efficient dimensionality reduction by multidimensional phasors that can be used for demixing three FRET pairs by spectral FLIM (either multicolor or spectrally resolved FLIM).
We show that phasors of higher dimensionality significantly improve demixing algorithms for both mcFLIM and sFLIM. We demonstrate that while spectrally resolved FLIM could provide advantages in demixing of more than three FRET probes, stateof-the-art fast time-correlated single-photon counting (TCSPC) still attains very high performances with significant implications for future developments of multiplexing time-resolved platforms.

Microscopy
We have developed a simple spectrally resolved FLIM (sFLIM) setup built with off-the-shelf components, including a 16-channel multi-anode GaAsP photomultiplier tube (PML-16-GASP16, Becker&Hick GmbH) placed at the de-scanned port of a SP5 Leica Confocal Microscope (Leica Microsystems United Kingdom, Ltd.). Spectral dispersion was achieved with a direct vision prism (G331120000, LINOS, GmbH) to provide a simple alignment and low optical losses. Notably, GaAsP photomultiplier tubes provide very high quantum efficiencies (45% at 500 nm) compared to the previous generation of photocathodes available (<20% for the bialkali PML-16-1), without compromising the instrument response function of the system significantly (220 vs 200 ps, nominal values provided by the manufacturer). The electrical signals from the photomultiplier assembly were routed to time-correlated singlephoton counting electronics (SPC-150 by Becker & Hick GmbH) utilizing a single arm of the hyperdimensional imaging microscopy electronics we have described previously [19]. Single confocal plane images were acquired with a 40x oil objective (Leica HCX PL APO CS 40.0 × 1.25 OIL UV), 256 × 256 pixel image size, and 120 s acquisition time. A simultaneous two-photon excitation of the FRET pairs was achieved with a Ti: Sapphire Laser Chameleon Vision II (Coherent Inc.) tuned at 860 nm.

Cell Culture And Time-Lapse Imaging
For time-lapse imaging, we used HeLa-CCL2 cell lines (European Collection of Cell Cultures #93021013) expressing the three sensors with the plasmid described in [27]. The sensors have been fully characterized in our former publication: TagBFP-sREACh, mAmetrine-msCP576, and mKeima-tdNirFP, fused with the flexible linkers containing the sequences VDTTD, DEVDR, and LEHD that are cleaved preferentially by caspase 2, caspase 3, and caspase 9, respectively. Cells were periodically mycoplasma-tested and STR profiled using the services of the CRUK Cambridge Institute. Cells were treated with the genotoxic drug cisplatin to induced cell death and imaged for 8 h at 1 h Frontiers in Physics | www.frontiersin.org May 2021 | Volume 9 | Article 637123 2 intervals in LabTek II glass-bottom chambered slides (Nunc, #1.5) containing a 400 µl Leibovitz (L-15) medium supplemented with 10% FCS, 100 μM Cisplatin, and 0.9% NaCl. Multidimensional phasor fingerprinting of individual components was performed just before the time-lapse experiment, with HeLa-CCL2 cells transiently transfected with donor fluorophores not fused with an acceptor, and a donor-acceptor fusion pair (known FRET), as shown in [27].  Figure S1 for green and red sensors) with the spectrally resolved time decay N (τ, λ) integrated over the mask of a single segmented cell is shown (inset). (D) Spectra of the blue (blue solid line), green (green solid line), and red (red solid line) positive controls (fusion constructs) and of the corresponding negative controls (only donors, dashed lines). Each curve was normalized to the total number of counts. The horizontal arrows indicate which spectral bins were used to emulate three-channel (ch.) multicolor FLIM (mcFLIM) data. Spectrally resolved time-decay N (τ, λ) (left panels in (E-H)) and corresponding spectrally resolved phasors (sTP → , right panel in (E-H)) for the blue donor (E,G) and positive control (F,H) for sFLIM (E,F) and emulated three-channel mcFLIM (G,H). Inserts are magnifications of the same graphs. (I) The same data from the blue positive (stars markers and label B) and negative (round marker and label b) controls are shown transformed also as time-(TP), spectral-(SP), and time-spectralphasor (TSP) space. Note, figure (I) is the same for spectrally resolved FLIM (sFLIM) and mcFLIM. In all panels, DCT-discrete cosine transform, DST-discrete sine transform; λ-marks spectral domain or spectral transformation (λDCT/λDST); τ-marks time-domain or time phasor transformation (τDCT/τDST).

Data Analysis
All the analyses were performed with custom-written Matlab script package freely available at https://github.com/inatamara/ NyxSense. Segmentation, tracking, and application of NyxSense to mcFLIM have been described previously [27]. Briefly, cell segmentation was performed on the intensity images (decay curves and spectral channels were integrated) with the combination of an active contour algorithm [9] and a manual curation of the mask. Subsequently, cells were tracked between two consecutive images using a nearest neighbor approach, and mistracked cells were manually reassigned. A spectrally and timeresolved measurement for each cell was then achieved by summing the two-dimensional TCSPC histograms within the cell mask. The latest version of NyxSense (used here) also provides the capability to analyze spectrally resolved FLIM data. The performance of NyxSense for spectrally resolved or multicolor FLIM was evaluated using the same datasets. mcFLIM was generated by binning to the sFLIM spectral channels 1-6, 7-13, and 14-16, for channels 1, 2, and 3, respectively.

Multidimensional Phasor-Based Demixing of Spectrally Resolved Fluorescence Lifetime Imaging Microscopy Data
In the phasor space, single exponential decays are mapped to points on a semicircle described by the equation (x-0.5) 2 + y 2 0.5 2 ( Figure 1A) [12,18]. The time phasor coordinates are defined by the real and imaginary parts of the Fourier transform of the exponential decay function, and calculated as a discrete cosine (DCT) and a discrete sine (DST) transform, respectively. Points lying inside this semicircle correspond to mixed exponentials, being either inherently multi-exponential or a mixture of single-lifetime components. Similarly, all possible spectral phasors lie on arcs bounded by a circle x 2 + y 2 1, resulting from the Fourier transform of pure Gaussian spectra ( Figure 1B) [22,23,23,25]. The two-dimensional time-spectral phasor (TSP) is a twodimensional Fourier transform, in which values are also bounded by a circle x 2 + y 2 1.
The phasor transform key advantage is additivity: a mixture of the spectral or lifetime components corresponding to a linear combination of phasors. This permits rapid demixing using a system of algebraic equations. A point in a phasor space corresponding to the combination of two lifetimes or spectra lies on the line connecting these two pure components. The distance to each pure component along connecting line segments translates directly to its fractional contribution within the mixture. In general, a phasor representing a mixture of n components is enclosed by a polygon with n vertices defined by the phasors of elementary components [23].
For every cell at each time point, we calculated multidimensional phasors: spectrally resolved time phasors (sTPs), spectrally integrated time phasors (TPs), spectral phasors (SPs), and timespectral phasors (TSPs). The time-spectral phasor is a twodimensional transformation along the time dimension followed by the transform along the spectral dimension.
where N TOT τλ N τλ is the total number of photons detected for a given cell, N τλ denotes spectrally resolved time decay, N τ λ N τλ is time decay, i denotes an imaginary unit, and i 2 −1, φ τ 2πnp(S i τ − 1/2)/S τ is a phase for time (τ) phasor computation, where S i τ is the ith time bin, S τ is a number of time bins used to compute the phasor transform (here we used 46 out of 64), and n is a harmonic number.
where N λ λ N τλ is a spectrum and φ λ 2πnp(S i λ − 1/2)/m is a phase for spectral (λ) phasor computation, where S i λ is the ith spectral bin, m is a number of spectral bins used to compute the phasor transform (16 for sFLIM and 3 for mcFLIM), and n is a harmonic number.
The fluorescence signatures of cells or reference samples were then characterized by the multidimensional phasors defined by the complex vector P x → TP SP TSP sTP 1 / sTP m , where the subscript x indicates the multidimensional phasors of the measurement or the reference fingerprints. The subscript m denotes the number of spectral channels for sFLIM (m 16) or mcFLIM (m 3).
Demixing of sFLIM can be achieved by minimization of a complex nonlinear multivariable constrained function (CF) with respect to fractional contributions of the six control signatures (C). At each minimization step, CF is computed as a squared residual between experimental phasors (P exp → ) and phasors estimated using the fractional contributions (P est → ): τλ is a column vector of spectrum and sTP k → is a column vector of spectrally resolved time phasors for the kth control signature. TP est k C k TP k , SP est k C k SP k , and TSP est k C k TSP k , where TP k , SP k , and TSP k are time, spectral, and time-spectral phasors for the kth control signature. However, the minimization of the complex CF renders undetermined system for certain phasor combinations (e.g., P est → [ TP est SP est TSP est ], three equations for six unknown variables, and most of mcFLIM P est → combinations). To assure that the system of equations is not underdetermined, to compare mcFLIM and sFLIM, we used real and imaginary parts of P exp → and P est → separately, that is, and the remainder is as described above. We note that the demixing results using complex CF and real CF (with twice as much equations) are almost the same even for the mcFLIM. In addition, for the minimization involving P est → [ TP est SP est TSP est ] for mcFLIM/sFLIM or P est → [ sTP est → ] for mcFLIM, the phasors were calculated at the first and the second harmonic. This assured that the number of equation is greater than the number of parameters to estimate, which was necessary to calculate the standardized residuals (see Eq. 10).
The minimization procedure was achieved using fmincon Matlab solver. The lower (LB) and upper bonds (UB) for the fractional contributions were constrained to 0 and 1, respectively. The initial values for the fractional contributions were typically 0 for all the control signatures.
The relative enzymatic activity (REA) for each FRET sensor (caspase) was calculated using the following equations: where f d and f unc are the fractional contribution of the donor-only and uncleavable sensor control signatures, respectively, E is FRET efficiency, and division by (1 − E) compensates for the change in brightness. To avoid the division by a very small number leading to large errors, REA was set to 0 for f d and f unc typically below 0.01-0.02. We note that in the specific case of proteolytic sensor, REA represents the cumulative enzymatic activity of the proteases as cleavage is irreversible (until new sensors are expressed de novo).
The standardized phasor residuals were calculated as a difference between experimental phasors and phasors calculated using the unmixed fractional contributions.
where σ is the estimated residual standard deviation and h ii is a leverage of the ith observation (i.e., ith element of the residual vector): σ n (P exp → − P est → ) 2 /(n − p),where n is the number of equations (number of elements in P est → or P exp → ) and p is the number of parameters (six control signatures).
The root mean square deviation (RMSD) for the REA was computed as follows: where P is a FRET pair (B, G, R), N t 8 is the number of the experimental time points, REA 0 is the known enzymatic activity (the ground truth), and REA is obtained from the demixing.

Simulating Spectrally Resolved Fluorescence Lifetime Imaging Microscopy Data
The following equation was used to generate time-and spectrally resolved emission for each FRET pair (n): where e D (λ), e A (λ) are the spectrally resolved normalized emission profiles, τ D , τ A are a lifetime of the donor and acceptor, respectively, and τ DA τ D (1 − E) is a donor lifetime in the presence of acceptor. f D is a fraction of free donors and r 0 is a fraction of the directly excited acceptors. The donor and acceptor absorption cross-sections, quantum efficiencies, the transition rates, and the donor-to-acceptor ratio were set to 1. Finally, ⊗ is a convolution operation. The final counts per pixel was calculated as follows: where N 1 , N 2 , N 3 are the total photon count for each FRET pair, respectively, which was set to 2,600 photons. The lifetime decays were modeled as a single exponential. The donor and acceptor emission spectra (e D (λ), e A (λ)) were modeled as the Lorentzian curves. The synthetic images containing simulated three FRET pairs were generated with the following parameters: the donors' lifetimes (τ D ) were 2.3, 2.5, and 2.7 ns for FRET pairs 1-3, respectively, the acceptors' lifetimes (τ A ) was set to 0.3 ns, and FRET efficiencies 0.35 for each FRET pair. The donors' spectra maxima were 470, 515, and 570 nm and FWHM 55 nm, and the acceptors' spectra maxima were 505, 550, and 605 nm with FWHM 55 nm for each fluorophore. The acceptor direct Frontiers in Physics | www.frontiersin.org May 2021 | Volume 9 | Article 637123 5 excitation relative to the donor excitation was set to 5 or 0%. The Poisson noise was added using the Matlab function imnoise and resulted in ∼12% noise.

Multidimensional Phasor Fingerprint Provides an Efficient Method for Dimensionality Reduction
To test the capabilities of the computational framework presented in the Methods section, we used the NyxBits sensor platform we have described recently to sense cleavage of different substrates (the peptides VDTTD, DEVDR, and LEHD) that are preferentially cleaved by caspase 2, caspase 3, and caspase 9, respectively [27]. In our former work, we demonstrated the capability of mcFLIM to demix the blue (labeled as B in all figures, caspase 2), green (G, caspase 3), and red (R, caspase 9) FRET pairs excited at the same wavelength ( Figure 1). Upon cleavage, each sensor yields two principal components with different lifetimes and spectra: an unquenched donor (labeled with the small b, g, and r letters in all figures) and the uncleaved donor-acceptor undergoing FRET (B, G, and R). Thanks to the large Stokes shifts of the probes and acceptor chromophores of a very low quantum yield, the free acceptors are excited with low efficiency and have a minor impact on our experiments. Here, we characterize the three FRET sensors and compare the performance of mcFLIM and sFLIM using a simple and optically efficient spectrally resolved FLIM (see Materials and Methods).
Each sFLIM image has two spatial dimensions (x, y-here 256 × 256), the time-resolved fluorescence decay histogrammed in 64 time bins (τ) and a spectral dimension represented with 16 spectral bins (λ). The spectrally and time-resolved fluorescence decay of each pixel can therefore be represented in an abstract space of high dimensionality (64 × 16 1,024 numbers, or photon-counts).
Multidimensional phasor transforms permitted us to project this space onto a space of lower dimensionality where the fluorescence characteristics detected in each pixel are described by a vector P exp → [ sTP → TP SP TSP] (see Eqs. 1-8) of 19 complex components (6 for mcFLIM). Although different combinations of phasor transforms have been used previously, here we maintain a higher dimensionality of P exp → than other works [23,52,61] to ensure sufficient features are preserved during dimensionality reduction. Aiming to limit acquisition time and phototoxicity that affect biologically relevant measurements, we have acquired typically 1,000-1,500 photons per pixel. Rather than on pixel basis, we perform cell-based demixing by segmenting and thresholding (pixels with typically minimum ∼200 photons are retained) individual cells and integrating photons collected within each cell. In the time-spectral domain (τλ), the biochemical state of a cell is thus described by a surface spanned by the number of photons (N), spectral information (λ), and time decay (τ) (Figures 1C,D).
The reference phasors were obtained by imaging cells expressing only one control signature, that is, only a donor or a sensor rendered non-cleavable by substituting the substrates with a proteolytically stable sequence. Figures  1E,F show the unquenched blue donor (b) and a blue uncleavable FRET pair (B) fingerprints, respectively, including spectrally resolved lifetime, decays (N[λτ]), spectrally resolved time phasors (sTPs), spectrally integrated time phasors (TPs), time-integrated spectral phasor (SP), and time-spectral phasors (TSPs) (see also Supplementary Figure S1). We compare the biochemical sensitivity of spectral FLIM to multicolor FLIM by binning the 16 spectral bins into three channels that numerically emulate multicolor detection ( Figures  1D,G,I, Supplementary Figure S1, see also Material and Methods). This strategy permitted us to compare the computational performance of the two methods, without having to account for differences in the detection efficiency of two detection systems that would be otherwise difficult to control experimentally.

Multidimensional Phasor-Based Demixing Minimizes Cross-Talks Between Sensors
Subsets of the components of P exp → (sTP → ) have been previously used to demix single FRET pair (donor and acceptor fluorescence, and interacting donor-acceptor pairs [25]). Different subsets of P exp → (TP, SP, and TSP) were applied to separate three fluorescence components using phasor plots [52], and blindly demix three signal components for contrast enhancement in tissue imaging [24,61]. Here, we used the full complement of the features described by P exp → (sTP → with the combination of TP, SP, and TSP) and experimental reference fingerprints to ensure robustness and reproducibility of the results. First, we tested our framework by demixing single FRET pair images containing only two reference components (b-B, g-G, or r-R; Figure 2). This approach permitted to evaluate false-positive detection of the four other components that were not present in a sample. For this, we recorded time-lapse sFLIM images of cells expressing individual sensors (B, G, or R) after exposure to the genotoxic drug cisplatin. Cisplatin induces irreparable DNA damage, leading to switch-like activation of caspases that execute apoptosis. Figure 2A shows that the biochemical trajectories of cells undergoing apoptosis ( Figure 2B) tend to trace a line connecting two control fingerprints (i.e., FRET and no-FRET). The FRET control corresponds to uncleaved sensors; that is, no caspase is activated. During the apoptosis, caspases are activated, sensors get cleaved, and the experimental phasors (black line, Figure 2A) approach no-FRET control phasors signatures. In Figure 3 and Supplemetary Figures S2, S3, we compare spectral demixing using different components of the multidimensional phasors P exp → for both sFLIM and mcFLIM. The single-cell traces of the blueand green-emitting biosensors displayed linear trajectories in a phasor space (Figure 2A). The demixing correctly detected the fractional contribution of control fingerprints (Figure 3, Supplemetary Figures S2, S3), which resulted in ∼80 and ∼50% final sensor cleavage (cumulative relative enzymatic activity (REA)) for a cell expressing the green (G) or the blue (B) sensor, respectively. However, we note that the red sensor (R) can display curved trajectories in the phasor space in several experiments, observation we attribute to non-idealities of the mKeima/tdNirFP, including residual fluorescence of tdNirFP and occasional photo-conversion of mKeima that can occur at higher excitation regimes. However, despite these non-idealities, demixing of the red FRET pair is also sufficiently robust.
For the single FRET pairs (demixing of two components), we observed that the different combinations of sTP → and TP, SP, and TSP components rendered comparable results as by the results shown in Figure 3, Supplementary Figures S2, S3 for both mcFLIM and sFLIM. To compare different demixing methodologies, we provide two figures of merits, the mean standardized phasor residuals, where the mean is taken over observations (i.e., the number of equations) (see Eq. 10, Materials and Methods), either summed over time (Figure 3) or timedependent (Supplementary Figures S2D, S3D). From Figure 3C, we see that the largest fitting errors for the sFLIM occur with the demixing using the smallest P exp → subset (TP, SP, and TSP) (see also Supplementary Figure S2D). Demixing with different combinations of sTP → with TP, SP, and TSP rendered similar results. In comparison, mcFLIM showed higher residuals for sTP → and SP than for sTP → or TP, SP, and TSP alone. Yet, when sTP → and SP demixing was calculated with the first and the second harmonic combined (as it is for sTP → and (TP, SP, and TSP)), the residuals were lower (data not shown). However, to avoid errors caused by overfitting, we utilized the second harmonics when essential to have a determined system of equations (see also Material and Methods). For mcFLIM sTP → or (sTP → , TP) alone, the wrongly detected components were not present in a sample ( Supplementary Figures S3A-C). We conclude that the representation of data with a multidimensional phasor P exp → improves the overall performance of demixing algorithms for both spectrally resolved and multicolor FLIM even for the single FRET pair demixing.

Multidimensional Phasor-Based Demixing Achieves Efficient Multiplexing of Forster Resonance Energy Transfer Biosensors
Next, we validated the proposed methodology by summing photon counts recorded from cells expressing the individual FRET pairs. This strategy permitted us to test the demixing algorithms in the presence of typical cross-talks, providing an experimentally valid ground truth (the known fractional contributions, calculated in Figure 3). For example, Figure 4 shows the same data displayed in Figure 3, Frontiers in Physics | www.frontiersin.org May 2021 | Volume 9 | Article 637123 8 summed, and then demixed either using sFLIM or mcFLIM. With this validation dataset and the analysis of the residuals, we observe that both spectrally resolved and multicolor FLIM provide efficient and comparable demixing, see also Supplementary Figure S4. Most importantly, the residual analysis of both mcFLIM and sFLIM showed that the smallest fitting errors were observed for the largest sTP → complement sTP → , TP, SP, and TSP and the highest for TP, SP, and TSP and sTP → or sTP → and TP. For both mcFLIM and sFLIM sTP → or sTP → and TP alone were not able to detect the small rise in REA R between 5 and 6 h (Supplementary Figure S4). We note that this rise is not an artifact as a clear apoptotic phenotype is seen at 5 h ( Figure 2B, red cell). As phasor residuals do not determine the accuracy in the estimation of relative enzymatic activities, we also estimated the root mean square deviation (RMSD) of the REA values when a ground truth can be estimated (see Eq. 11, Materials and Methods) by subtracting REA obtained with a single FRET pair denoted as REA 0 (Figure 3) from the REA obtained with triple FRET demixing (REA, Figure 4B). In Figure 4D, we show that the RMSD of sFLIM and mcFLIM is very similar.
Next, we applied the proposed methodology to the analysis of the experimental data with cells co-expressing all three sensors. Figure 5A shows the time evolution of biochemical traces of two cells (Cell-1 and 2) that exemplify the different responses we have described previously [27]. Cell-1 exhibits a robust activation of all three caspases, while Cell-2 does not. Figures 5B, C show the demixing results using sTP → , TP, SP, and TSP phasors (see also Supplementary Figure S5) for sFLIM and mcFLIM, respectively. The comparison between the 16 spectral-channel sFLIM and the three-channel mcFLIM demixing suggests that both modalities can efficiently retrieve the six signatures we analyzed. As for the semisynthetic data shown in Figure 4, the addition of spectral phasors and, more generally, the use of the multidimensional phasors P exp → [ sTP → TP SP TSP] improve fluorescence lifetimebased multiplexing of multiple FRET pairs. The residual analysis showed that sTP → , TP, SP, and TSP produce the smallest and TP, SP, and TSP and sTP → or sTP → and TP produce the highest fitting errors. In addition, with multicolor probes and detection channels optimized for FRET multiplexing, mcFLIM performs almost as well as spectrally resolved FLIM. However, comparing red traces in Figures 4B,C, 5B,C, as well as Supplementary  Figures S4A,B, S5A,B, suggests that sFLIM may be more robust for low fractional contributions of the individual components (typically below 2-5%). We conclude that the efficient demixing of three FRET pairs required the combination of sTP → and spectral phasors SP and TSP for both mcFLIM and sFLIM and that, more generally, P exp → [ sTP → TP SP TSP] provides robust demixing. Multidimensional Phasor-Based Demixing is Necessary for Demultiplexing High Cross-Talk Data Taken together, the results shown suggest that multicolor FLIM is already an efficient technique for demix of three FRET pairs when analyzed with multidimensional phasors although it requires optimized FRET pairs. With the engineering of faster and efficient spectral FLIM system, we envisage that sFLIM might provide a significant advantage for demixing, for example, with FRET pairs not specifically designed for heavily multiplexed detection. We therefore investigated how more significant spectral overlaps might affect the performance of mcFLIM and sFLIM. We generated a fully synthetic triple FRET pair images (see  Figure S6). Further increase in the overlap between the green and the blue FRET pair still led to correct demixing with sFLIM and mcFLIM, Supplementary Figure S7.
Once again, multidimensional phasors provided more robust demixing performing optimally in all the tested conditions. With low-to-medium spectral overlap, mcFLIM performs almost as

DISCUSSION
Time-or spectrally resolved imaging has become very accessible, thanks to commercial confocal microscopes frequently installed in core facilities that can support these applications. Such technologies offer opportunities to match the increasing demands for enhanced multiplexing of fluorescent markers. The quantitative characterization of biochemical network dynamics is an invaluable application of multiplexed fluorescence to formulate or test hypotheses at the interface between cell and systems biology. Nevertheless, this application is comparatively immature, at least concerning more quantitative approaches extended to small biochemical networks. Demultiplexing several biochemical activities from complex photophysical datasets is still a challenge that requires breaking new grounds to enable the robust characterization of biochemical networks at single-cell resolution [7,15,30,62]. This challenge is made harder by the nonideal conditions necessary to image living cells. For example, the need to minimize photon toxicity leads to low photon counts and spurious signals. In addition, the use of several fluorescent proteins that are determined by the need for multiplexing constraints (spacing between excitation/emission spectra) might exacerbate issues such as improper maturation of acceptor chromophores and brightness, photobleaching, and photochromism. Therefore, here we contribute to this endeavor with a framework based on multidimensional phasor transforms, representing efficient and intuitive methods for demixing three FRET pairs excited at a single wavelength. Building on work published by us and others (e.g., [21-23, 25, 27, 61]) we extended this computational framework, including higher phasor dimensionality. We demonstrate its efficacy to demix three FRET pairs imaged at a single excitation wavelength that we previously optimized for the multiplexing and by simulating synthetic triple FRET pair images with more spectral overlap. We provide a description of our methodology and the extension to spectrally resolved FLIM of the NyxSense computational platform that we had briefly described only for mcFLIM [27]. This platform is available in the public domain (https://github. com/inatamara/NyxSense) and could be used by the community to test, further improve, or simply use the methodology we proposed. Spectrally resolved FLIM is readily available commercially, and several bespoke implementations aimed to make available cost-effective and user-friendly solutions have been also published (e.g., [6,42,43,51,56,58,61,63,76], promising increased availability of such sophisticated assays in the near future. We showed that the combination of spectrally resolved time phasors (sTP → ) with the spectral phasors (SP or TSP) permitted efficient demixing of three FRET pairs, presenting a low level of direct acceptor excitation using only six control signatures. Interestingly, the results for sFLIM (16 spectral bins) and mcFLIM (three spectral channels) were very similar for the data discussed here, that is, six main control signatures (donoronly and uncleavable sensors) with a low level of direct acceptor excitation. However, the sFLIM is more robust to unmix lower fractional contributions. Spectrally resolved FLIM will therefore be an essential tool either to demix common FRET pairs with large spectral overlaps or to further expand our capability to multiplex more than three biochemical reactions from single living cells.
However, we suggest that readily available equipment dedicated to multicolor FLIM, particularly instruments capable of fast detection, can already perform such complex experiments efficiently. Therefore, the innovation of detection technologies of both scanning and wide-field microscopes can make biochemical multiplexing a routine technique in the future [11,32,47,59,60]. We show that excellent demixing results can be achieved with the open-source toolbox NyxSense for both sFLIM and mcFLIM. NyxSense implements the multidimensional phasor transforms that facilitate the projection of complex multidimensional photophysical data onto biochemical spaces of lower dimensionality to represent the biochemical trajectory of single cells in response to stimuli.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.