Systematic Enzyme Mapping of Cellular Metabolism by Phasor-Analyzed Label-Free NAD(P)H Fluorescence Lifetime Imaging

In the past years, cellular metabolism of the immune system experienced a revival, as it has become clear that it is not merely responsible for the cellular energy supply, but also impacts on many signaling pathways and, thus, on diverse cellular functions. Label-free fluorescence lifetime imaging of the ubiquitous coenzymes NADH and NADPH (NAD(P)H-FLIM) makes it possible to monitor cellular metabolism in living cells and tissues and has already been applied to study metabolic changes both under physiologic and pathologic conditions. However, due to the complex distribution of NAD(P)H-dependent enzymes in cells, whose distribution continuously changes over time, a thorough interpretation of NAD(P)H-FLIM results, in particular, resolving the contribution of various enzymes to the overall metabolic activity, remains challenging. We developed a systematic framework based on angle similarities of the phase vectors and their length to analyze NAD(P)H-FLIM data of cells and tissues based on a generally valid reference system of highly abundant NAD(P)H-dependent enzymes in cells. By using our analysis framework, we retrieve information not only about the overall metabolic activity, i.e., the fraction of free to enzyme-bound NAD(P)H, but also identified the enzymes predominantly active within the sample at a certain time point with subcellular resolution. We verified the performance of the approach by applying NAD(P)H-FLIM on a stromal-like cell line and identified a different group of enzymes that were active in the cell nuclei as compared to the cytoplasm. As the systematic phasor-based analysis framework of label-free NAD(P)H-FLIM can be applied both in vitro and in vivo, it retains the unique power to enable dynamic enzyme-based metabolic investigations, at subcellular resolution, in genuine environments.


Introduction
The importance of cellular metabolism in immunology has become increasingly recognized in the last years. Although, in the last decades, molecular immunology mainly focused on signaling pathways within various immune cell subtypes and their supporting cells, such as stroma cells, it has become evident that cellular metabolism does not merely provide energy necessary for transcription and biosynthesis in immune cells, but is tightly interconnected with their function [1,2]. For instance, the process of selection and further differentiation of germinal B cells was shown to be directly linked to cellular metabolism [3]. These correlations between metabolism and immune cell function and dysfunction still require extensive investigations, as oxygen consumption measurements related to enzymatic activity are mainly performed ex vivo and not at the single-cell level. Only by using synthetic probes, such as dichlortris(1,10-phenanthroline) ruthenium(II) hydrate (Ru(Phen)) [4], the oxygen concentration in cells and tissues can be measured at the level of single cells [5]. However, these approaches allow only an indirect link to enzymatic (metabolic) activity.
Fluorescence lifetime imaging (FLIM) was developed three decades ago, and since then it has been applied to study cellular function in a quantitative manner [6]. FLIM generates images in which contrast is obtained by the excited-state lifetime τ of fluorophores instead of their intensity, therefore, having negligible experimental bias. The combination of FLIM with two-photon microscopy-the technology of choice when performing dynamic deep tissue imaging in living mammals [7]-allowed it to reach its full potential under genuine conditions and to study tissue and cellular function and dysfunction in vivo. In this way, using FLIM, molecular mechanisms underlying neuronal damage and oxidative stress have been understood in chronic neuroinflammation [8][9][10][11][12], neurodegeneration [13], and glioblastoma [14,15]. Furthermore, the shift from differentiation towards uncontrolled proliferation in cancer cells has been elucidated in several primary and metastatic tumor cells [16], and the physiologic process of differentiation has been unraveled in cells and organs [17][18][19][20][21].
A multitude of approaches have been developed over the years to measure fluorescence lifetimes in bulk solution, without spatial resolution. Later on, they have been applied to microscopy to perform FLIM. Both frequency-domain and time-domain methods have been developed to meet different needs in biosciences, biomedicine and medicine [22][23][24][25][26][27][28]. Frequency-domain technologies are particularly adequate for fast, high-throughput measurements [22,29,30], as well as for retrieving multiple differing fluorescence lifetimes from complex samples using several modulation frequencies, ranging from kHz to MHz [31,32]. Among the time-domain technologies, time-correlated single-photon counting (TCSPC), which requires pulsed excitation as delivered by two-photon microcopy, is the method of choice to comprehensively acquire the molecular complexity within living cells in deep tissue, despite its rather slow speed (1-10 s/frame) [11,33]. TCSPC directly measures the fluorescence decay of all contained fluorophores; however, its thorough analysis remains a challenge especially in complex microenvironments as found in cells.
The fluorophores first imaged by FLIM were the coenzymes nicotinamide adenine dinucleotide-NADH and nicotinamide adenine dinucleotide phosphate-NADPH [6] (hereafter NAD(P)H), i.e., ubiquitous coenzymes governing energy production in all cells across species, as well as being responsible for various biosynthetic processes and influencing diverse signaling pathways [34]. NAD(P)H-FLIM [35] has been extensively employed in cancer research, revealing shorter NAD(P)H fluorescence lifetimes in tumor cells as compared to controls, presumably due to a shift from oxidative phosphorylation (OxPhos) towards glycolysis [16,[36][37][38][39][40][41]. By enabling the differentiation between the two major cellular metabolic pathways, namely, glycolysis and OxPhos, it has also revealed the tight link between metabolism and differentiation in other cell types, for example, in neuronal stem cells [18], in mesenchymal stromal cells during their differentiation towards various fates such as adipocytes, in chondrocytes, or osteoblasts [21,[42][43][44], and in enterocytes in the small intestine [17]. Additionally, by dynamically monitoring the activation of NADPH oxidases, this technique has been applied to monitor the oxidative burst and massive oxidative stress production during infection in the small intestine [45], as well as apoptosis in cells [46]. Both fluorescence lifetime spectroscopy and FLIM are able to differentiate between the states of free and enzyme-bound coenzymes NADH and NADPH. In the free state (unbound to enzymes), the fluorescence of the adenine is intramolecularly quenched by nicotinamide, resulting in a short fluorescence lifetime of around 450 ps [29,47]. Upon enzymatic binding, this quenching is prevented due to steric hindrance caused by the specific binding site of the coenzyme, resulting in enzyme-specific prolongation of the fluorescence lifetime towards 2000 ps [48]. Additionally, the fluorescence lifetime of both free and enzyme-bound NAD(P)H, as with the lifetimes of the majority of fluorophores, is also influenced by other factors such as refractive index [49], solvent polarity [29], pH [50,51], ion concentration [52], or the degree of freedom for the diffusional rotation of the fluorescent molecules. Gregorio Weber in 1970, and then Joseph R. Lakowicz in 1992, showed that there are differences in the fluorescence lifetime of NADH and NADPH depending on the binding site of the enzyme [6,29]. This principle has been used for a number of isolated enzymes. However, in living cells, multiple pathways employing NAD(P)H-dependent enzymes are active at any given time, making it hard to interpret NAD(P)H-FLIM data derived from living cells with respect to the proportionate usage of certain enzymatic reactions. No framework allowing a systematic interpretation of the NAD(P)H-FLIM data with respect to both metabolic activity and specific, predominant enzymatic activity is available yet. As such, there are apparently contradictory interpretations, indicating that the fluorescence lifetime of NAD(P)H bound to enzymes depends on the coenzyme alone [53] is caused by other molecular species, e.g., oxidized lipids [54], or has a fixed value of~2000 ps [23,39,48,55,56] or~3400 ps [54,57].
Here, we established a general framework retrieving two kinds of information from NAD(P)H-FLIM data: (i) the fraction of free vs. enzyme-bound NAD(P)H present in cells, i.e., a measure of immediate metabolic activity, and (ii) which enzymes are the major contributors to the NAD(P)H-dependent metabolic activity in a certain cellular compartment, at a particular time point. Our approach is based on the measurement of the fluorescence lifetime of NADH and NADPH free and bound to enzymes, which are highly and ubiquitously expressed by mammalian cells. We performed the measurements of homogeneous mixtures of NAD(P)H and single enzymes in buffered media resembling the pH, ion concentration, solvent polarity, and refractive index of cells to ensure that only the binding state of the coenzymes affected the fluorescence lifetime. Typical images acquired by NAD(P)H-FLIM on cells and tissues have a voxel size of 500 nm × 500 nm × 1500 nm. Considering the mean concentration of proteins in a cell [58], we calculated that a number in the range 10 8 -10 9 protein molecules are present within a voxel. Of these proteins, 18% are metabolic enzymes [58]. Even if the number of enzyme molecules may vary between different cell populations, a large number of enzymes are competing to bind NAD(P)H. Therefore, we expect that only the activity of highly abundant NAD(P)H-dependent enzymes will have a major impact on the NAD(P)H fluorescence lifetime. Unpublished RNA-Seq data from our colleagues showed that out of 191 NAD(P)H-dependent enzymes expressed by mesenchymal stromal cells, 16 ubiquitously expressed metabolic enzymes, their including isoforms are contained within the 50 most abundantly expressed genes. These findings are in line with other studies, i.e., cancer cells [58], and presumably hold true for many other cell types. From these highly expressed metabolic enzymes we only excluded complex I, due to difficulties in maintaining the activity of its flavin mononucleotide-binding domain in solution. The NADH-binding unit of complex I is one of the less abundant among these 16 enzymes. We focused on these highly abundant NAD(P)H-dependent metabolic enzymes (covering the metabolic enzymes but excluding their isoforms), i.e., malate dehydrogenase (MDH), lactate dehydrogenase (LDH), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), glucose-6-phosphate dehydrogenase (G6PDH), pyruvate dehydrogenase (PDH), alcohol dehydrogenase (ADH), C-terminal binding protein 1 (CTBP1) [57], and isocitrate dehydrogenase (IDH); complemented by the still abundant enzymes, i.e., hydroxyacyl-coenzyme-A dehydrogenase (HADH) and sorbitol dehydrogenase (SDH); and enzymes playing distinct roles in the cells, i.e., inducible nitric oxide synthase (iNOS) and the NADPH oxidases family (NOX1-4 and DUOX1,2) [59], which play a major role in the catalysis of oxidative burst. As the quaternary structure of the NAD(P)H binding site and the related catalytic function in these NAD(P)H-dependent enzymes is highly conserved across species and in various cellular organelles, it was sufficient to investigate only one variant for each enzyme.
Our framework to evaluate NAD(P)H-FLIM data by the phasor approach, i.e., a model-free analysis tool of FLIM data [59,60], additionally uses for normalization the signal-to-noise ratio of the NAD(P)H intensity image for signal quality control. It decouples the metabolic activity level, i.e., the fraction of free to enzyme-bound NAD(P)H, from the type of predominant enzymatic activity, as indicated by the angle similarity to the phase vector of NAD(P)H bound to pure enzymes. We validated the performance of our analysis framework on homogeneous mixtures of NADH and enzymes as well as on a stromal-like 3T3-L1 cell line. Our systematic framework meets the increasing demand to interpret in a reliable manner, the metabolism of immune and stromal cells, at subcellular resolution, by NAD(P)H-FLIM, both in cell culture [61] and in living tissues and organisms [8,9,13].

Benchmarking of NAD(P)H-FLIM Data Evaluated Using the Phasor Approach
To define the image and signal properties to ensure reliable data interpretation, we performed NADH-FLIM measurements in standardized NADH solutions of various concentrations, as well as in NADH mixtures with metabolic enzymes, i.e., malate dehydrogenase (MDH) or lactate dehydrogenase (LDH).
We evaluated the FLIM data acquired in time-domain by time-correlated single-photon counting in a two-photon microscope using the phasor approach. The raw data consist of a fluorescence decay curve of the fluorophore, in our case NADH, in each pixel of the image. The decay curve is monoexponential if all fluorophore molecules in the volume corresponding to a pixel in the image have the same molecular environment. If the molecules have various molecular environments, the decay curve became multiexponential. The characterizing parameter of the fluorescence decay curve is the fluorescence lifetime, i.e., the time period the fluorescence intensity decreases to 1/e (where e is Euler's number) of its maximum value immediately after the excitation pulse. In the phasor plot, the lifetimes follow a naturally logarithmic distribution on the half-circle (Figure 1(Ai)). After adding the LDH inhibitor, FX11, the experimental data cloud migrates back to free NADH with increasing FX11 concentration. The asterisks indicate the data sets used to calculate the SNR ( Figure S3) as well as calculate τ 2 to mark the half circle. The cyan lines represent the trajectories on which the data cloud migrates if the enzyme activity increases.

Low Signal-to-Noise Ratio of the Image Affects the NAD(P)H-FLIM Data
Both mono-and multiexponential decays are convolved with shut noise, i.e., an undamped multifrequency oscillation, which is also transformed to the phase domain and lies within the area defined by the half circle. To quantify the role of shut noise, we first investigated the impact of the signal-to-noise ratio of the image on the results of the phasor-evaluated NAD(P)H-FLIM data by varying either the concentration of NADH while the laser power stayed the same at 100 mW ( Figure 1A) or by varying the excitation laser power while the NADH concentration was fixed at 50 µM (Video S1 and Figure S1).
The signal-to-noise ratio (SNR) is defined as the ratio of the mean background-free fluorescence signal (µ signal − µ BG ) and the noise of the background (shut noise), i.e., the standard deviation of the background intensity distribution (σ BG ): An increasing NADH concentration results in increasing signal-to-noise-ratio (SNR), identifiable by decreasing width of the phasor cloud, i.e., the distribution of the phase vectors measured in each pixel of the image, in the phasor plots. Figure 1(Av) shows the SNR vs. the vector length, which is the ratio of the vector "(0|0) to free NAD(P)H position onto the half circle" (green line in Figure 1(Aii-v)) and the vector connecting (0|0) with the center of the phasor cloud, in percent. Increasing the laser power while the NADH concentration is fixed shows the same effect.
Independently of the way we modified the SNR in the image, we found that for SNR values above 5, the quality of the fluorescence signal delivered results very similar to those predicted by the theory of the phasor approach. No experimental data are able to reach the accuracy of theoretically calculated values due to the shut noise contained therein. In the range of SNR values less than 5, the vector length, i.e., Re 2 + Im 2 (where Re is the real part and Im the imaginary part of the complex number), decreases linearly with decreasing SNR, whereas the direction of the vector given by the angle to the abscise remains constant.

Validation of Phasor-Evaluated FLIM Data on Mixtures of NADH and MDH or LDH
The theory of the phasor approach predicts that, for a mixture of fluorophore molecules experiencing two different molecular environments, the vector will point onto the segment connecting the positions on the half circle corresponding to fluorophore molecules experiencing only one of these molecular environments. Therefore, the phase vectors of a FLIM measurement of a mixture of NADH and malate dehydrogenase (MDH) point onto the segment connecting the position of free NADH onto the half circle and that of NADH completely bound to MDH.
As the NADH:MDH mixture was a homogeneous solution containing only two fluorescent species, the time domain decay curve was averaged over all pixels in order to smooth it for proper biexponential fitting (least-square sum method). This fit results in τ 1 = 450 ps and τ 2 = 1240 ps (R 2 = 0.99956). τ 1 is the lifetime of free NADH and is in good agreement to known data [47]. τ 2 indicates the lifetime of NADH bound to MDH and was transferred to the phase domain to mark the position of MDH onto the half circle ( Figure 1(Bi-iii)).
In Figure 1(Ci-iv), we investigated mixtures of NADH and MDH with various relative concentrations of coenzyme and enzyme given by the ratio MDH:NADH. We performed the FLIM experiments in such a way that the SNR remained the same for all measurements. SNR was calculated for panel (iii) (marked by asterisk) and was 39.6 ( Figure S3). As expected, the position of the phasor cloud migrates from the free NADH position, towards the position of NADH fully bound to MDH with increasing MDH:NADH ratio, along the connecting vector. FLIM, especially when transformed to phase domain, is very sensitive to noise. Noise has an infinite fluorescence lifetime, which is positioned around the origin (0|0) in the phasor plot. Although our SNR is very high, in the case of MDH and LDH, respectively, noise still represents a third "lifetime component". For this reason, the phasor cloud does not lie exactly on the connecting line, but migrates slightly parallel to it. In the following analysis (especially later in the enzyme assignment), this circumstance was taken into account by the fact that the vector free-enzyme starts at the center of the measured phasor cloud and does not start from the marked position on the half circle ( Figure 3A). Figure 1(Cv) shows the MDH:NADH ratio vs. the vector length, which is the ratio of the vector "free NADH to the center of the center of the phasor cloud" to "free NADH to NADH fully bound to MDH position on the half circle" (red line), in percent.
NADH-FLIM experiments on mixtures of NADH and lactate dehydrogenase (LDH) performed at a SNR far above 5 generally revealed the same valid behavior. Due to the low LDH concentration (10 µM, LDH:NADH ratio = 0.2), the phasor cloud did not reach the position corresponding to NADH fully bound to LDH. By adding the LDH inhibitor FX11 at various concentrations, we could reverse the binding of NADH to LDH, which resulted into a shift back to free NADH in the phasor plot ( Figure 1D and Figure S2). The smearing of the experimental data at high FX11 concentrations in panel (Figure 1(Dv)) results from enzyme denaturation since the intensity image becomes in this case heterogeneous. The SNR remains constant at 22.3 ( Figure S3).

Enzyme-Based Reference System to Interpret Label-Free NAD(P)H-FLIM Data in Cells
The basic concept of our approach for interpreting NAD(P)H-FLIM data takes into account the most likely states in which NADH and NADPH can be found in cells and tissues. Except for the free state, we find the coenzymes bound to diverse enzymes, in a mixture containing~350 NADHand 300 NADPH-dependent enzymes. However, many of these enzymes are so rare that they have a negligible effect on the overall fluorescence decay of the coenzymes. We identified from RNA-Seq data the most abundant and, thus, relevant NADH and NADPH-dependent enzymes and completed the list with NADPH oxidases, i.e., enzymes responsible for oxidative burst generation in certain cell types. The phasors of free NAD(P)H and fully bound NAD(P)H to these enzymes, respectively, are depicted in Figure 2A. The corresponding fluorescence lifetimes are listed in Figure 2B. Using FLIM, we measured mixtures of NADH or NADPH and several of these pure enzymes ( Figure 2C) and evaluated the data using the phasor approach as described in Figure 1B for the NADH/MDH mixtures. The list is completed by our previously published data on other relevant enzymes measured extracellularly in solution [8]. The only exception represents the NADPH Oxidases family (here NOX), which has a fluorescence lifetime and position in the phasor plot and was determined intracellularly using chemical inhibition, activation, and knockout strategies [33,55]. We paid particular attention to perform all FLIM measurements on NAD(P)H solutions and NAD(P)H-enzyme mixtures in buffered media resembling similar pH, ion concentrations, and refractive index as the cellular environment to avoid lifetime artifacts caused by these parameters. As shown by the good agreement of our results of NADH/MDH and NADH/LDH mixtures with the findings of other groups [47] (red lines in Figure 2C), our reference system is a generally valid system that can be applied to any cell or tissue type. The fact that we were not able to measure NADH bound to functional complex I, neither extracellularly nor intracellularly, is a limitation of our reference system.  Figure 1B, fluorescence lifetimes calculated from single images represented as gray dots, average values as black lines, standard error of the mean (s.e.m.) values over all pixels in all measured images per coenzyme/enzyme mixture as black boxes. The list was completed by our previously measured fluorescence lifetimes [6] measured for NADH and NADPH bound to SDH, NADH bound to iNOS, NADH bound to ADH, and NADPH bound to enzymes of the NOX family (red lines). For LDH and MDH the literature data of other groups are indicated as red lines. The color coding in panels (A-C) is consistent.

Analysis Framework of NAD(P)H-FLIM Data: Validation on Homogeneous NAD(P)H-Enzyme Mixtures
Based on the previously generated reference system of coenzymes bound to pure enzymes, we created a systematic analysis framework of NAD(P)H-FLIM data. First, we determined the SNR of the NAD(P)H fluorescence intensity image, which is a quality criterion of the time-domain signal. After calculating the phasor in each pixel of the image, we selected the data points that were located within the circular area in the phasor plot identified as "free NAD(P)H" area ( Figure 3(Ai)). The radius r and position of this circular area was determined from the measured distribution of free NADH ( Figure  S4). We attribute those data points to "free NAD(P)H" and consider that a cell would not show any metabolic activity in those image areas where these data points originate from. To assign the data point to enzymes of our reference system, we calculated the angles α i and α data (Figure 3(Aii)) in the remaining pixels. α i is the angle of the vector connecting the position of measured free NAD(P)H and the positions of NAD(P)H fully bound to the pure enzyme of our reference system with the index i (e.g., i = 2 for MDH, represented as dark green, solid line), and α data is the angle of the vector connecting the position measured free NAD(P)H and each experimental data point (violet, dashed line) (Figure 3(Aii,Aiii)). Furthermore, we calculated the absolute value (expressed by |. . .|) of angle similarity q i to each enzyme of the reference system as follows, Further, we normalized the angle similarities q i for all enzymes in each pixel and obtained the weights w i of the angle similarities q i as follows, Equation (5) gives the assignment probability for each enzyme of our reference system. After multiplying by 100, it may reach from 0 to 100 and the data point is assigned to the enzyme for which this value reaches its maximum. Repeating this procedure for each pixel of an image generated an enzyme map as shown in (Figure 3(Ci-iii)). At the same time, all weights w i may be retrieved as a 3D matrix (2D × 13), as depicted in Figure 3(Fii), for an exemplary 20 × 1 pixel line. In the case that the maximum weight w i , on the basis of which the enzyme assignment is performed, is similar to the neighboring enzymes, the comparison of the weights in the 3D matrix becomes relevant for data interpretation.
In Figure 3B-D, we validated our analysis framework on homogeneous mixtures of free NADH (i), NADH and MDH (ii), and of NADH and PDH (iii), respectively. The enzyme maps ( Figure 3C) of free NADH, as well as NADH mixed with MDH, show a high accuracy of the analysis framework, substantiated by the histograms shown in Figure 3D. The histograms show that 99.9% of the pixels are correctly assigned to free NADH (Figure 3(B,Ci)) and 97.7% of the pixels in Figure 3(B,Cii) are correctly assigned to MDH (enzyme index = 2). The SNR of free NADH image was 39.6 and for MDH:NADH mixture image 22.3. In contrast, the histograms of the enzyme map of NADH mixed with PDH showed that 37.4% of the data points were assigned to ADH, 35.1% to iNOS and only 17.4% to PDH/CTBP1. This is partially due to the high similarity of the lifetimes of NADH if bound to PDH/CTBP1, iNOS, and ADH (see Figure 2C), which makes an accurate assignment difficult. In such a case, the expression level of these enzymes and their subcellular distribution need to be taken into account. Furthermore, the SNR in the PDH:NADH mixture image was with a value of 9.99 lower than that for free NADH and MDH:NADH mixture. As becoming evident from the broad phasor cloud in the corresponding phasor plot (Figure 3(Biii)), in a pixel-based evaluation the SNR criterion needs to be set higher for a sufficient accuracy.
The island-shaped structures in the enzyme maps, especially in the case of the mixture NADH and PDH, do not result from a badly solved enzyme solution, but rather originates from the Gaussian blurring of the time-resolved fluorescence data; one of the very first steps of the analysis. By spatially blurring the time domain raw data, the time resolution is improved at the cost of spatial information. The Gaussian blur with a radius of two pixel (σ = 2) is a good compromise between time and space accuracy ( Figure S5). In the case of homogeneous solutions, the spatial resolution is irrelevant, but this becomes important when performing label-free NAD(P)H-FLIM of cells and tissues (Figure 4). In this case, we additionally mapped the ratio length of the vectors we used to calculate angle α i and α data as defined above. This ratio represents the fraction of free NAD(P)H vs. bound NAD(P)H to the assigned enzyme in each pixel of the image, and thus is a hint for the metabolic activity in this area of the cell.

Validation of the NAD(P)H-FLIM Analysis Framework on Stromal-Like 3T3-L1 Cells
We applied our systematic analysis framework on NAD(P)H-FLIM data acquired in 3T3-L1 cells ( Figure 4A)-a cell line which has mesenchymal stromal cell-like properties and can differentiate into adipocytes when appropriately stimulated. In the bone marrow, stromal cells, the heterogeneity of which is still poorly understood, are known to fulfill various functions. Particularly, in the deep marrow cavity of long bones, they, together with the vasculature, form crucial components of the survival niche for various immune cell populations such as long-lived plasma cells [61], as well as for hematopoietic stem cells [62]. Both in homeostasis and during bone regeneration, stromal cells may differentiate towards osteoblasts, chondrocytes or adipocytes. In different phases of life, the balance between these differentiation pathways and the nondifferentiated state is expected to shift, e.g., in aged individuals, or due to metabolic syndrome, differentiation towards adipocytes is more dominant [63]. The various differentiation pathways have been previously linked to differences in the NAD(P)H metabolism as measured by FLIM [20,42,47].
In particular, differences between the mean NAD(P)H fluorescence lifetime in the nucleus vs. cytoplasm have been reported. However, it is not yet clear which metabolic mechanisms underline these observations.
The SNR in the image ranged between 28 ± 6 in nuclei and 48 ± 16 in cytoplasm, thus our measurements lies within a reliable range. Note that these SNR values could be only reached by summing up subsequently acquired five time-domain images of the cells. The NAD(P)H fluorescence signal is dim in cells and an increase in excitation power would entail the risk of photodamage [64]. Our systematic framework for NAD(P)H-FLIM in 3T3-L1 cells revealed similar metabolic activity as substantiated by the distance maps in Figure 4(B,Cii), i.e., similar fraction of free NAD(P)H. Although there was a certain degree of heterogeneity among the metabolic activity, we observed no significant differences between cytoplasm and nucleus in this respect. In contrast, we found, in the enzyme map, predominant activity of different NAD(P)H-dependent enzymes in nuclei as compared to the surrounding cytoplasm (Figure 4(B,Ciii)). These data are also supported by the corresponding weights in the table in Figure 4(Dii). They reveal that different enzymes from a group of enzymes with presumably similar NAD(P)H binding sites are active in nuclei as compared to cytoplasm. G6PDH was identified as the mainly predominant active enzyme in the nuclei, the probabilities of enzymatic activity of LDH, SHD if bound to NADH and of GAPDH are at similar levels, due to the high similarity of the phasors of the enzymes. We expected to find CTBP1 as a NAD(P)H-dependent enzyme in the nucleus [57], but only found evidence for the presence of LDH, G6PDH, GAPDH, and SDH in there. As these are cytosolic enzymes, we expect the yet uncharacterized nuclear enzymes to have a similar fluorescence lifetime. In the cytoplasm, we found a high predominance of PDH, but also of iNOS and ADH. The contribution of CTBP1 (having a phasor identical to that of PDH) is rather improbable, as this is to the best of our knowledge mainly located in the cellular nucleus [57].

Discussion
Cellular metabolism has a strong impact on cell functions-this holds true for the majority of cell types, across species. In immunology, the need for appropriate methods for the analysis of metabolism has increased since it became evident that there is a strong cross-link between signaling pathways and the cellular metabolism, beyond the mere supply of energy. For example, the activation of B and T lymphocytes coincides with metabolic reprogramming towards a preferential usage of glycolysis [1,65,66], whereas their transition into memory cells goes along with an increase in oxidative phosphorylation. In the bone marrow, the birth place of immune cells and the site of immunological memory, various immune subsets find special microenvironments, which secure their survival and support their function. These special microenvironments are denoted as survival niches. The stable components of these niches are built by mesenchymal stromal cells (MSC), supported by the dense marrow vasculature [67,68]. Except for providing nutrients to diverse immune cell subsets, MSC may differentiate under certain conditions, e.g., during bone growth or bone regeneration after injury, into osteoblasts or chondrocytes. Under other conditions, e.g., in aged individuals or due to metabolic syndrome, they may differentiate into adipocytes [63]. Although these differentiation pathways are possible, the microenvironmental constraints which favor one or the other pathway are not yet fully understood. By using NAD(P)H-FLIM, others have provided evidence that cellular metabolism and differentiation stage of MSC are closely linked [21,42,47]. The interpretation of the NAD(P)H-FLIM data in these studies is hampered by either strongly simplified or too complex models of NAD(P)H fluorescence lifetime. If a monoexponential model is assumed, the decision whether a shorter mean fluorescence lifetime is caused by a lower NAD(P)H consumption (more free NAD(P)H) or by a higher activity of MDH or LDH cannot be made. Biexponential and multiexponential models, which facilitate this decision, are numerically instable and do not take into account the enzymatic landscape within the cell. The phasor approach proposed by Enrico Gratton and coworkers is a model-free analysis approach that circumvents the majority of numerical artifacts [69]. Nevertheless, data interpretation without a reference system of pure, relevant enzymes is also in this case not reliable.
To establish the systematic framework for a generally valid label-free NAD(P)H-FLIM data analysis, we first identified the most relevant, i.e., abundant NADH-and NADPH-dependent enzymes. These are responsible as catalyzers of biochemical reactions within the cell engine for energy supply, reductive biosynthesis, and other vital cell functions such as oxidative burst during phagocytosis. We measured the fluorescence decay of the coenzymes specifically bound to the respective enzymes, completing and validating published data from our and other laboratories. Although our reference NAD(P)H-dependent enzyme system for NAD(P)H-FLIM data is generally valid, an extension with other NAD(P)H-dependent enzymes relevant to specific cell populations or in certain tissue, organs, or pathologies is easily possible and desirable. Our systematic analysis framework of NAD(P)H-FLIM data relies on the following steps. (i) Calculating the SNR of the fluorescence intensity image; (ii) performing the phasor analysis of the time-domain NAD(P)H-FLIM data and identifying the pixel containing only free NAD(P)H; (iii) based on our reference enzyme system, calculating the enzymatic activity probability for all considered enzymes and assigning the predominantly activated enzyme to each pixel; and (iv) calculating the fraction of free to enzyme-bound NAD(P)H, as an index of metabolic activity.
We first validated our systematic analysis framework of NAD(P)H-FLIM data on homogeneous mixtures of NAD(P)H and pure enzymes. These data revealed the fact that the SNR value has a strong impact on the quality of the phasor data and their interpretation. If focusing only on the accuracy of the phase vector averaged over an image representing the central position of the phasor cloud, our results indicate that up to SNR 5, a linear correction of its modulation length, i.e., vector length, is necessary while the phasor angle does not need any correction. Above a SNR value of 5, both length and angle of the phase vector averaged over the whole image are numerically reliable and the effects on the position of the phasor cloud are caused by biologically relevant phenomena. Therefore, although it is often mentioned that FLIM is independent of intensity-time and not intensity is the imaging criterion-this holds true only for higher SNR values (above 5) for which the decay curve properties remain unchanged. Additionally, the SNR value in an image also strongly influences the width of the phasor cloud. The width of the phasor cloud has an immediate impact on the pixel-based calculation of enzymatic activity probabilities, enzyme activity assignment, and on the fraction free to enzyme-bound NAD(P)H. Our FLIM data on pure NADH solution and on NADH:MDH mixtures revealed that for high accuracy above 97%, SNR values above 20 are necessary. These results are in line with our finding concerning the interdependency of GFP fluorescence lifetimes on the SNR values [9]. Our data on the NADH:PDH mixture revealed that lower SNR values and high similarity of the phase vectors of several enzymes, in this case of PDH as compared to CTBP1, iNOS, and ADH, lead to a lower accuracy of enzyme assignment. Although the reference enzyme system cannot be changed, the SNR values in an image can be increased by increasing fluorophore concentration or the excitation power, smoothing the raw data in time or in space, or increasing the image acquisition time; whereas, in cells, a concentration increase is not possible, as this would change the cell physiology, all the other strategies to improve the SNR are linked to experimental drawbacks. An increase of excitation power may lead to photodamage [64], spatial smoothing of the data leads to loss of subcellular information, whereas temporal smoothing of the time-domain data decreases the accuracy of fluorescence decay curves and, thus, the accuracy of the phasor analysis. Increasing the acquisition time per image lowers the repetition rate of time-lapse imaging. In living cells and tissues, fast changes may be overseen in this way. Therefore, a combination of these procedures to improve SNR is required if FLIM is to be performed in living cells or tissues.
Taking this into account, we applied our systematic analysis framework to NAD(P)H-FLIM data on stromal-like 3T3-L1 cells, which develop into adipocytes upon stimulation. We decided to increase the SNR to highly reliable values by increasing the acquisition time of the imaged cells, because these cells do not migrate. In this way, we found shorter mean NAD(P)H fluorescence lifetimes in the nuclei as compared to the cytoplasm, in line with previous results in differentiating myoblast cells [19]. However, we found no differences in metabolism at the subcellular level, but a different predominant enzymatic activity in nuclei as compared to the surrounding cytoplasm, leading to this change in the NAD(P)H fluorescence lifetime. Although G6PDH was identified as the predominant active enzyme in the nuclei, the probabilities of enzymatic activity of LDH, SHD if bound to NADH, and GAPDH are at similar levels, due to the high similarity of the phase vectors of the enzymes. We expected to find CTBP1 as a NAD(P)H-dependent enzyme in the nucleus [57], but did not find any hints for the presence of LDH, G6PDH, GAPDH, and SDH there. In the cytoplasm, we found a high predominance of PDH, but also of iNOS and ADH. The contribution of CTBP1 (having a phase vector identical to that of PDH) is rather improbable, as this is mainly located in the cellular nucleus.
Our unexpected results regarding the metabolic activity in the nuclei of stromal-like cells led us to the conviction that our systematic analysis framework of NAD(P)H-FLIM retains the potential of detecting yet unknown enzymatic mechanisms related to cellular metabolism. This will have a tremendous impact on the way we will interpret the impact of cellular metabolism of immune and stromal cell populations, at a subcellular level and under in vivo conditions [61].

Two-Photon Microscope Setup Adequate for FLIM
Two-photon fluorescence imaging experiments were performed as previously described [33], using a specialized laser-scanning microscope based on a commercial scan head (TriMScope II, LaVision BioTec, Bielefeld, Germany). A near-infrared laser (Ti:Sa, Chameleon Ultra II, Coherent, Duisburg, Germany) tuned at 760 nm, repetition rate 80 MHz, and pulse width 140 fs was used as excitation source. The linearly polarized Ti:Sa beam was scanned over the sample by two galvanometric mirrors. A water-immersion objective lens (20×, NA 1.05, Apochromat, Olympus, Hamburg Germany) was used to focus the laser beam into the sample. The laser power was controlled by combinations of λ/2 wave-plates and polarizers. The ultrashort pulses of the laser were compressed using external compressor. NADH and NADPH fluorescence was collected in the backward direction using a dichroic mirror (775, Chroma, Marlborough, MA, USA), passed through an interference filter (466 ± 30 nm) and was detected by a GaAsP PMT (Hamamatsu, Herrsching, Germany) connected to previously described TCSPC electronics (LaVision BioTec). The TCSPC data were collected at a time resolution of 55 ps, over at least 9 ns and with a Gaussian-shaped instrument response function of 250 ps FWHM ( Figure S6). In all imaging experiments, we used an average maximum laser power of 10 mW to avoid photodamage. The acquisition time for an image with a field-of-view of 200 µm × 200 µm and a digital resolution of 512 × 512 pixel was 472 ms.

Phasor Analysis of Time-Domain NAD(P)H-FLIM Data
Fluorescence lifetime data were measured and analyzed as previously described [33,60,69]. The phasor approach transforms the time-domain data (the fluorescence decay curve), to a virtual, normalized phase domain by calculating the discrete Fourier transformation numerically (modulations frequency = 80 MHz). The transformation leads to a complex number the real and imaginary parts of which give the coordinates of the vector in the phase domain ("phasor"). That vector originated in (0|0) and points towards the half circle (centrum at (0.5|0), radius = 0.5), due to the exponential nature of the original time domain data. Because of the value normalization, the real part of the phasor reaches from 0 to 1, and those of the imaginary part from 0 to 0.5. In this way, short fluorescence lifetimes of homogeneous fluorophores (monoexponential decay) are located on the half circle at large real values, whereas with increasing lifetime, the real value decreases.